1. Introduction

Floral trait variations are widespread across flowering species, with distinct morphologies often adapted to specific environmental conditions. These phenotypic variations are sometimes associated with a change in the plant mating system (Gómez et al., 2008; Goodwillie et al., 2010; Toräng et al., 2017; Zeng et al., 2023). One of the most common evolutionary transitions is from outcrossing to self-fertilization (Stebbins, 1957; Lloyd, 1965; Husband and Barrett, 1993; Igic et al., 2008; Zhong et al., 2019). This transition is typically accompanied by a reduced floral display size, decreased pollen production, and altered stamen–pistil arrangements, a suite of traits known as the "selfing syndrome" (Sicard and Lenhard, 2011; Zhong et al., 2019; Yuan et al., 2025). Environmental gradients modulate the spatial distribution of mating systems and floral traits (Goodwillie et al., 2010; Abdusalam et al., 2022; Ma et al., 2023). For example, selfing populations have been observed at the peripheries of distributions, high latitudes, and high elevations and are characterized by pollinator scarcity and suboptimal conditions for outcrossing (Richards, 1997; Grossenbacher et al., 2015; Sorokhaibam et al., 2024). Under these conditions, selfing has evolved as a reproductive assurance mechanism, facilitating colonization and persistence in challenging environments (Baker, 1955; Stebbins, 1974; Jain, 1976; Morgan and Wilson, 2005).

Empirical studies have consistently revealed a correlation between reduced pollinator visitation and increased selfing, and outcrossing populations exhibit lower fruit and seed set under limited pollination (Lloyd and Schoen, 1992; Ashman et al., 2004; Kalisz et al., 2004; Christopher et al., 2021). However, it remains unresolved whether selfing lineages outperform outcrossing lineages in specific habitats, as observational data alone cannot separate the effects of intrinsic mating system traits from environmental filtering. Direct comparisons of fitness outcomes or phenotypes have rarely been evaluated between selfing and outcrossing populations under standardized conditions, largely due to the geographic segregation and local adaptation of these ecotypes (Hereford, 2010). Transplant experiments have emerged as an indispensable methodology for addressing this limitation. By reciprocally introducing individuals with divergent mating systems into native and novel ecological contexts, researchers can rigorously quantify the context-dependent mating success of different ecotypes (Boberg et al., 2014; reviewed by Johnson, 2025). For example, transplant experiments have shown that self-compatible morphs of Leavenworthia alabamica exhibit fitness advantages over their self-incompatible counterparts in both core and marginal habitats (Busch, 2005), whereas distinct ecotypes of Arabis alpina have shown adaptive divergence, outperforming others in their native environments (Toräng et al., 2015).

Mating portfolios quantify the relative contributions of selfing and outcrossing to the genetic diversity and fitness of offspring, complementing transplant experiments. Research on plant mating often neglects the importance of the paternal contribution and seldom examines male mate diversity (Tomaszewski et al., 2018; Christopher et al., 2019; Aljiboury and Friedman, 2022). Mating portfolios provide a comprehensive framework for evaluating the mating patterns of both individuals and populations. A mating portfolio encompasses the spectrum of mating outcomes, including selfing, outcrossing, and mate diversity, through both female and male functions (Barrett and Harder, 2017; Christopher et al., 2019; Yuan et al., 2023; Yang and Zhang, 2025). Determining the allelic contributions from both reproductive roles can provide a more nuanced understanding of the selective pressures shaping mating system evolution.

The evolution from heterostyly to homostyly is a classic transition from predominant outcrossing to selfing, as reviewed by Barrett (2019). Heterostyly is a genetically controlled floral polymorphism characterized in populations by two or three flower types (Darwin, 1877; Charlesworth and Charlesworth, 1979; Ganders, 1979; Barrett and Shore, 2008). These floral morphs commonly differ in style length and stamen height, promoting outcrossed pollen transfer. However, homostyles have been widely reported in many heterostylous families (Ganders, 1979; Barrett, 1992), which indicates that distyly is breaking down, leading to evolutionary transitions from predominant outcrossing to varying levels of selfing (Piper et al., 1984; de Vos et al., 2018; Yuan et al., 2023). Both distyly and homostyly can be found in some heterostylous species, providing opportunities to identify the ecological factors governing mating.

Primula oreodoxa is a well-studied heterostylous species in which both distylous and homostylous populations naturally occur (Yuan et al., 2017, 2019). Distylous and homostylous populations predominantly occupy low and high elevations, respectively, with a general demarcation at 1600–1700 m (Yuan et al., 2023). Previous studies have found that the origin of homostyly in P. oreodoxa occurred several times and that reproductive assurance is the primary driving force for the transition from distyly to homostyly (Yuan et al., 2017). Although mating system estimates based on the maternal mating success have confirmed the shift from outcrossing to selfing in natural populations (Yuan et al., 2023), little is known about the ecological and genetic performance of distylous and homostylous morphs in the same population across ecological gradients. This study addresses this gap through transplant experiments designed to directly compare the fitness of distylous and homostylous plants across different elevations.

