Primula wannanensis usually grows on shaded slopes under or at the edge of deciduous, broad-leaved forests between 50 and 1,000 m a.s.l. (He et al., 2021). Plants mostly possess 10–30 pinnate compound leaves and 3–10 (up to 15) scapes, each bearing 6–9 flowers in 2–3 layers. The flowering period is long, with inflorescences blooming in succession from late February to mid-May. Flowers are pollinated primarily by the bee species, Bombylius major and Anastoechus chinensis (Shao et al., 2008). Controlled pollination indicates that the expression of self-incompatibility is not strict and varies extensively among individuals and populations; all distylous populations contain self- or partly self-compatible individuals, although their proportions are significantly higher in peripheral populations than in central populations (Shao et al., 2012, 2019). Primula wannanensis has two distinct populations, distylous and homostylous, based on the morph structure. Distylous populations contain almost equal proportions of L- and S-morphs, and are mainly distributed at (or near) the geographical distribution center (e.g., in and around the Huangshan mountain) and in marginal areas (e.g., near the Qidu township). Homostylous populations have monomorphic long or short styles (with stamens and stigma both at the mouth or middle of the corolla tube) and occur only in the most marginal areas of the distribution range (He et al., 2021). Four distylous populations of P. wannanensis were used in this study (Supplementary Table 1): two populations in the central distribution area (GCN and HDY), and two populations in the marginal distribution area (QDX and YLD). The GCN and HDY populations were found growing under temperate conditions in the mixed deciduous forest in the Huangshan mountain, whereas QDX and YLD populations grew along the roadside or near the edge of the forest on small hills near the Qidu township.

In mid-early March and mid-late April 2017, approximately 150 newly opened flowers were sampled from each population, with at least 0.5 m distance between two flowers. Because of differences in population area and plant density, the number of samples collected from each population varied from 61 to 365. In total, 1,242 flowers were sampled including 797 early opening flowers and 445 late-opening flowers. To minimize the influence of inflorescence structure on flower morphology, all flowers were picked from the first layer. Each flower was sliced open longitudinally and digitally photographed. The following six floral traits were measured or calculated: stigma height (SH), anther height (AH), corolla tube length (CTL), herkogamy (HE, i.e., stigma–anther separation), relative herkogamy (RHE; the rate of HE/CTL reflect whether the change in HE is proportional to the overall change in CTL), and corolla diameter (CD) (Supplementary Figure 1). All traits were measured to an accuracy of 0.01 mm using the ImageJ software (Rueden et al., 2017).

Before the start of florescence (early March) in 2018, approximately 60 plants were marked in each population, and covered by a dense net (mesh diameter < 0.5 mm) to exclude pollinators (Shao et al., 2008). Based on the preliminary observation that the degree of herkogamy does not vary among flowers within a layer and does not change after flower opening (unpublished data), one flower from the first layer was picked from each marked plant at the early florescence stage in mid–early March for floral trait measurement, while another flower in the same layer was marked for testing the autonomous selfing capacity. One month later, late-opening flowers were observed and marked as described above. Floral traits were measured according to the above description. Fruits of marked flowers were collected in early and late May, respectively, and the number of seeds and undeveloped ovules per fruit was counted under a stereomicroscope (Shao et al., 2019). Additionally, to evaluate the seed setting ability of plants under natural conditions (with the participation of pollinators), approximately 60 plants were marked in each population, and their floral traits and seed setting ability were investigated in the same manner as that for plants covered with a net.

To obtain a general overview of the spatial and temporal variation in floral traits, principal component analysis (PCA) of the six floral traits was performed separately for each morph type at different sites and different flowering phases. To test whether the overall floral traits of distylous flowers differed between the central and marginal sites and between early and late flowering phases, multivariate analysis of variance (MANOVA) was performed, with the site, flowering phase, and site × flowering phase interaction as fixed factors. Subsequently, univariate analyses of variance were performed to identify the floral traits that differed significantly between the two sites and the two flowering phases. Tukey’s honestly significant difference (HSD) test for multiple comparisons was performed to determine the significance of differences in floral traits between the two sites in the same flowering phase or between the two flowering phases at the same site. Pearson correlation coefficient was calculated to determine the correlation among the floral traits of different morphs.

