Intraspecific variation of self-incompatibility in the distylous plant Primula merrilliana

Heteromorphic self-incompatibility can prevent self- and intramorph fertilization while favouring intermorph mating and the maintenance of morph-ratio stability in heterostylous populations. However, variation in the expression of self-incompatibility intraspecies has seldom been assessed. Through hand pollinations and microsatellite markers, the variation in the expression of self-incompatibility and genetic diversity were studied in distylous plant Primula merrilliana. We discovered that the strength of self-incompatibility varied extensively among individuals and populations, from pronounced to weaker self-incompatibility in distylous populations, all the way to strong self-compatibility in homostylous populations. Each distylous population included self-incompatible (SI), partly self-compatible (PSC) and self-compatible (SC) individuals, with the index of self-compatibility (ISC) ranging from 0.07 to 0.68 across populations. Self-compatible populations (ISC > 0.25) were not genetically clustered but were more closely related to populations with high SI and SC individuals were mixed with SI individuals within populations. The ISC and the proportions of SC and PSC individuals were higher in peripheral than in central populations, but no decrease of genetic diversity and no deviations of floral morph ratio from isoplethy were detected in peripheral populations. Additionally, the expression of self-incompatibility was stronger in long-styled flowers than in short-styled flowers. The variation in the strength of self-incompatibility documented in P. merrilliana cautions against the estimation of ISC from a few individuals or populations in distylous species and provides a more complex and nuanced understanding of the role of self-incompatibility in heterostyly.