In this study, we established six populations along an altitudinal gradient to investigate the fitness differences between distylous and homostylous Primula oreodoxa. Based on the natural distribution of P. oreodoxa, regions below 1700 m were designated as low-elevation habitats, and those above 1700 m were identified as high-elevation habitats. Both distylous and homostylous morphs were subsequently transplanted into these differentiated habitats to evaluate their adaptability and reproductive success in diverse environments. We estimated the mating portfolios for each population in these two habitats, considering both maternal and paternal perspectives. Our key research questions were as follows: (1) How do fruit set and the number of seeds vary between distylous and homostylous flowers across low- and high-elevation habitats? We hypothesized that the distylous morph would have lower fertility than the homostylous morph at high elevations because homostyles are capable of autonomous self-pollination, which provides reproductive assurance under pollinator-limited conditions (Yuan et al., 2017). At low elevations, we predicted that homostyles would exhibit comparable reproductive success relative to distylous morphs, as low-elevation populations likely do not experience pollinator limitations. (2) How does the selfing rate vary among floral morphs, populations, and habitats? We hypothesized that homostylous flowers would have a higher selfing rate than distylous flowers and that the average selfing rate would increase from low to high elevations, consistent with increasing pollinator limitations. (3) How do mate diversity and assortative mating patterns differ across populations at low- and high-elevation habitats? (4) Are there correlations among the population selfing rate, mating diversity, and pollinator visitation frequency? Through manipulative experiments, we aimed to determine how pollinator availability across ecological gradients shapes the mating characteristics and divergence of mating systems in P. oreodoxa.

2. Materials and methods 2.1. Study species

Primula oreodoxa is a perennial herbaceous plant native to western Sichuan Province, China (102°–104°E, 28°–31°N). It blooms from March to May, producing pink flowers, and its fruit ripens from May to June. Previous research has shown that both the distylous and homostylous morphs of this plant are fully self-compatible (Yuan et al., 2017, 2019). However, only the homostylous morph is capable of autonomous self-fertilization, and the distylous morph requires pollinators for successful reproduction.

2.2. Transplant experiment design

Transplant experiments were conducted along an elevation gradient, with six experimental populations of Primula oreodoxa established at different altitudes: MTG (1269 m), ZLP (1623 m), JCC (1668 m), CSQ (1702 m), XXC (2099 m), and LDP (2420 m). Based on previous field observations showing that distylous populations are generally found below 1700 m and that homostylous populations are found above this elevation, we classified the first three transplanted populations as the low-elevation habitat group and the last three as the high-elevation habitat group for subsequent data analysis (Fig. S1).

A total of 15–48 mature plants of both distylous and homostylous morphs were randomly transplanted into each population from glasshouses at the Biological Resources Research Station at E'mei Mountain (Sichuan, China). The mature plants in the greenhouse were the descendants of two natural distylous (WWS and DWS), one natural homostylous (TTS), and two natural mixed populations (HZG and QLP). We included the original maternal plants as a covariate in our statistical analysis to minimize the influence of potential genetic background differences among transplanted individuals. Among the transplanted populations, JCC had locally occurring distylous plants, and CSQ, XXC, and LDP had locally occurring homostylous plants. We dug up and replanted the local plants in their original location to minimize the transplant effect. Except for watering at planting, we did not perform any human intervention. The transplant design aimed to determine how ecological variation (pollination environment) across elevations affects reproductive outcomes.

Mature plants of both distylous and homostylous morphs were transplanted into six experimental populations during the autumn of 2017. Due to heavy damage to MTG and ZLP populations before the fruit harvest in 2018, we re-transplanted some mature plants of the same origin as these two populations during the autumn of 2018. Our estimations of mating patterns were conducted in 2020, except for population JCC, in which we estimated the mating patterns in 2018 (Table S1). In 2020, we marked all flowering plants in each transplanted population and conducted experiments to compare the performance of distylous and homostylous plants within and among groups. Because the transplanted distylous and homostylous plants belonged to different populations, their flowering periods differed, potentially affecting the mating pattern phenotypes observed in transplant gardens. Although we did not record the flowering periods of individual plants in all gardens, the flowers of distylous plants in MTG, ZLP, XXC, and LDP gardens opened about one week earlier than those of homostylous plants. Therefore, when estimating the mating patterns, we used plants and flowers of the two floral morphs with overlapping flowering periods. In the CSQ and JCC populations, the distylous and homostylous flowers bloomed almost simultaneously.

2.3. Pollinator visitation