The indices developed by Sánchez et al. (2008) were used to compare the degree of spatial matching between the reciprocal arrangements of reproductive organs in different populations. These indices compare each organ position with every other organ position in the sample and provide a single combined metric of reciprocity (R) for the lower- and upper-level reproductive organs. These indices can be interpreted as a measure of the average population-level deviation from perfect reciprocity, with a value of 1 indicating perfect reciprocity, and values > 0.5 indicating distyly (Sánchez et al., 2008, 2013).

A generalized linear mixed (GLM) model, with Gaussian distribution and restricted maximum likelihood (REML), was used to test whether the seed set rate of plants covered with or without a net differed significantly between sites (central vs. marginal), flowering phases (early vs. late) and morphs (L- vs. S-morph). In this model, space, time, morph, and space × time were entered as fixed factors. Bonferroni correction was performed for multiple comparisons to determine the significant differences in seed set rate among different experimental treatments. A linear least regression model was used to test the relationship between herkogamy and the seed set rate, assuming a normal distribution of the errors. All statistical analyses were performed in SPSS version 19.0 (IBM Inc., Armonk, New York, United States).

The overall morphological traits of S- and L-morphs of P. wannanensis showed significant spatial variation [S-morph: Wilks’ λ = 0.32, F_{(6, 630)} = 227.86, P < 0.0001; L-morph: Wilks’ λ = 0.37, F_{(6, 594)} = 169.66, P < 0.0001] and temporal variation [S-morph: Wilks’ λ = 0.48, F_{(6, 630)} = 114.53, P < 0.0001; L-morph: Wilks’ λ = 0.44, F_{(6, 594)} = 125.35, P < 0.0001], along with significant space × time interaction (Table 1), indicating that site-specific differences in floral traits vary between different flowering phases. The results of PCA also showed obvious differences in floral traits of S- and L-morphs between the central and peripheral sites and between early and late flowering phases; the floral traits of the S-morph in marginal populations changed significantly during the late flowering phase (Figures 1, 2). In addition, the first and second principal components (PC1 and PC2, respectively) of all floral traits, except CD, contributed much to these spatiotemporal variations (Supplementary Tables 2, 3). Univariate tests indicated that all floral traits differed significantly between the central and peripheral sites (P < 0.05) and between early and late flowering phases (P < 0.05), and all traits, except CD of the S-morph, were significantly affected by the site × flowering phase interaction (Supplementary Table 4). Compared with central populations, both floral morphs in peripheral populations showed a significantly higher positioning of stigmas and anthers, with longer corolla tubes in the early opening flowers, but without larger herkogamy (conversely smaller in the L-morph), resulting in significantly smaller RHE value in peripheral populations (Figure 3). Compared with early opening flowers, the late-opening flowers showed significantly reduced CD, CTL, and SH in both morphs, and this declining trend was more striking in peripheral populations than in central populations (Figure 3). However, pistil showed significant differences between the two floral morphs; the SH of the L-morph decreased significantly, especially in marginal populations, whereas that of the S-morph decreased in central populations but increased significantly in marginal populations (Figure 3). Consequently, in central populations, the HE value of late-opening flowers decreased from 1.65 to 1.45 mm in the S-morph and from 1.53 to 1.27 mm in the L-morph, but the RHE value showed no significant decline in the S-morph (P > 0.05) and only a marginally significant decrease in the L-morph (P = 0.011). By contrast, in peripheral populations, the HE value of late-opening flowers decreased significantly from 1.35 to 0.58 mm in the L-morph and from 1.62 to 0.33 mm in the S-morph, resulting in a significant decline in the RHE value (Figure 3).

