3 RESULTS
3.1 Frequency variation of floral morphs in natural populations of Primula vulgaris
The frequency of homostyles ranged from 0% to 100% across the 22 sampled populations of P. vulgaris in Somerset, England (Table 1; Figure 2A). Ten populations were dimorphic heterostylous (L- and S-individuals; D01–D10), one was dimorphic with L- and H-individuals (D*11), ten were trimorphic (L-, S- and H-individuals; T01–T10), and one was monomorphic for the H-morph (M01). Morph ratios in dimorphic heterostylous populations did not deviate from isoplethy except in D02 and D08, which had an excess of L- and S-plants, respectively (G = 7.13, P = 0.007 and G = 6.11, P = 0.01, respectively; Supplementary material Table S3). In trimorphic populations, L-individuals were significantly more frequent than S-individuals except in T03, where S- and L-morph frequencies did not differ significantly (G = 0.13, P = 0.71; Supplementary material Table S3). Congruently, the ternary plot suggested a skewed trajectory towards a reduction of S-individuals when H-individuals were present (Figure 2B).
3.2 Molecular basis of the transition to homostyly
We sequenced ~97% of the 1,587 bp comprising the five exons of CYPᵀ from 17 S- and 27 H-individuals of P. vulgaris from seven and ten populations, respectively. In total, we detected seven CYPᵀ alleles (CYPᵀ-1 to -7; Figure 3A). All 17 S-individuals shared the same allele (CYPᵀ-1), identical to the functional copy of CYPᵀ previously reported in P. vulgaris (GenBank: KT257665.1; Li et al., 2016). Unexpectedly, the putatively functional CYPᵀ-1 allele also occurred in six H-individuals. In the 21 remaining H-individuals, we detected six CYPᵀ alleles (CYPᵀ-2 to -7), comprising altogether a single synonymous mutation, four nonsynonymous mutations, and two deletions. Specifically, CYPᵀ-4 contained two private mutations (one synonymous and one nonsynonymous) and five CYPᵀ alleles contained a single private mutation (three different nonsynonymous mutations in CYPᵀ-2, -3, -7 and two different deletions in CYPᵀ-5 and -6). The nonsynonymous mutation in CYPᵀ-2 introduced a premature stop codon in exon 2, whereas a 31 bp deletion in exon 5 of CYPᵀ-5 and an 8 bp deletion in exon 1 of CYPᵀ-6 introduced a frameshift leading to an early stop codon. Moreover, a mutation in exon 5 of CYPᵀ-3 changed a non-polar (phenylalanine) to a polar amino acid (serine), whereas a mutation in exon 3 of CYPᵀ-4 changed an amidic (asparagine) to a hydroxylic (serine) amino acid. Finally, the arginine to histidine mutation in CYPᵀ-7 caused no change in amino acid side-chain polarity. Overall, nucleotide diversity in homostyles was πS = 0.0005 and πNS = 0.001, while S-morph individuals showed no genetic variation in CYPᵀ (both πS and πNS were zero).
In the haplotype network, CYPᵀ-1, the known functional allele, occupied a central node and was connected by a single mutational step to all other alleles except for CYPᵀ-5, which differed by two mutational steps (Figure 3B). Additionally, our phylogenetic analysis recovered lack of structure in the backbone of the inferred tree but well-supported clades where individuals sharing the same CYPᵀ allele are grouped together (Figure 4). Specifically, CYPᵀ-2 was shared among homostyles from different populations (T03, T04, T07, T10, and D*11; Figure 4); CYPᵀ-3 was unique to population T05, and CYPᵀ-5 was unique to population T09. Finally, two CYPᵀ alleles were detected in each of three populations: CYPᵀ-4 and CYPᵀ-6 were found in population T06, CYPᵀ-2 and CYPᵀ-4 in population T10, and CYPᵀ-2 and CYPᵀ-7 in population D*11 (Figure 4).
3.3 Evolutionary consequences of transitions to homostyly
Mating system in homostyles - Outcrossing rates estimated from progeny arrays were more variable and significantly lower in homostylous than heterostylous families (single locus [ts]: mean ± SE; 0.24 ± 0.10 [homostylous] and 0.83 ± 0.05 [heterostylous]; multi-locus [tm]: 0.14 ± 0.06 and 0.98 ± 0.02, respectively; Wilcoxon signed-rank test W = 63, P = 0.002 and 0.001; Figure S1). As expected, the frequency of homostyles within populations was positively associated with population-level estimates of selfing rate (pseudo-R² = 0.447, P < 0.001; Table 1, Figures 5A). Seed germination rate was significantly lower for homostyles than heterostyles (0.11 ± 0.03 vs. 0.26 ± 0.05; W = 308, P = 0.04) resulting in an estimated inbreeding depression (δ) of 0.58.
Patterns of genetic diversity – Eighteen out of 22 populations (except for D03, D06, D08 and D10) had significant deviations from Hardy-Weinberg equilibrium due to a heterozygote deficit (Supplementary Table S4). Population-level estimates of allelic richness (Ar) varied from 1.58 to 3.78 (mean ± SE, 3.19 ± 0.09), expected heterozygosity (He) from 0.14 to 0.53 (0.42 ± 0.01), and observed heterozygosity (Ho) from 0.11 to 0.46 (0.37 ± 0.02). Inbreeding coefficient (Fis) ranged from -0.08 to 0.39 among populations (0.13 ± 0.02; Table 1) and had a significant positive relationship with the frequency of homostyly ( = 0.405, P = 0.03), but not with population size (P = 0.27; Figure 5B).
Global values of genetic differentiation among populations were 0.083 (95% CI: 0.074–0.093) and 0.164 (95% CI: 0.147–0.181) for Fst and G’st, respectively (Table S4), and were strongly correlated with each other (pseudo-R² = 0.974, P = 0.001). Pairwise population values of Fst and G’st are shown in Supplementary Material Table S5. As expected, pairwise homostyle frequencies showed a significant positive relationship with pairwise G’st (Mantel test: = 0.481, P < 0.001).
3.4 Genetic structure and gene flow
InStruct analyses recovered two clusters (Figure 6A) and showed that T08 differs from the remaining 21 populations. Identical results were found with Structure (results not shown). Accordingly, the DAPC analyses separated T08 from the rest of the populations (Figure 6B). Furthermore, we found no evidence of isolation by distance (Mantel test, P = 0.457). Together, these results suggest that genetic differentiation among populations is moderate, corroborating Fst results presented above. Using the frequency of private alleles, we estimated that the effective number of migrants per generation (Nm) was 1.96. Comparison of different demographic scenarios lent higher support to models with inter-population gene flow than to the isolation model. Among the models with varying migration rates (m), the one with m = 0.2 had the highest posterior probability, indicating moderate gene flow (Figure S2).