In each transplanted population, insect visits to flowers were documented during the flowering season (2018–2020). To accomplish this, a 5 × 2 m plot was identified during peak flowering within each population. We recorded all pollinators that visited flowers on sunny days between 9.00 and 17.00 h (referred to as a "day"). All insects observed during this study have been reported in our previous studies (Yuan et al., 2017, 2023). Pollinators were divided into two categories, namely long- and short-tongued insects, based on their proboscis length, following the classification of Yuan et al. (2017). For each day of observation, we recorded the number of visits by both pollinator types and the total number of blooming flowers. We observed 2–8 sunny days for each transplanted population. Following the approach of Yuan et al. (2023), we constructed generalized linear mixed models (GLMMs; Stroup, 2013) using the ‘glmmTMB’ (v.1.1.4; Brooks et al., 2017) and ‘emmeans’ packages (v.1.7.5) in R (v.4.2.0; R Core Team, 2021) to investigate the relationship between pollinator visits and population elevation. The number of pollinator visits was modeled with a log link function using a zero-inflated negative binomial distribution. We considered pollinators as fixed effects and populations as random effects. Sampling year was used as a categorical factor, and the population elevation and the ln (number of flowers per observation plot) were considered continuous independent variables. We integrated pollinator visitation over multiple years, including that previously published by Yuan et al. (2023), to enhance statistical robustness.

2.4. Fruit and seed set

During the flowering seasons from 2018 to 2020, we identified flowering plants in each of our transplanted populations and compared the fruit and seed set among floral morphs, populations, and habitats. We randomly marked the inflorescences of each plant and counted the number of flowers and buds. After 6–8 weeks, we counted the number of fruits on the marked inflorescences and calculated the fruit set ratio. We also randomly selected two or three fruits from each plant and counted the number of seeds in each fruit. We estimated the fruit set ratio for 653 individuals, with 121, 174, and 358 from L, S, and homostylous plants, respectively. Additionally, we sampled 98/253, 140/271, and 199/376 plants/fruits from the L, S, and homostylous plants, respectively, to count the number of seeds. The summary of fruit and seed set is provided in Table S2.

To investigate the effects of floral morphs, populations, and their interaction on fruit and seed set, we constructed GLMMs using the ‘glmmTMB’ and ‘emmeans’ packages in R. For fruit set, we assumed a Gaussian distribution, and for seed set, we used a Conway–Maxwell–Poisson distribution to account for overdispersion. In both models, population, floral morph, their interaction, and year were treated as fixed effects. Estimated marginal means (±SE) were calculated using the ‘emmeans’ package for post hoc comparisons.

2.5. Sampling, DNA extraction, and genotyping

All flowering plants were marked in each transplanted population, and fresh leaf material was collected and dried in silica gel for subsequent genotyping. Mature fruit was harvested from late May to late June 2020, and seeds from two fruits per plant were mixed and sown on soil-filled rectangular flats. The seedlings were cultivated in a glasshouse at the Biological Resources Research Station at E'mei Mountain. Based on the number of seedlings emerged per number of sown seeds, the average germination rate was 40%. From October to December 2020, we randomly selected 8–12 seedlings (averaged 11.5 ± 0.1) from each maternal family and dried their fresh leaves in silica gel. The sample sizes for seed families and progeny per floral morph for each transplanted population are provided in Table S1. We extracted the total genomic DNA from dried leaf tissues of parental plants and progeny using a modified cetyl trimethyl ammonium bromide protocol (Doyle, 1991) and used multiplex PCR with 11 simple sequence repeat (SSR) markers for DNA amplification. Among the 11 SSR markers, eight were selected from previously developed SSR markers for Primula oreodoxa (Yuan et al., 2018), and the other three SSR markers were first reported by Yuan et al. (2023). The SSR markers were divided into three groups, with each of the four primer pairs co-amplified in a single PCR tube. The amplified 10 μL of mixture for SSRs comprised 5 μL of 2 × Multiplex PCR Mastermix (Mei5bio, Beijing, China), 1 μL of primers (3 or 4 primers were mixed in advance), 2 μL of deionized water, and at least 2 μL of genomic DNA. We then ran the PCR using the following cycling conditions: initial denaturation at 95 ℃ for 10 min; 8 cycles of 95 ℃ for 30 s, 60 ℃ for 30 s with an increment of −1 ℃ per cycle, and 72 ℃ for 30 s; 32 cycles of 95 ℃ for 30 s, 52 ℃ for 30 s, and 72 ℃ for 30 s; and ending with extension at 72 ℃ for 10 min. We scanned PCR products using an ABI PRISM 3100 Genetic Analyzer (Invitrogen) with internal size standard GeneScan^TM 500 LIZ and conducted allele binning and calling using GeneMarker v.3.0.0 (SoftGenetics LLC, State College, PA, USA).

2.6. Genetic diversity and paternity analysis