Floral trait variation during the early and late flowering phases in central populations showed discrete dimorphism, and all flowers showed a certain level of herkogamy, which changed from 0.57 to 3.15 mm (Figure 4). However, floral trait variation in late-opening flowers in peripheral populations appeared continuous, based on the obscure dimorphism variation pattern, and herkogamy disappeared in nearly one-third individuals (34.05% in YLD, and 26.09% in QDX; Figure 5). Consequently, the global indices of reciprocity (R), calculated according to Sánchez et al. (2013), were high (0.89–0.92) during the early and late flowering phases in central populations (HDY: 0.88, GCN: 0.87; Figure 4) but significantly lower during the late flowering phase in peripheral populations (YLD: 0.79, QDX: 0.77, t₆ = 6.50, P < 0.001; Figure 4). The analysis of morphological data measured in 2018 (to investigate the influence of floral traits on seed setting ability) revealed that the spatiotemporal variation pattern in 2018 was consistent with that in 2017 (Supplementary Tables 6, 7 and Supplementary Figure 2), indicating that these spatial and temporal patterns of floral trait variation observed in P. wannanensis are stable.

Regardless of the presence or absence of effective pollinators, the seed set rate was significantly affected (P < 0.001) by space (central vs. peripheral site) and time (early vs. late flowering phase) but not by the morph (L- vs. S-morph) (Table 2). In addition, the interaction between space and time significantly (P < 0.001) affected the seed set rate only when plants were covered with a net to exclude pollinators (Table 2), indicating site-specific differences in seed set rate vary between different flowering phases in the absence of pollinators. When effective pollinators were excluded, the seed set rate in central populations was very low (average range: 0.05–0.08) and showed no significant differences between morphs (L- vs. S-morph) and between flowering phases (early vs. late) [F_{(3, 274)} = 0.41, P > 0.05; Figure 6]. The seed set rates of both morphs during the early flowering phase were also very low in peripheral populations (average range: 0.08–0.09), similar to those in central populations [F_{(5, 381)} = 0.60, P > 0.05]. However, the seed set rates of both morphs were quite high during the late flowering phase (average: 0.52 in the S-morph, 0.53 in the L-morph), which were significantly higher than those during the early flowering phase or those in central populations [S-morph: F_{(1, 107)} = 94.77, P < 0.001; L-morph: F_{(1, 102)} = 101.63, P < 0.001; Figure 6]. Under natural conditions, the seed set rates of both morphs in the central and peripheral populations were relatively high (means changed from 0.77 to 0.89; Figure 6). In central populations, the seed set rates of both morphs during the early flowering phase were similar to those during the late flowering phase [F_{(3, 156)} = 2.42, P > 0.05]. The seed set rates of both morphs during the early flowering phase in peripheral populations (S-morph: 0.89 and L-morph: 0.87) were significantly higher than those in central populations [S-morph: 0.80, F_{(1, 78)} = 15.37, P < 0.001; L: 0.77, F_{(1, 78)} = 23.72, P < 0.001] and those during the late flowering phase [S-morph: 0.80, F_{(1, 78)} = 12.07, P < 0.01; L-morph: 0.79, F_{(1, 78)} = 10.55, P < 0.01].

Under natural conditions, the seed set rates of both morphs showed no correlation with the herkogamy value (Figure 7). However, when effective pollinators were excluded with a net, the seed set rates of both morphs showed a significant negative correlation with the herkogamy value (L-morph: r = 0.62, P < 0.001; S-morph: r = 0.78, P < 0.001), and this negative correlation was stronger in peripheral populations (L-morph: r = 0.72, P < 0.001; S-morph: r = 0.83, P < 0.001) than that in central populations (L-morph: r = 0.31, P < 0.001; S-morph: r = 0.22, P < 0.05) (Figure 7). In peripheral populations, the average seed set rate was relatively high (range: 0.55–0.71; Figure 8) when the herkogamy value was less than 1.0 mm but decreased remarkably when the herkogamy value increased beyond 1.0 mm; average seed set rates of 0.18 and < 0.05 were observed at herkogamy values of 1.0–1.2 mm and > 1.2 mm, respectively (Figure 8).

In general, relatively sufficient pollinators and resources required for seed set are available at the peak or earlier flowering time, which is very important for most plants to execute sexual reproduction (Elzinga et al., 2007; Marshall et al., 2010). In our study system, variation in the position of sexual organs in both morphs in central and peripheral populations showed obvious discrete dimorphism (Figure 4) at the early flowering stage, with high global R indices (>0.87), indicating a good match between male and female organs. Furthermore, the seed set rate less than 10% of these flowers when effective pollinators are excluded with a net (Figure 6), indicating they were lack of effective automatic pollination mechanism and their effective pollen transfer normally requires pollinators activity. Consequently, regardless of the site, the early flowers of P. wannanensis, similar to those of most other distylous species, showed typical reciprocal herkogamy characteristics that promote disassortative pollination and reduce sexual interference (Ganders, 1979; Keller et al., 2014; Zhou et al., 2015; Barrett, 2019).