Introduction
Heterostyly, a term referred to plant species comprising either two (distyly) or three (tristyly) self-incompatible floral morphs that differ reciprocally in the positions of stigmas and anthers (i.e. reciprocal herkogamy), is considered as one of the most effective mechanisms to avoid selfing and promote outcrossing (Darwin 1877;Barrett 1992). This complex floral syndrome, which occurs in at least 28 families of flowering plants (Ganders 1979;Lloyd and Webb 1992), represents a remarkable example of convergent evolution. Therefore, explaining the origin, function and genetic basis of heterostyly has attracted considerable attention ever since Darwin's classic book on polymorphic sexual systems in plants ( Darwin 1877;Barrett 1992).
Reciprocal herkogamy is usually associated with a sporophytic self-incompatibility system that prevents self-and intramorph (assortative) mating and promotes intermorph (disassortative) mating, while maintaining morph-ratio stability in populations (Ganders 1979;Barrett 1992; Barrett and Cruzan 1994). Joint investigations into the variation of physiological self-incompatibility among populations (individuals) and the genetic relationships are crucial to enhance our understanding of the evolution of the heterostylous syndrome (Barrett and Cruzan 1994), since the two most noteworthy evolutionary models explaining its origin differ in hypothesizing that heteromorphic incompatibility either preceded or followed the evolution of reciprocal herkogamy (Charlesworth and Charlesworth 1979;Lloyd and Webb 1992). Studies of seed set following controlled pollinations have revealed that not all heterostylous plants possess strict self-incompatibility and wide variation exists in the degree of self-incompatibility (Barrett and Cruzan 1994). For example, in heterostylous Rubiaceae, the expression of selfincompatibility exhibits considerable variation, ranging from conventional diallelic self-incompatibility (Carlson et al. 2008) to cryptic self-incompatibility (Wu et al. 2010) and complete self-compatibility (Riveros et al. 1995). Even in Primula, the classical model for heterostyly, self-incompatibility varies greatly among distylous species, from the more common self-incompatibility to the less frequent self-compatibility, as documented, for example, in P. conchifolia, P. calderiana, P. pulverulenta, P. chungensis and P. oreodoxa (Wedderburn and Richards 1990;Yuan et al. 2017;Zhou et al. 2017). This wide variation in the expression of self-incompatibility between related taxa may offer opportunities for re-evaluating previously hypothesized models on the evolution of heterostyly inferred from phylogenetic analysis (Mast et al. 2006;Pérez-Barrales et al. 2006; Barrett and Shore 2008;de Vos et al. 2014). A prerequisite to correctly interpret the evolutionary value of interspecific variation in self-incompatibility is to first determine whether it is stable within species, because if self-incompatibility varies considerably among or within populations in a species, its strength quantified from only one or two populations may not correctly reflect that of the species, hence invalidating estimates of its variation between taxa. The transition from self-incompatibility to full selfcompatibility concomitant with the morphological transition from heterostyly to homostyly is relatively well documented (Crosby 1940;Dowrick 1956;Piper et al. 1984;Boyd et al. 1990;Wadi and Richards 1993;Mast et al. 2006;Carlson et al. 2008). Controlled pollination experiments have shown that a large amount of variation can also exist in the expression of trimorphic incompatibility among individuals and morphs within a species (Ornduff 1972;Barrett and Anderson 1985;Barrett and Shore 1987). However, to our knowledge, there are few empirical studies on the variation of selfincompatibility among populations within distylous species, even though distyly is more common than tristyly (but see Sobrevila et al. 1983). Unlike the case for heteromorphic self-incompatibility, variation in homomorphic self-incompatibility has been comparatively well investigated within species. Empirical studies suggest that the major pathway of variation within species, similar to evolutionary trends between species, is the transition from self-incompatibility, which enforces outcrossing, to self-compatibility, which enables selfing (Busch 2005;Mable et al. 2005;Mable 2008;Koelling et al. 2011;Griffin and Willi 2014). Besides genetic transmission advantage (a selfing genotype transmits twice as many copies of its genes to offspring as an outcrossing genotype), reproductive assurance is the most widely accepted ecological force favouring selfing (Fisher 1941;Jain 1976). Because peripheral populations of a species are more prone to suffer from scarcity of pollinators or conspecific plants, reproductive assurance is more likely to play a role in peripheral environments; hence, self-compatible individuals and populations are mostly located at the ecological edge of a species distribution (Brown 1984;Busch 2005;Griffin and Willi 2014).
In heterostylous plants, the pistil and stamen are separated and reciprocally positioned, so pollen transfer normally requires pollinator activity (Keller et al. 2014;Armbruster et al. 2017). Furthermore, heteromorphic self-incompatibility restricts mating opportunities to one-half (distyly) or two-third (tristyly) of the plants in a population. Compared to homomorphic species, heterostylous species hence are more sensitive to scarcity of pollinators and conspecific individuals, and may more strongly need the benefits of reproductive assurance in peripheral populations. This consideration also explains the fact that self-compatible populations with homostylous flowers often occur at the geographical margins of a heterostylous species range (Barrett 1985;Barrett and Shore 1987;Barrett et al. 2009;Yuan et al. 2017) or that homostylous species recolonize peripheral areas after glacial retreat (Guggisberg et al. 2006(Guggisberg et al. , 2009Theodoridis et al. 2013;de Vos et al. 2014). Thus, we hypothesize that if self-incompatibility varies among distylous populations within a species, the loss (or partial loss) of self-incompatibility will predominantly occur in peripheral populations. The transition of self-incompatibility system in homomorphic species usually accompanied by change of mating pattern (from outcrossing to selfing) and the decrease of genetic diversity (Mable et al. 2007). In heterostylous plants, however, selfing may not easily occur even with the loss of self-compatibility due to separation of pistil and stamen morphology (de Vos et al. 2018).
Thus, it is not clear whether the loss of self-incompatibility in heterostylous plants will lead to a transition towards self-fertilization.
Primula merrilliana (Primulaceae) is a distylous herb with short-lived life history (8-9 months) and can be relatively easily cultivated in the greenhouse. Thus, it is a suitable material to detect whether self-compatibility varies within a species through pollination experiments. In this study, we attempt to clarify the following questions through manual pollination experiments and population genetic analyses with nuclear microsatellite (simple sequence repeat, nSSR) markers: (i) Does self-incompatibility vary among individuals, populations and between morphs within P. merrilliana? (ii) If so, does the loss of self-incompatibility occur mostly in peripheral populations? (iii) Does the variation in self-incompatibility result in departure from 1:1 morph ratios? (iv) What patterns of genetic relationships between SC and SI individuals and populations?

Study species and sampling
Primula merrilliana plants are characterized by 10-30 pinnate compound leaves and 20-50 white or lilac flowers in 3-10 inflorescence stalks. Flowering usually begins in late February and ends in the middle of May. Dehiscent capsular fruits ripen in May-June. Primula merrilliana usually occurs on shaded slopes under or at the edge of deciduous, broad-leaved forests between 50 and 1000 m a.s.l. (Hu and Kelso 1996;Shao et al. 2008aShao et al. , 2012Chen 2009). Extensive field surveys in recent years show that P. merrilliana is now restricted to Huangshan mountains (extending westward to Mt. Guniujiang) and its surrounding areas in southern Anhui Province (eastern China), of which the Huangshan mountains are the geographical distribution centre (Fig. 1).
To measure the strength of self-incompatibility, during October-November 2011, 30-60 seedlings per population from 11 populations representing the entire distribution area of the species were randomly selected for transplantation into a greenhouse for controlled pollination experiments (Fig. 1). Only seven individuals were collected from the QD population due to its small census size of ~100 individuals, as opposed to ~700-5000 individuals in the remaining populations. Among these sampled populations, six populations (TJQ, THF, TJ, GCD, GCN and GCI) were classified as central populations due to they are nearer to the distribution centre Huangshan mountains than to the distribution boundary, and the remaining five populations (QYS, PLD, YLD, QD and TL) as peripheral populations with smaller distance to the edge of the distribution area than to the distribution centre.

Variation of self-incompatibility
Pollinations were performed of the newly opened flowers during March-April in 2012. From each of the 11 populations mentioned above, 30-53 (only seven in QD) plants (for a total of 427 plants) were selected for manual pollinations, which were carried out in the greenhouse without any pollinators. Each plant was subjected to two treatments: self-pollination and intermorph cross-pollination. Self-pollinations were performed on at least four flowers per plant: using tweezers, we gently pulled apart the tube of each flower and brushed the pollen on the stigma of the same flower. Intermorph cross-pollinations were carried out by collecting pollen from one or two randomly chosen plant donors with opposite morph (i.e. L-morph × S-morph or S-morph × L-morph) within populations; each treated flower had been emasculated by gently removing the anthers with tweezers prior to treatment. Intermorph pollinations were performed on at least two flowers per plant. However, the anthers and stigma of two populations (TL and QD) were at the same position of corolla tube, both occurred at the mouth (in TL) or at the middle (in QD). Their flowers can complete self-fertilization early in anthesis, and it is very difficult to emasculation prior to anthesis. Thus, the plants in these two populations were only to be treated by self-pollinations.
Ripe capsules were collected just prior to dehiscence. Full seeds, abortive seeds and unfertilized ovules were counted under stereomicroscope. Seed-set rates for self-and intermorph pollinations were calculated as the mean number of full seeds/mean number of ovules (full seeds + abortive seeds + unfertilized ovules) per plant. Given that the abortive seeds were probably caused by problems such as inbreeding depression or undernutrition after fertilization (Mahy and Jacquemart 1999;Yang et al. 2003), the index of self-compatibility (ISC) was calculated by dividing the mean number of seeds (full and abortive seeds) produced from self-pollinations over that produced from cross-pollinations per treated plant (Good-Avila and Stephenson 2002). Plants were classified as strongly selfincompatible (SI) if ISC was <0.20, as partly self-compatible (PSC) if ISC was between 0.20 and 0.80 and as fully self-compatible (SC) if ISC >0.80, following Good-Avila and Stephenson (2002).

Microsatellite genotyping
During flowering, one or two leaves were collected from each plant transplanted from the field into the greenhouse and dried in silica gel. Subsequently, 10-12 plants per population were selected for genetic analysis based on the results of controlled pollination experiments, as follows. If only one to five selfcompatible (SC) individuals were detected in each population, all of them were included in genetic analyses (except in YLD population, where six plants were randomly selected out of the 22 SC plants); if more than five partially self-compatible (PSC) and self-incompatible (SI) individuals, respectively, were detected in each population, two to six plants of each type were randomly selected. Total genomic DNA was extracted from leaf tissue using the CTAB method (Doyle 1991). Genotypes of DNA samples were scored using 10 pairs of microsatellites labelled with fluorescent-dye PCR primers (Table S1 in Shao et al. 2015), and amplified by polymerase chain reaction (PCR) following protocols described in Peng et al. (2009). PCR products were separated on an ABI3730 DNA Analyzer (Applied Biosystems, Foster City, CA, USA). Alleles were scored manually with the aid of Genemarker software version 1.91 (Applied Biosystems, Foster City, CA, USA).

Genetic diversity and structure analysis
For each population, mean number of alleles (R o ), observed (H o ) and expected heterozygosity (H e ), and inbreeding coefficient (F IS ) were estimated with FSTAT 2.9.3 (Goudet 1995). Using the same software, the Hardy-Weinberg equilibrium (HWE) of each locus and genotypic disequilibrium for all locus pairs were tested in each population by randomization, and the obtained P-values (<0.05) were adjusted applying a sequential Bonferroni correction to avoid false positives. Bayesian clustering implemented in STRUCTURE version 2.3.1 was used to assign individuals to genetic clusters (K) and to estimate admixture proportions for each individual (Pritchard et al. 2000). The number of K was set to vary from 1 to 15. For each value of K, we performed 10 independent replicates with a burn-in of 50 000 and a run length of 100 000 iterations. The analyses were run using the admixture model with independent allele frequencies.
The posterior probability of the data [ln P(D)] for a given K was computed using the runs with highest probability for each  Table 1. K (Pritchard et al. 2000). We then used structure harvester to summarize the results, including the identification of optimal K through the method of Evanno et al. (2005). Results across replicate runs were combined using the program CLUMP (Jakobsson and Rosenberg 2007) and visualized with DISTRUCT v1.1 (Rosenberg 2004). We also performed a phylogeographic analysis using the NeighbourNet method implemented in SplitsTree4 (Huson and Bryant 2006) based on the uncorrected p-distance between individuals to reveal the overall genetic pattern of relationships among populations and individuals, and constructed a consensus Neighbour-Joining tree based on pairwise estimates of Nei's genetic distance among populations, using 5000 bootstrap trees and random input order (PHYLIP v3.69;Felsenstein 2009). In addition, to visualize the genetic relationships among all individuals, their Euclidean distance matrix was subjected to a principal coordinate analysis (PCoA) using GENALEX 6 (Peakall and Smouse 2006). Analysis of molecular variance (AMOVA) in ARLEQUIN v 3.5.2.2 (Excoffier and Lischer 2010) was used to quantify the partitioning of genetic variance among and within populations. The significance of each variance component was tested through 10 000 permutations. We also calculated genetic differentiation (F ST ) between populations using ARLEQUIN. All statistical analyses were performed in SPSS version 13.0 (SPSS Inc., Chicago, IL, USA).

Population characters
The six central populations were located at higher altitude (from 245 to 850 m a.s.l.) than the five peripheral populations (from 100 to 210 m a.s.l.) (Table 1); the census sizes of the former populations were relatively high (from 700 to 5000 individuals, with mean above 2000), while those of the latter were comparatively lower (from 100 to 2000 individuals, with mean <1000), but the difference between them was not significant (Mann-Whitney U-test: Z = −1.203, P = 0.229). Two peripheral populations TL and QD were homostylous (i.e. exclusively consisted of long homostylous plants in TL and short homostylous plants in QD), while the remaining nine populations were distylous with L-morph ratios ranging from 0.34 to 0.60, but not significantly deviating from 0.5 (G-test: all P > 0.05).
In all distylous populations, the strength of selfincompatibility varied among individuals and all included SI, PSC and SC individuals (Figs 1 and 3). However, the frequency of SI, PSC and SC individuals differed significantly among populations (chi-square test: χ 2 = 112.615, P < 0.001), with >80 % self-incompatible individuals in TJQ, THF and GCI (of the central populations) and <30 % self-incompatible individuals in QYS and YLD (of the peripheral populations). In general, central populations harboured higher proportions of SI plants and lower proportions of SC plants (accounting for 72 and 3 % of the total number of studied individuals, respectively) than peripheral populations (χ 2 = 88.567, P < 0.001), which included a mean of 24 % SI and 20 % SC individuals, respectively (Figs 1  and 3).
For all investigated distylous populations, the number of ovules per flower and seeds per fruit resulting from intermorph pollinations did not significantly differ between floral morphs (Table 2; see Supporting Information-Table S1). However, following self-pollination, the mean number of seeds Table 1. List of sampled populations with the results of manual pollination experiments and analyses of genetic diversity. Note: N, populations size; D, distyly; SH, short homostyly; LH, long homostyly; n, sample size for individuals cultivated in the greenhouse and used for subsequent analyses; L-ratio, long-styled morph ratio in randomly selected individuals under cultivation; self, mean (± SD) seed set rate of self-pollinations; cross, mean (± SD) seed set rate of intermorph cross-pollinations within population; ISC, index of self-compatibility inferred from seed-set results; R o , mean number of alleles; H o , observed heterozygosity; H e , expected heterozygosity: both indexes were inferred from microsatellite data; F IS , inbreeding coefficient. In columns 'self' and 'ISC', numbers followed by different letters differ significantly among treatments (P < 0.05, Scheffé's post hoc test), in F IS column, asterisks denote significantly deviated from HWE (*P < 0.05; **P < 0.01; ***P < 0.001).  Table  S1). The proportions of SI, PSI and SC individuals between two morphs within populations were roughly equal and not significantly different from each other across populations (P > 0.05), except that the proportion of PSC individuals was higher in the S-morph and the proportion of SI individuals was higher in the L-morph for the GCN population (χ 2 = 16.718, P < 0.001; Fig. 3).

Genetic diversity and structure
The total number of alleles per locus in our sample of 123 individuals across 11 populations ranged from 3 to 23, with a total of 146 alleles distributed among the 10 microsatellite loci. and H e ) and population characteristics (i.e. population altitude, population size, HSI and ISC) were found. The inbreeding coefficients (F IS ) of the two homostylous populations (QD and TL) were high (0.859 and 1.000, respectively), while those of the distylous populations ranged from 0.058 to 0.391. Except for TJQ, all the studied populations significantly deviated from HWE, and the mean inbreeding coefficient of the peripheral populations was not significantly higher than that of central populations (Z = −1.643, P = 0.100). Analysis of molecular variance revealed that about half of the total genetic variance (48.08 %, P < 0.001) of P. merrilliana was among populations and F ST values between pairs of populations were high (mean 0.366, ranging from 0.127 to 0.798, all being significantly different from zero after sequential Bonferroni correction; see Supporting Information-Tables S2 and S3). When using the Bayesian clustering approach of Pritchard et al. (2000), the inference of the number of gene pools K was not straightforward, as log-likelihood values for the data conditional on K, lnP(D), increased progressively from K = 1 to K = 15 [see Supporting Information- Fig. S1a]. ΔK values (Evanno et al. 2005) computed for all K classes indicated a relatively strong signal for K = 11 [see Supporting Information- Fig. S1b], K = 2 and K = 6, with two other not very significant signals [see Supporting Information- Fig.  S1b]. When K equalled 11, the 11 sampled populations were relatively well separated from each other, with less than six individuals with probability above 0.5 of belonging to other populations, implying clear genetic differentiation among the sampled populations [see Supporting Information- Fig.  S2]. When K equalled 2, 6 or 11, SC and PSC individuals did not form a single or two clusters, but were intermixed with SI individuals in each population, which was also confirmed by the results of PCoA [see Supporting Information- Fig. S3]. The NeighbourNet network analysis (Fig. 4) showed that highly  Population codes are identified in Table 1. self-compatible populations (ISC > 0.25, i.e. TL, QD, YLD, QYS, TJ and PLD) did not form a single cluster, but occurred in different branches together with high SI populations (THF, TJQ and GCI), and also confirmed that SC individuals were intermixed with SI and PSI individuals.

Variation in the expression of self-incompatibility among populations
We documented extensive variation of SI in distylous populations of P. merrilliana, and even the transition to full SC and homostyly in two peripheral populations of this species. Although homostylous populations have been documented in several heterostylous species of Primula (Crosby 1940;Dowrick 1956;Piper et al. 1984;Wadi and Richards 1993;Yuan et al. 2017), these homostyles usually are confined to long homostyles, while short homostyles are very rare. It is noteworthy that, in P. merrilliana, both long homostyly and short homostyly occur, which may provide a unique opportunity for us to study the origin and maintenance of homostylous (especially short homostylous) populations within a distylous species. These two SC homostylous populations (TL and QD) were in different branches of the Neighbour-Net network and not closely related to each other (Fig. 4), suggesting they were the outcome of two independent evolutionary events. Additionally, the two homostylous populations of P. merrilliana occurred at the edge of its distributional range (Fig. 1), suggesting that reproductive assurance may play a major role in the evolution and maintenance of homostyly. In fact, self-compatibility and the intra-floral proximity between anthers and stigma of homostylous plants enable autonomous self-pollination also when pollinators and mates are scarce or absent, as typical in peripheral populations. Indeed, the self-pollination experiments performed under controlled conditions in the greenhouse confirmed that the two homostylous populations of P. merrilliana have a significantly higher ability to produce selfed seed than any of the distylous populations of the same species (Table 1). Also as expected, the values of observed genetic diversity (H o ) for TL and QD were significantly lower (0.000-0.043) than those of the distylous populations (0.250-0.558; Table 1), indicating that the genetic diversity of these homostylous populations has decreased noticeably, most likely due to their ability to self (Mable and Adam 2007;Willi and Maattanen 2011). Besides the loss of self-incompatibility accompanied by the morphological transition to homostyly documented for populations TL and QD, the expression of self-incompatibility varied remarkably also among distylous populations of P. merrilliana. Despite the crucial importance of selfincompatibility for the maintenance of balanced morph ratios in distylous species, the potential variation in the expression of SI has seldom been investigated. Our study reveals that the strength of self-incompatibility can in fact vary among populations of a distylous species. In P. merrilliana, the mean index of self-compatibility (ISC) ranged from 0.07 (in THF) to 0.68 (in YLD), and that of distylous peripheral populations was higher (ranging from 0.28 to 0.68, mean 0.48) than that of central populations (ranging from 0.07 to 0.30, mean 0.14).
The loss of SI has been documented in peripheral populations of species with a homomorphic type of SI (Busch 2005;Koelling et al. 2011;Griffin and Willi 2014), and reproductive assurance has been proposed as the major ecological driver for the transition to self-compatibility in peripheral populations also in those species. In populations of P. merrilliana characterized by heteromorphic SI and reciprocal herkogamy of anther and stigma positions, previously published field pollination experiments showed that the stigma inevitably comes into contact with self-and intramorph pollen, and the mean number of intermorph pollen grains deposited on the stigma was lower than that of its own ovules, especially in small and peripheral populations, but the total number of pollen grains was greater than that of its own ovules (Shao et al. 2008b(Shao et al. , 2011. Furthermore, studies in P. merrilliana (Shao et al. 2012) and other heterostylous species demonstrated that intramorph pollinations are always more compatible than self-pollinations (Sobrevila et al. 1983;Barrett et al. 2000b;Wu et al. 2010;Zhou et al. 2012), indicating that the transition to self-compatibility also implies intramorph compatibility. Thus, higher self-and intramorph compatibility  Table 1. Table 2. Two-way ANOVA for testing the effects of population and morph (L-morph vs. S-morph) on ovules per flower and seeds per fruit in self-and intermorph pollination. Note: df, degrees of freedom; MS, mean squares; asterisks denote significant differences (***P < 0.001). in peripheral, distylous populations of P. merrilliana should increase the number of potential mates, because a plant can mate with self, a plant of the same morph and a plant of the other morph, thus enriching the capacity for uniparental reproduction in peripheral colonizing situations and promoting reproductive success even under conditions of intermittently reliable pollinator service (Pannell et al. 2015). Usually, the transition from self-incompatibility to selfcompatibility and concomitant transition from outcrossing to selfing is accompanied by loss of genetic diversity (Mable and Adam 2007;Willi and Maattanen 2011), a prediction confirmed by the lower genetic diversity detected in the homostylous populations TL and QD. However, the mean genetic diversity of peripheral, distylous populations was not significantly lower than that of central, distylous populations, although the former have higher ISC than the latter. A possible explanation is that the more frequent occurrence of plants capable of producing selfed seeds in peripheral distylous populations might not immediately decrease genetic diversity, because such plants may still outcross at a rate sufficient to maintain high levels of genetic diversity. For example, in the homostylous, selfcompatible P. halleri, outcrossing can be achieved via approach herkogamy, as demonstrated by both pollination experiments and microsatellite analyses of genetic variation and progeny arrays (de Vos et al. 2012(de Vos et al. , 2014(de Vos et al. , 2018. Furthermore, previous manual pollination experiments in P. merrilliana revealed that the speed and rate of pollen germination on the stigma and speed of pollen tube growth in the style decrease progressively in intermorph, intramorph and self-pollinations (Shao et al. 2011); hence, SC and PSC individuals may self only when their stigmas do not receive sufficient quantities of inter-and intramorph pollen from other plants. Thus, the greater frequency of PSC and SC individuals in peripheral vs. central distylous populations might favour intramorph cross-fertilization over selfing, explaining the lack of significant differences in levels of heterozygosity between populations at the margin vs. the centre of the species distribution.

Variation in expression of self-incompatibility among individuals within populations
The expression of self-incompatibility varied remarkably also among individuals within populations and all distylous populations contained SI, PSC and SC individuals, but the proportions of SC and PSC individuals in peripheral populations were significantly higher than those in central populations.
Limited knowledge is available about the physiological and biochemical basis of incompatibility in heterostylous plants (Barrett and Cruzan 1994). The expression of self-incompatibility in the dimorphic P. merrilliana varied remarkably among individuals and populations, as documented in some trimorphic or homomorphic species (Barrett and Anderson 1985;Good-Avila and Stephenson 2002;Mable et al. 2005;Arroyo et al. 2012;Griffin and Willi 2014;Costa et al. 2017). Preliminary experiments to investigate the inheritance of self-incompatibility in P. merrilliana showed that the self-bred progenies of SC individuals included also SI and PSC individuals. And their proportions were similar to those of their parental population (unpublished data by J. Liu and J. W. Shao). Furthermore, the genetic relationship analysis showed that highly SC (PSI) populations and individuals were intermixed with SI populations and individuals (Fig. 4). These two phenomena suggest that self-incompatibility in P. merrilliana is probably a complex trait controlled by multiple genes, a conclusion that deviates from the classical, simple model of diallelic Mendelian control of the traits associated with distyly, including self-incompatibility, at the heterostyly supergene.

Morph-specific difference of self-incompatibility and morph-ratio maintenance
Morph-specific differences in the expression of selfincompatibility have been documented in several heterostylous species (Barrett and Cruzan 1994;Riveros et al. 1995;Larson and Barrett 1998;Nishihiro and Washitani 2011;Huang et al. 2015). But no consistent morph-associated patterns of variation for self-incompatibility were found by wider surveys (Wedderburn and Richards 1990;Barrett and Cruzan 1994). Our controlled pollination experiments revealed that the expression of selfincompatibility in L-and S-morphs varied among individuals and populations, but the proportions of SC, PSC and SI individuals were not significantly different between the two morphs within population, except in GCN, where they did differ significantly (Fig. 3). Across all distylous populations, however, the mean seed set from self-pollination in S-morph was significantly higher than in L-morph, differently from what proposed by Lewis (1942) ( Table 2; see Supporting Information- Table S1).
It is assumed that strong sporophytic incompatibility is necessary to prevent selfing and, coupled with reciprocal herkogamy, promote mating between floral morphs (Ganders 1979;Zhou et al. 2012). In turn, disassortative mating, in combination with heterozygous (Lewis and Jones 1992) or hemizygous genetic control of heterostyly (Huu et al. 2016;Kappel et al. 2017), ensures the maintenance of isoplethic equilibrium between morphs, due to negative frequency-dependent selection. If both self-and intramorph mating are permitted, floral morph ratios deviate from isoplethy, depending on the relative amounts of assortative mating in populations (Ganders 1979;Sobrevila et al. 1983;Barrett and Cruzan 1994;Barrett et al. 2000a, Barrett andShore 2008;Zhou et al. 2012). However, in P. merrilliana, all sampled populations have some PSC and SC individuals, especially in YLD, QYS, PLD and TJ populations. As discussed earlier, assortative mating (especially intramorph mating) may occur in these populations, but their morph ratios do not significantly differ from 1:1 (Table 1, equilibrium morph ratio; also see in Shao et al. 2008a). Cryptic incompatibility system and selective abortion of embryos were proposed to regulate distylous morph ratios (Weller and Ornduff 1977;Ganders 1979;Beerli 2008). Establishing whether these two processes contribute to maintaining equilibrium morph ratio in distylous populations of P. merrilliana requires further experiments. Here, we additionally propose that higher seed set in S-than L-morph possibly contributes to keeping the stability of equal morphs ratio. In Primula, it has been widely accepted that there exists a S-locus linkage group (supergene) governing distyly, with the L-morph being homozygous recessive (ss) and the S-morph being governed by a dominant S-allele (Ss) (Charlesworth and Charlesworth 1979;Lewis and Jones 1992), or with the S-locus being present exclusively in the genome of the S-morph, but not in that of the L-morph (hemizygous S-locus; Huu et al. 2016;Li et al. 2016;Kappel et al. 2017). If the heterozygous model of the S-locus is valid in P. merrilliana, the L-and S-morph ratios of the progenies from selfed or intramorph crosses in short-styled flowers will be 1:2 (when SS genotypes are lethal) or 1:3, while the progenies from selfed or intramorph crosses in long-styled flowers will all be L-morph. However, the higher seed set from self-and intramorph pollinations of S-flowers can help buffer the increase in the proportion of L-flowers in the population resulting from self-and intramorph pollinations of L-flowers, contributing to the temporary maintenance of isoplethy and distyly even when self-incompatibility becomes weaker in some populations.

Conclusion
Our study revealed that the expression of self-incompatibility can vary extensively among individuals and populations within a distylous species. In P. merrilliana, high intraspecific variation of self-incompatibility is displayed not only through the transition to full self-compatibility accompanied by the morphological shift from distyly to long or short homostyly, but also through high variation of pollen compatibility among individuals and populations that retain distyly. Therefore, the mosaic of variation in self-compatibility of this annual herb offers an ideal opportunity to investigate the mechanisms and causes of transitions between mating systems in distylous and derived homostylous plants. Importantly, our results show that estimates of self-incompatibility drawn from a few individuals in one or two populations should be regarded with extreme caution.

Sources of Funding
This work was supported by the National Natural Science Foundation of China (31570336 31170317) and the grant of the Key Laboratory of Biotic Environment and Ecological Safety in Anhui Province.

Conflict of Interest
None declared.

Supporting Information
The following additional information is available in the online version of this article- Figure S1. The results of genetic structure analyses by STRUCTURE soft. (a) Plot of mean posterior probability LnP(D) values of each K; (b) The corresponding ΔK statistics calculated according to Evanno et al. (2005). Figure S2. Bayesian clustering results of the STRUCTURE analysis for nSSR data of 123 selected individuals (11 populations) of Primula merrilliana. Histogram of the STRUCTURE analysis for the model with K = 2, 6 and 11 (showing the relatively higher ΔK). Each thin bar represents a single individual, which may be partitioned into K colours depending on the estimated multilocus membership in each of K clusters, where each colour represents the posterior probability of that individual belonging to a cluster. Individuals marked by circle dots were SC (circle solid dot) or PSC (circle hollow dot), and unmarked represent SI individuals. Population codes are identified in Table 1. Figure S3. Distribution of 123 Primula merrilliana individuals sampled from 11 populations along the first two axis of a factorial correspondence analysis (with inertia of 17.0 and 11.7 %, respectively) based on microsatellite genotypes. The symbols in big, middle and small size represent SC, PSC and SI individuals, respectively. Population codes are identified in Table 1. Table S1. Ovules and seed set of morph-specific differences in the sampled distylous populations. Table S2. Analysis of molecular variance (AMOVA) within/ among populations of Primula merrilliana. Table S3. Pairwise population differentiation among sampled populations.