Statistical software GenAlEx v.6.5 was used to calculate the genetic parameters, such as the number of observed alleles (N_A), expected heterozygosity (H_E), observed heterozygosity (H_O), Shannon's information index (I), and fixation index (F). The fixation index was calculated using the following formula: F = 1−H_O/H_E. These parameters were estimated separately for each locus, each floral morph, maternal plants, and progeny. Generalized linear models were then used to identify differences in these parameters between floral morphs and habitats for progeny. Floral morphs were considered fixed effects, and populations and loci were treated as random effects. The estimated genetic diversity of maternal plants was also included as a covariate variable to control for the potential influence of maternal genetic background on the offspring. We hypothesized that the progeny of homostylous plants would have a lower N_A, H_O, and I and a higher F than the progeny of distylous plants, as homostylous plants are expected to have a higher selfing rate than distylous plants.

We also conducted paternity analysis for the six transplanted populations using COLONY 2.0.6.5, employing maximum likelihood to identify the putative parents for each offspring within each transplanted population. In the software, we used the following settings as recommended by Wang et al. (2012): polygamy for both female and male parents, with inbreeding, no clonality, hermaphrodite, diploid, full-likelihood method, medium length run, medium precision, and allele frequencies not updated. Locus-specific error rates were calculated for each population, and the average genotype error rate was set to 0.01. Progeny from the same maternal plants was coded as known maternal sibships with the paternal sibships unknown.

2.7. Mating portfolio estimation

We estimated the mating portfolios using a set of parameters, including the proportion of selfing, outcrossing number of mates, mate diversity, and intra-versus intermorph mating frequency. All calculations were based on progeny with a parental assignment probability of ≥ 95% and obtained from COLONY analysis. Specifically, the selfing rate was defined as the proportion of progeny assigned to the maternal plant. The number of mates was calculated as the number of distinct outcrossing partners (fathers or mothers) per plant. Mate diversity was estimated using Hill's diversity indices (Hill number 1D, reflecting the effective number of common mates; see Supplemental Methods) based on the mate identity and frequency. Inter- and intramorph mating were determined by comparing the floral morphs of the assigned maternal and paternal plants. Considering the relative morph frequency in each transplanted population, we calculated assortative versus disassortative mating following the method described by Yuan et al. (2023). We set the probability threshold at 95% to predict whether a genotyped seed was produced by selfing and could be assigned a known father. Statistical analysis was performed using subsets of the 2932 genotyped seeds. We used generalized linear mixed-effects models to determine the mating portfolios for maternal and paternal families, with the selfing rate, inter-morph mating frequency, number of mates, and ¹D index considered response variables. Fixed factors included the floral morphs, habitats, populations nested in habitats, and their interactions, with the families in each population considered random variables. Covariate variables were included when analyzing the number of mates, using log-transformed total outcrossed seeds produced by outcrossed mates. We conducted GLMM analyses using the same R packages described in Section 2.3, selecting appropriate distributions and link functions as the dependent variables.

We determined whether pollinator availability influences mating patterns differently among floral morphs. To assess this, we used Pearson correlation in R (v.4.1.1) to identify relationships between visitation frequency of long- and short-tongued pollinators and mating parameters (i.e., selfing rate and mating diversity) separately for each morph. The visitation frequency data were obtained from the adjusted average visits per flower per day, and the selfing rate and mating diversity data were obtained from the adjusted selfing rate and number of mates.

3. Results 3.1. Pollinator visitation rates

The six transplanted Primula oreodoxa populations were visited by diverse pollinators, including long-tongued bees, bee flies, hawkmoths, butterflies, small bees, and syrphid flies (Fig. S2). As reported in a previous study (Yuan et al., 2017), the first four groups were long-tongued pollinators, with long-tongued bees being the most abundant, and the last two groups were short-tongued pollinators. Using all pollinator visit data, we observed a decrease in the adjusted visits of long-tongued pollinators with increasing population elevation (Fig. 1, Z = −6.45, P < 0.0001). However, we did not find a significant correlation between the visits of short-tongued pollinators and population elevation (Fig. 1, Z = −0.20, P = 0.8407). The adjusted number of long-tongued visitors was positively correlated with the number of flowers per plot (Fig. S3, Z = 7.76, P < 0.0001), but the adjusted number of short-tongued pollinators was not correlated with the number of flowers per plot (Fig. S3, Z = 0.91, P = 0.36).

3.2. Fruit and seed set