Highly significant correlations were detected between CTL and the positions of sexual organs in both morphs (Supplementary Table 5), as documented in other distylous Primula species, suggesting that corolla tube elongation plays an important role in the spatial arrangement of sexual organs in the Primula genus (Kálmán et al., 2007; Zhu et al., 2009; Keller et al., 2012; Brys and Jacquemyn, 2015). Interestingly, compared with central populations, marginal populations exhibited longer corolla tubes and higher positioning of stigmas and anthers in both morphs (Figure 3). Preliminary transplanting experiments in the greenhouse showed that this site-specific variation was stable under the same cultivation conditions (Zhang and Shao, unpublished data), indicating that these differences were a result of long-term local adaptive selection and represent heritable variation. Floral traits are largely attributed to phenotypic selection in plant–pollinator interactions (Boberg et al., 2014; Huang et al., 2016; Wu and Li, 2017). Previous studies showed that the two long-tongued bee species, B. major and A. chinensis, are effective pollinators of P. wannanensis, and the proboscis of the former bee species (average length: ∼5.78 mm) is significantly shorter than that of the latter (average length: ∼7.90 mm) (Shao et al., 2008). In the current study, we could not determine whether pollinator communities varied between central and peripheral populations. Further investigation is needed to determine whether the observed site-specific variation could be attributed to interactions with local pollinators, although coevolutionary between some key floral traits (such as corolla tube or stigma position) and local pollinators has been documented (Moeller, 2006; Paudel et al., 2015, 2016), even including some primroses (Brys and Jacquemyn, 2015; Wu and Li, 2017). It should be noted that the relatively narrow corolla tube, as in P. wannanensis (<2 mm diameter), may restrict the entry/exit path of the proboscis of effective pollinators foraging for nectar secreted at the bottom of the corolla. Higher position of sex organs in a floral tube may lead to a larger area of contact with the proboscis of pollinators, thus increasing pollination efficiency (Nishihiro et al., 2000). Thus, in peripheral populations, where pollinators were usually limited and unpredictable, higher CTL associated with higher stigma and anther positions may be selected to improve the efficiency of pollination, which may explain why the seed set rate of early flowers under natural conditions was slightly higher in peripheral populations than in central populations in 2018 (Figure 6).

Although most Primula species have several inflorescences opening in succession, and thus have a relatively long flowering period (Hu and Kelso, 1996; Richards, 2003), floral trait variation in these distylous species during the flowering period and its reproductive implications have not yet been reported. Here, we documented a strong temporal difference in flower morphology in P. wannanensis. Compared with early opening flowers, the late-opening flowers showed significantly smaller reproductive organs (Figure 3), consistent with the major trend variation revealed in other sequentially blooming plants (Guitián and Navarro, 1996; Diggle, 2003; Bateman and Rudall, 2006). Furthermore, this temporal variation pattern of flowers showed conspicuous spatial differences. In the central populations, the measured floral traits of distylous such as (SH,AH, CTL, HE, and CD) all most scaled down in late opening flowers, thus causing the herkogamy near proportionally declined with floral tube (Figure 3). However, in peripheral populations, the SH variable did not scale down with other floral traits (AH, CTL, and CD), but the style was elongated in the S-morph and shortened in the L-morph (Figure 3), resulting in significantly decreased values of both HE and RHE (Figure 3). In addition, at the early flowering stage, the flowers outside the net could form seeds, while those covered with a net hardly formed seeds because of the exclusion of effective pollinators. However, at the late flowering stage, we did not document that the size of the flowers inside the net was significantly larger (P > 0.05) than that outside the net. Thus, whether the early opening flowers bear seeds or not had little effect on the size of late-opening flowers, suggesting that architectural effects, rather than phenotypic plasticity, could explain the decline in the size of late-opening flowers (Stephenson, 1981; Diggle, 1995, 2003).