Across all six transplanted populations, each inflorescence (n = 653) produced 6.9 ± 0.1 flowers, on average, with 6.4 ± 0.1 setting fruit. The fruit set ratio per inflorescence ranged from 0.29 to 1.000 (mean ± SE: 0.93 ± 0.01). The GLMM showed that the fruit set ratio differed among populations (Table S3, X₅² = 24.452, P < 0.001), whereas no significant difference was observed in the fruit set ratio among the three floral morphs in the low-elevation habitat (Fig. 2A, P > 0.05). In the high-elevation habitat, the fruit set ratio of the S morph was significantly lower than that of the H morph (Fig. 2A, P < 0.01); this pattern was particularly evident in populations XXC and LDP (Fig. S4A). However, no significant difference was observed between the L and S morphs or between the L and H morphs (L vs. S, P = 0.298; L vs. H, P = 0.712; Fig. 2A). Furthermore, the fruit set ratio of the S morph was significantly reduced in high-elevation habitats compared with low-elevation habitats (P < 0.001, Fig. 2A). Significant differences were observed in seed set among populations (X₅² = 19.023, P < 0.01) and floral morphs (H morph: 304 ± 5.5; L morph: 280.2 ± 5.6; S morph: 258.4 ± 5.6; X₂² = 41.248, P < 0.001, Table S3). In low-elevation habitats, seed set did not significantly differ among populations or floral morphs (Fig. 2B, all P > 0.05). However, in high-elevation habitats, the number of seeds per carpel was significantly lower in the L and S morphs than in the homostylous morph (Fig. 2B, all P < 0.001), especially for populations CSQ and XXC (Fig. S4B, P < 0.01).

3.3. Genetic diversity and paternity assignment

We analyzed 2932 offspring and amplified 6–28 alleles for 11 SSR markers, with an average of 16.6 ± 2. Populations XXC and LDP, which were at higher elevations, exhibited the lowest H_O, H_E, and I values (Table S4). The F values were greater than zero for all populations (Table S4), indicating inbreeding in all of these populations. In low-elevation habitats, there were no differences in these genetic parameters among the three floral morphs (Fig. 3). However, in high-elevation habitats, H_O, H_E, and I were significantly lower in the H morph than in the L and S morphs (Fig. 3). In progeny produced by S morphs, the genetic diversity decreased with the increase in elevation (Fig. 3).

We sampled 650 mature individuals, 256 maternal families, and 2932 offspring from 6 transplanted populations, with an average of 11.5 ± 0.1 progeny per maternal family (range: 4–15, with most families having 12 offspring). Paternity analysis assigned 2125 progeny a father with over 95% probability. Among the populations, the proportion of successfully assigned progeny varied from 51.20 to 90.30%, with an average of 71.85 ± 5.52% (Table S5).

3.4. Selfing rates

To determine the proportion of self-fertilized and outcross-fertilized seed, paternity analysis was conducted, differentiating between female and male contributions. Generally, the number of outcrossed seeds produced by female parents was significantly higher than that produced by male parents (Fig. S5A). Both female and male selfing rates significantly varied across habitats, populations, and floral morphs both within and between populations and habitats (Table S6). In low-elevation habitats, populations MTG and JCC showed higher selfing rates than population ZLP (Fig. 4A, P < 0.01). However, in high-elevation habitats, populations XXC and LDP showed higher selfing rates than population CSQ (Fig. 4A, P < 0.01). In low-elevation habitats, the selfing rates of the L, S, and H morphs were 0.107 ± 0.031, 0.139 ± 0.039, and 0.414 ± 0.081, respectively, with the H morph having a significantly higher selfing rate than the L or S morph (Fig. 4B, P < 0.01). The selfing rates of dimorphic morphs significantly increased from lower elevations to higher elevations (all P < 0.01). However, in high-elevation habitats, there were no differences in selfing rates among the three floral morphs (Fig. 4B, P > 0.05). The pattern observed for male selfing rates (Fig. 4C and D; Table S6) was consistent with that for female selfing rates, and a positive correlation was observed between female and male selfing rates at the population level (Fig. S5B). The female selfing rate was negatively correlated with long-tongued pollinator visitation in the L morph (Fig. 5A, P < 0.05), and the male selfing rate showed a similar negative correlation with long-tongued visitation in the S morph (Fig. 5E, P < 0.05). No significant correlations were found in the H morph for either sex function (Fig. 5C and F). Visitation by short-tongued pollinators showed no significant associations with selfing rates in the L or S morph but was positively correlated with the male selfing rate in the H morph (Fig. S6F, P < 0.05).

3.5. Mate diversity

In low-elevation habitats, there were no significant differences in the number of female mates among the three populations (Fig. 6A). However, in high-elevation habitats, population CSQ had a higher number of male mates than populations XXC and LDP (Fig. 6A, P < 0.01). The number of male mates was generally consistent with the number of female mates (Fig. 6B). Moreover, there was a positive correlation between the number of male and female mates (Fig. 6C). The number of male and female mates was negatively correlated with female (Fig. 6D) and male selfing, respectively (Fig. 6E). The patterns of the ¹D index for female and male parents were generally consistent with the number of mates (Fig. S7A and B). We found that the visitation frequency of long-tongued pollinators was positively correlated with the number of both female and male mates in the L and S morphs, but no significant associations were detected in the H morph for either sex function (Fig. 5G–L). In contrast, visitation by short-tongued pollinators was only negatively correlated with the number of female mates in the H morph (P < 0.05; Fig. S6I), with no significant relationship observed in the other morphs (Fig. S6G, H and J–L).

3.6. Inter- and intra-morph mating

The frequency of female inter-morph mating was significantly higher in low-elevation habitats than in high-elevation habitats (Fig. S8A, 0.755 ± 0.042 vs. 0.388 ± 0.101, P < 0.01). In low-elevation habitats, population ZLP had a significantly higher frequency of inter-morph mating than population JCC (Fig. S8A, population ZLP, 0.894 ± 0.038, population JCC, 0.627 ± 0.082, P < 0.01). Populations MTG and ZLP mainly exhibited male inter-morph mating, whereas other populations mainly displayed random mating or intra-morph mating, especially in high-elevation habitats (Fig. S8B). Disassortative mating was mainly observed in the three populations in low-elevation habitats, whereas the populations in high-elevation habitats did not deviate from random mating. Considering the relative frequency of other floral morphs in the populations to estimate the proportion of female inter-morph mating, we found that most populations exhibited random mating for all three floral morphs, especially in the S morph and homostylous plants (Fig. 7).

4. Discussion

Our study supports that pollinator availability influences mating system transitions in Primula oreodoxa, with homostylous individuals having a reproductive advantage in pollinator-limited environments. Through transplant experiments, we demonstrated that population selfing rates and the number of mates were negatively and positively correlated with pollinator visitations, respectively, along an elevation gradient (Fig. 5). Additionally, male and female parents functioned consistently within populations (Fig. 4, Fig. 6). The prevalence of random mating in high-elevation populations (Fig. 7) signals the breakdown of the distyly polymorphism. Thus, in P. oreodoxa, as the elevation increases, the distyly polymorphism may be disrupted by the fixation of homostyly and a shift from predominant outcrossing to selfing. Moreover, our findings highlight the importance of considering both maternal and paternal contributions when assessing mating diversity to provide a more comprehensive perspective on mating system evolution.

4.1. Fecundity differs among floral morphs and habitats

In our experiments, fruit and seed set analyses showed that distylous and homostylous plants had similar reproductive performance in low-elevation habitats (Fig. 2). Specifically, the comparable fruit and seed production among floral morphs in these regions suggest that the reproductive capacities of distylous and homostylous plants converge under unlimited pollinator visitation. In contrast, in high-elevation habitats, the fruit set ratio of the S morph was significantly reduced compared with that of homostylous plants (Fig. 2A). This differential performance indicates that the ability of homostylous individuals to self-fertilize confers a reproductive advantage. Our findings thus align with the classical reproductive assurance hypothesis (Baker, 1955) and are supported by previous studies demonstrating how autonomous self-pollination serves as an adaptive mechanism to overcome pollination deficits in non-heterostylous (Kalisz et al., 2004; Moeller, 2006) and heterostylous species (Piper et al., 1986; Carlson et al., 2008; de Vos et al., 2018; Jiang et al., 2018). Although both floral types are self-compatible in P. oreodoxa, controlled experiments have shown that only homostylous morphs have effective autonomous self-pollination (Yuan et al., 2017).

The factors causing the observed reduced fruit set ratio in short-styled morphs are probably associated, in part, with the relative frequency of short-versus long-tongued pollinators visiting flowers, their morph-specific foraging behavior, and differing influences on pollen transport (see Keller et al., 2014; Yuan et al., 2017). As elevation increased, the visiting frequency of long-tongued pollinators decreased significantly (Fig. 1). The loss of "high-quality" long-tongued pollinators may have reduced pollen transfer in S morph plants. The seed number per fruit of homostyles was significantly higher than that of dimorphic plants (Fig. 2B), which can be explained as follows: First, limited pollinator activity at higher elevations restricted the seed set of dimorphic plants, and homostyles can guarantee seed set through autonomous self-pollination. Second, pollen transfer is not restricted in high-elevation habitats, and the higher number of seeds per fruit in homostyles results from the increased number of ovules in these plants, although previous ovule counts have not shown significant differences among floral morphs (Yuan et al., 2025).

In addition to pollinator limitations, other environmental constraints, particularly resource limitations, such as soil moisture and nutrient availability, may contribute to the observed variation in fecundity among floral morphs and habitats. High-elevation environments are often characterized by water and nutrient limitations, influencing plant reproductive success both directly, through physiological stress, and indirectly, by altering plant–pollinator interactions. Recent studies (Wu et al., 2023; Cha et al., 2025) have shown that these abiotic factors significantly reduce reproductive output by influencing the expression of floral traits and modifying the behavior of floral visitors.

4.2. Varied selfing rates and mate diversity