The site-specific temporal variation in floral traits also led to a remarkable variation in the reproductive strategy of late-opening flowers. In central populations, the late-opening flowers maintained typical distylous floral characteristics (with relatively high reciprocal herkogamy, R > 0.86; Figure 4), and could not bear seeds when pollinators were excluded (Figure 6), implying that the floral traits of late-opening flowers promoted outcrossing. By contrast, in peripheral populations, the late-opening flowers did not show typical distylous floral traits, and could produce seeds by automatic self-pollination when pollinators were excluded because of the strong reduction or disappearance of herkogamy (Figure 6). Notably, the arrangement of sexual organs in the most late-opening flowers in marginal populations has changed to favor automatic selfing. The occurrence of self-pollination in only peripheral populations is consistent with the prediction that peripheral populations are more prone to suffer from pollinator scarcity or conspecific plants, and reproductive assurance is more likely to play a role at the ecological edge of the species distribution range (Brown, 1984; Busch, 2005; Griffin and Willi, 2014). Among a variety of self-pollination modes, delayed autonomous selfing is generally regarded as an adaptive strategy, as it provides reproductive assurance without incurring pollen and/or seed discounting costs (Schoen and Brown, 1991; Lloyd, 1992). Therefore, many studies report that species possess herkogamous flowers with sundry mechanisms to bring the reproductive organs together usually toward the end of the floral life cycle, e.g., by corolla dragging, which involves movement of the stigma toward the anthers or vice versa (Fenster and Martén-Rodríguez, 2007). In our current study, late-opening flowers were facilitated to undergo autonomous self-pollination, which differs from the delayed self-pollination that generally occurs at the end of the flower life cycle, the ecological significance of both mechanisms is the same. For P. wannanensis, the seed set rate of open-pollinated flowers (especially early opening flowers) in peripheral populations was high (Figure 6) and didn’t lower than that in central populations, raising doubts about the supposed limited pollinator service in peripheral populations (at least in 2018). Previous manual pollination experiments also revealed that outcrossing pollen grains possess a higher germination rate and pollen tube growth than that of selfing (Shao et al., 2011). Moreover, peripheral populations still contained relatively high genetic diversity and normal flower morph ratio (Shao et al., 2008, 2019). All these indicated that most of the seeds set by late flowers in peripheral populations might well result from outcrossing under natural conditions, although the actual contribution to autonomous pollination to seed set in late flowers in peripheral populations will require further molecular analysis of the offspring. Thus the autonomous selfing of late-opening flowers in peripheral populations serves as only a potential (or standby) reproductive assurance mechanism. Only when pollinators are scarce throughout the flowering period, the seeds produced by this self-pollinating system can be promoted, which provides a reproductive assurance and maintains population stability. Thus, the autonomous self-pollination capacity of late-opening flowers in marginal populations of P. wannanensis is another delayed reproductive assurance mode.

In typical distylous species, mating patterns are rarely affected by herkogamy variation because of their strict sporophytic self-incompatibility, although small changes in reciprocal herkogamy significantly influence the efficiency of legitimate pollen deposition on stigmas (Ganders, 1979; Brys and Jacquemyn, 2015; Liu et al., 2015). According to our current observations, we conform this prediction, i.e., the selfing seed set rate is low and effected by herkogamy variation is very limited in central populations (Figure 5), which possess a relatively strict self-incompatibility system (Shao et al., 2019). However, it is largely unknown whether and to what extent the variation in herkogamy affects opportunities for autonomous selfing and reproductive assurance in heterostylous plants with self- or partially self-compatible systems, although the self-incompatibility of quite a few distylous species is not strictly (Riveros et al., 1995), even in the typical distylous genus Primula (Wedderburn and Richards, 1990; Barrett and Cruzan, 1994; Riveros et al., 1995; Yuan et al., 2017; Zhou et al., 2017; Jiang et al., 2018; Shao et al., 2019). In the marginal populations of P. wannanensis, self-incompatibility is fairly weak (Shao et al., 2019), and we found that the degree of herkogamy directly affected the autonomous selfing capacity of plants in the absence of pollinators (Figure 7). Further analysis revealed that the average seed set rate via selfing was relatively high (> 0.55 for all herkogamy classes; Figure 8) and showed no significant difference among herkogamy classes with separation less than 1.0 mm; however, the seed set rate sharply decreased to less than 0.10 with an increase in separation (Figure 8). This sharp decline in the seed set rate could be explained based on the assumption that the HE value of 1.0 mm possibly represents an approximate threshold, consistent with previous studies on the homostylous species P. halleri (de Vos et al., 2014, 2018) and distylous species Primula chungensis (Jiang et al., 2018), suggesting that this threshold likely has wider applicability, although experimental confirmation is needed.