We found that both the female and male selfing rates varied significantly among habitats, populations, and floral morphs in the six transplanted populations (Table S6 and Fig. 4). In low-elevation habitats, the female and male selfing rates of homostyles were significantly higher than those of dimorphic morphs (Fig. 4B–D), which was consistent with our expectation that the dimorphic morphs mainly have outcrossing and that homostyles have a higher frequency of self-fertilization (Yuan et al., 2023). Differences in selfing rates between dimorphic morphs and homostyles are largely determined by the proximity of the stigma and anthers, as selfing increases with the decrease in the distance between them (Karron et al., 1997; Takebayashi et al., 2006; Opedal, 2018; Ma et al., 2021). Population CSQ in the high-elevation habitat had a significantly lower selfing rate than population JCC in the low-elevation habitat. A small degree of herkogamy was preliminarily observed in the homostylous plants in population CSQ, although it was not formally quantified. This may have contributed to the increased proportion of outcrossing in this population (de Vos et al., 2012, 2018).

Both the female and male selfing rates increased as the elevation increased (Fig. 4A–C), suggesting that selfing rates are strongly influenced by the environmental factors of the population in plants with mixed mating systems (Kalisz et al., 2004). High selfing rates with reduced pollinator visitation may imply a reduction in the investment in floral morphology and pollinator attraction (Goodwillie et al., 2010; Zeng et al., 2023). The balance between selfing and outcrossing influences genetic diversity and the population structure, which, in turn, affects the mating opportunities for individuals and the mating system within populations (Wright et al., 2008, 2013; St Onge et al., 2011). Generally, the mating consequences of homostyly differ from those of distylous morphs, as the genetic diversity of the progeny produced by homostyles was significantly lower and the fixation index was higher than those of the L and S morphs (Fig. 3). Our study estimated the male selfing rate among populations and found that maternal and paternal selfing rates were positively correlated, with paternal selfing rates slightly higher than maternal selfing rates (Fig. S5).

Few studies have examined the component of mate diversity arising from a plant's paternal success in siring seeds on other plants (see Christopher et al., 2019). In this study, we determined the number of mates that an individual had through female and male functions and calculated the ¹D index, providing additional insights into hermaphrodite mating success. Male fitness is influenced by the amount of self-fertilization that occurs, as it affects pollen competition and the opportunity to sire offspring (Bernasconi, 2003; Karron et al., 2006). Notably, there were no differences in the number of male or female mates among floral morphs in both low- and high-elevation habitats (Table S6), suggesting that the plant's mating success as a female and male parent is relatively consistent among floral morphs. Analysis of the ¹D index yielded similar results. As it considered both maternal and paternal contributions, our paternity analysis showed a strong positive correlation between female and male selfing rates across populations. This parallel pattern reinforces that P. oreodoxa individuals function consistently as both female and male, further emphasizing that environmental constraints drive mating system evolution in this species.

4.3. Divergence of mating systems and associations with pollinator visitation

The heterostyly polymorphism facilitates inter-morph pollen transfer (disassortative mating) in plants (as reviewed in Ganders, 1979; Lloyd and Webb, 1992; Barrett and Shore, 2008; Keller et al., 2014; Zhou et al., 2015). In our study, none of the six transplanted populations significantly deviated from random mating at the population level, but we found a varying degree of assortative mating among floral morphs within most populations (Fig. 7). The three populations in low-elevation habitats did not significantly deviate from disassortative mating, with population JCC exhibiting the lowest proportion of inter-morph mating. In contrast, at high elevations, random mating was dominant in all populations (Fig. S8). Furthermore, the proportion of disassortative mating of dimorphic morphs was significantly lower in high-elevation habitats than in low-elevation habitats (P < 0.001). Considering the relative frequency of morphs in each transplanted population, the L morph plants achieved mainly inter-morph mating in most populations; however, the S morph and homostylous plants did not deviate from random mating in most populations. These results suggest that maintaining the distyly polymorphism through sufficient disassortative mating was difficult in these transplanted populations, especially in populations at higher elevations.

The observed increase in the selfing rate along the elevation gradient and divergent fertility among floral types are linked to the frequency of flower visitation by pollinators, particularly long-tongued pollinators. We found that long-tongued pollinator visitation was significantly negatively correlated with the female selfing rate in the L morph and the male selfing rate in the S morph (Fig. 5A–E), supporting context-dependent mating system variations (Kalisz et al., 2004; Brunet and Sweet, 2006). Our study provides the first evidence to date of the negative correlations between pollinator visitations and both female and male selfing. Long-tongued pollinator visitation was also positively correlated with the number of mates in the L and S morphs (Fig. 5G, H, J), underscoring their critical role in enhancing outcrossing opportunities. This aligns with the spatial separation (herkogamy) between reproductive organs in these morphs, where long-tongued pollinators are more efficient at mediating pollen transfer between plants compared to short-tongued pollinators (Pérez-Barrales and Arroyo, 2010). In contrast, short-tongued pollinator visitation showed no significant relationship with selfing rates or the number of mates in the L and S morphs (Fig. S6). However, for the homostylous morph, where stigmas and anthers are at the same height, short-tongued pollinator visitation exhibited a positive correlation with selfing rate and a negative correlation with the number of mates (Fig. S6). The underlying reason for this pattern remains unclear, but it may stem from the floral morphology of homostyles. The foraging behavior of short-tongued pollinators, particularly their tendency to linger near the corolla opening, could facilitate self-pollination in these flowers. Nevertheless, further evidence is required to confirm this hypothesis.

In Primula oreodoxa, high-altitude environments consistently limit pollinator availability, conferring a reproductive advantage to homostylous over distylous plants, especially short-styled morphs. Despite the reduced pollinator visitation at high-elevations, our results indicate that the H morphs achieved comparable outcrossing siring success to the L and S morphs. The number of mates and mate diversity were not significantly different among morphs (Table S6 and Fig. S7). This pattern implies that the H morphs rely on occasional but successful outcrossing events and compensate for limited pollinator services through autonomous selfing to ensure seed production. Additionally, homostyle individuals outcrossed to some extent and replaced short-styled individuals in mating with long-styled individuals. This is because the anthers of homostyles and the S morph are positioned at the same height, allowing feasible pollen transfer from homostyles to long-styled individuals. Homostylous plants not only facilitate pollen transfer to the L morph but also ensure seed production through self-fertilization. Consequently, short-styled plants may be replaced by homostylous plants over time, which has been documented in other Primula species (Crosby, 1949; Curtis and Curtis, 1985; Mora-Carrera et al., 2023). This potential shift in floral morph frequency has profound implications for the genetic structure and evolutionary trajectory of P. oreodoxa populations.

4.4. Limitations and future directions

Several limitations should be acknowledged despite the robustness of our findings. First, the study's temporal scope was limited. Many studies examining pollination and fruit set in selfing versus outcrossing populations have been confined to a single year. Although our research compared fruit and seed set from 2018 to 2020, data were not recorded for all populations in all three years. Long-term comparisons are essential to reliably observe changes in the proportions of floral morphs within populations as well as variations in fruit set and seed production over time (reviewed by Johnson, 2025). Second, the evolution of selfing often results in inbreeding depression (Lande and Schemske, 1985; Charlesworth and Willis, 2009), with the advantages of selfing over outcrossing depending on the balance between inbreeding depression and reproductive assurance (Kalisz et al., 2004; Delgado-Dávila and Martén-Rodríguez, 2021). Although we did not conduct extensive seed germination tests, transplant experiments and the survival rates of homostylous individuals showed that they were well suited to high-elevation habitats and likely did not experience severe inbreeding depression. To address this, subsequent work should incorporate detailed fitness assessments, such as controlled germination trials, to quantify potential depression more precisely. Finally, some transplanted populations still contained native individuals. Although we removed and replanted native individuals to eliminate transplant effects, they may still have been better adapted to the current elevation than other individuals, which is a common issue in transplant experiments (reviewed by Johnson, 2025). Our results showed that homostylous plants had a higher survival rate at higher elevations, whereas distylous individuals had a higher survival rate at lower elevations, suggesting adaptive floral morph transitions. However, to mitigate such effects, future experiments should control for native presence more rigorously, for example using isolated plots or genetic tracking, and account for interannual variations in fecundity and mating patterns. Overall, these limitations highlight opportunities to refine our understanding of evolutionary adaptations in dynamic environments.

5. Conclusion

Our transplant experiments in Primula oreodoxa showed that homostyles and distyles had comparable performance at lower elevations, but homostyles demonstrated advantages at higher elevations. This suggests an evolutionary trend toward selfing in high-altitude environments. The variation in the selfing rate and number of mates among different morphs can serve as a model for understanding mating system transitions. Our previous studies on the evolutionary breakdown of distyly in P. oreodoxa identified populations exclusively composed of selfing homostyles, which may have originated from mixed populations (Yuan et al., 2017). Using transplant experiments and mating parameters, the current study suggests that homostyles achieve a similar mate diversity as distylous morphs. Together with the higher fruit set of homostyles in high-elevation habitats, this could facilitate their spread and replacement of the dimorphic morphs, particularly short-styled plants, in these environments. Furthermore, our study highlights the importance of considering both female and male reproductive success when examining mating patterns. Overall, this research contributes to our understanding of mating system evolution in flowering plants and may have broader implications for plant reproductive ecology.

Acknowledgement

We thank Mingsong Wu, Tong Zeng, Xiaojuan Li and Cehong Li for assistance with field studies. We thank Mrs. Yunxiao Liu for the drawings of Mount E'mei. We also extend our gratitude to Professor Spencer C.H. Barrett for his insightful comments and valuable suggestions regarding the manuscript. This research was funded by grants from the National Natural Science Foundation of China (31800314, 32370239, U160323); the foundation of South China Botanical Garden, Chinese Academy of Sciences (QNXM-06) to SY and the Doctoral Research Foundation of China West Normal University (412994).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2025.08.005.