Autonomous selfing is a common strategy employed by plants to accommodate an unpredictable pollinator environment because selfed seeds can provide reproductive assurance (Lloyd, 1992; Barrett, 2002). The spatial and/or temporal separation of female and male organs within flowers of self-compatible plants is generally assumed to have a direct impact on the degree of selfing and the capacity to autonomous self-pollination (Barrett, 2002; Nishihiro et al., 2000; Brys and Jacquemyn, 2011; de Vos et al., 2014, 2018). In distylous species, a commonly exhibited heteromorphic incompatibility system prevents self- and intramorph mating, and a significant component of phenotypic disassortative mating is maintained under the interaction of pollinator and reciprocal herkogamy in distylous species exhibiting self- and/or intramorph compatibility (e.g., Luculia pinceana; Zhou et al., 2015). Most heterostylous groups contain species that are monomorphic for style length, with anthers and stigmas positioned close together within a flower. These plants are generally self-compatible and predominantly exhibit selfing as a result of autonomous self-pollination (Ganders, 1979). As the main form of selfing variants, homostylous species are evolutionarily derived from distylous ancestors. In Primula, phylogenetic analyses demonstrated a single origin of distyly but numerous independent transitions to homostyly (Mast et al., 2006; de Vos et al., 2014; Zhong et al., 2019). Investigations of the genetic basis of homostyly in Primula indicated that the rapid transition from outcrossing to selfing was caused by loss-function mutations in a single-gene (CYP734A50), which controls the stigma height in the S-morph (Huu et al., 2016; Kappel et al., 2017). Consistent with the theoretical expectations of the genomic selfing syndrome (Cutter, 2019), a recent analysis of genome-wide molecular evolution revealed the maladaptation of autonomous selfing in homostylous species with strongly reduced genetic diversity, decreased purifying selection efficacy, and low adaptive evolution rates (Wang et al., 2020).

Delayed autonomous selfing has not been previously reported in distylous taxa, although studies have suggested different levels of autonomous selfing in derived lineages of homostylous species including Primula (Yuan et al., 2017; Zhong et al., 2019). Here, we report a novel example of delayed autonomous selfing occurring in the marginal populations of distylous P. wannanensis. Previously, Shao et al. (2019) reported that self-incompatibility was largely weakened in the marginal habitat of P. wannanensis; however, it is unclear what mode of selfing has been adopted to cope with the unstable pollination environment, i.e., pollinator facilitated selfing/geitonogamy, prior autogamy, simultaneous autogamy or delayed autogamy. All these options may provide the advantage of reproductive assurance, as long as the fitness of selfed progeny (ω_s) to that of outcrossed progeny (ω_o) exceeds 0.5 (ω_s/ω_o > 0.5), i.e., inbreeding depression (δ = 1–ω_s/ω_o) is less than 0.5. In this study, late-opening flowers in marginal populations showed a high seed set rate under pollinator-excluded conditions (Figure 6); however, the early opening flowers maintained a high level of reciprocal herkogamy and consequently pollinator-mediated disassortative mating (i.e., outcrossing) (Figure 4). This implies a transition of mating strategy from obligate outcrossing in central populations to mixed-mating in marginal habitats, further suggesting that the selfing component of mixed-mating is more likely to occur in late-opening flowers, which is thus referred to as delayed autonomous selfing. This transition was distinctly different from the most common events of the mating system shift from outcrossing to selfing associated with the breakdown of distyly to homostyly. Predominant self-fertilization of homostyly may result in the accumulation of genetic load, leading to a “dead-end” over the long run (Wang et al., 2020). However, the ability to exhibit delayed autonomous selfing may represent a floral strategy that limits the most harmful genetic consequences of selfing and increases the long-term adaptive potential of the species.