\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Different molecular changes underlie the same phenotypic transition:
origins and consequences of independent shifts to homostyly within
species}



\author[1]{Emiliano Mora-Carrera}%
\author[1]{Rebecca Stubbs}%
\author[1]{Barbara Keller}%
\author[1]{Étienne Léveillé-Bourret }%
\author[2]{Jurriaan de Vos}%
\author[1]{Peter  Szövényi}%
\author[1]{Elena Conti}%
\affil[1]{University of Zurich}%
\affil[2]{University of Basel}%


\vspace{-1em}



  
  \date{April 28, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
The repeated transition from outcrossing to selfing is a key topic in
evolutionary biology. However, the molecular basis of such shifts has
been rarely examined due to lack of knowledge of the genes controlling
these transitions. A classic example of mating system transition is the
repeated shift from heterostyly to homostyly. Occurring in 28 angiosperm
families, heterostyly is characterized by the reciprocal position of
male and female sexual organs in two (or three) distinct, usually
self-incompatible floral morphs. Conversely, homostyly is characterized
by a single, self-compatible floral morph with reduced separation of
male and female organs, facilitating selfing. Here, we investigate the
origins of homostyly in \emph{Primula vulgaris} and its \emph{}
microevolutionary consequences by integrating surveys of the frequency
of homostyles in natural populations, DNA sequence analyses of the gene
controlling the position of female sexual organs (\emph{CYP?}\selectlanguage{english}), and
microsatellite genotyping of both progeny arrays and natural populations
characterized by varying frequencies of homostyles. As expected, we
found that homostyles displace short-styled individuals, but long-style
morphs are maintained at low frequencies within populations. We also
demonstrated that homostyles repeatedly evolved from short-styled
individuals in association with different types of loss-of-function
mutations in \emph{CYP?}. \selectlanguage{english}Additionally, homostyly triggers a shift to
selfing, promoting increased inbreeding within and genetic
differentiation among populations. Our results elucidate the causes and
consequences of repeated transitions to homostyly within species,
enabling a likely explanation for the fact that homostyly has not become
fixed in \emph{P. vulgaris}. This study represents a benchmark for
future analyses of losses of heterostyly.%
\end{abstract}%



\sloppy




Article type: Original Article

\textbf{Different molecular changes underlie the same phenotypic
transition: origins and consequences of independent shifts to homostyly
within species}

Short title: Independent shifts to homostyly within species

Mora-Carrera E\textsuperscript{1*}, Stubbs RL\textsuperscript{1}, Keller
B\textsuperscript{1}, Leveille-Bourret E\textsuperscript{1,2}, de Vos
JM\textsuperscript{3}, Sz\selectlanguage{ngerman}övényi P\textsuperscript{1} and Conti
E\textsuperscript{1*}

\textsuperscript{1}Department of Systematic and Evolutionary Botany,
University of Zurich, Zollikerstrasse 107, Zurich 8008,
Switzerland;\textsuperscript{2}Institut de recherche en biologie
végétale, Département de sciences biologiques, Université de Montréal,
Québec, Canada;

\textsuperscript{3}Department of Environmental Sciences -- Botany,
University of Basel, Schönbeinstrasse 6, 4056 Basel, Switzerland;

*Corresponding
authors:\emph{emiliano.mora@systbot.uzh.ch},\emph{elena.conti@systbot.uzh.ch}

\textbf{ABSTRACT}

The repeated transition from outcrossing to selfing is a key topic in
evolutionary biology. However, the molecular basis of such shifts has
been rarely examined due to lack of knowledge of the genes controlling
these transitions. A classic example of mating system transition is the
repeated shift from heterostyly to homostyly. Occurring in 28 angiosperm
families, heterostyly is characterized by the reciprocal position of
male and female sexual organs in two (or three) distinct, usually
self-incompatible floral morphs. Conversely, homostyly is characterized
by a single, self-compatible floral morph with reduced separation of
male and female organs, facilitating selfing. Here, we investigate the
origins of homostyly in \emph{Primula vulgaris} and its \emph{}
microevolutionary consequences by integrating surveys of the frequency
of homostyles in natural populations, DNA sequence analyses of the gene
controlling the position of female sexual organs (\selectlanguage{english}\emph{CYP?}), and
microsatellite genotyping of both progeny arrays and natural populations
characterized by varying frequencies of homostyles. As expected, we
found that homostyles displace short-styled individuals, but long-style
morphs are maintained at low frequencies within populations. We also
demonstrated that homostyles repeatedly evolved from short-styled
individuals in association with different types of loss-of-function
mutations in \emph{CYP?}. \selectlanguage{english}Additionally, homostyly triggers a shift to
selfing, promoting increased inbreeding within and genetic
differentiation among populations. Our results elucidate the causes and
consequences of repeated transitions to homostyly within species,
enabling a likely explanation for the fact that homostyly has not become
fixed in \emph{P. vulgaris}. This study represents a benchmark for
future analyses of losses of heterostyly.

\textbf{KEYWORDS}

heterostyly, intra-specific, loss-of-function mutations, mating system,
\emph{Primula,} selfing

\textbf{1 INTRODUCTION}

{} A central topic in evolutionary biology is the repeated transition
from outcrossing to selfing (Cutter, 2019). A main advantage of selfers
vs. outcrossers is the increase in the number of gene copies transmitted
to the next generation (i.e., automatic selection advantage; Busch \&
Delph, 2012). Additionally, selfers can reproduce when mates or
pollinators, in the case of plants, are scarce (i.e., reproductive
assurance; Darwin, 1876). These advantages are counteracted by fitness
reduction in selfed progeny due to inbreeding depression (Schemske \&
Lande, 1985). In plants, transitions to selfing are enabled by the loss
of self-incompatibility (SI). The genetic basis of shifts from SI to
self-compatibility (SC) has been investigated in systems where the genes
responsible for self-pollen recognition have been characterized (Shimizu
\& Tsuchimatsu, 2015). Repeated transitions to SC \emph{via} the
inactivation of SI genes are often studied by comparing closely related
species with contrasting mating systems (Vekemans, Poux, Goubet, \&
Castric, 2014). Although independent losses of SI have also been
reported within species (Busch, Joly, \& Schoen, 2011; Foxe et al.,
2010; Goodwillie 1999), only a few studies have investigated the
molecular basis of these losses (Chantha, Herman, Platts, Vekemans, \&
Schoen, 2013; Shimizu, Shimizu-Inatsugi, Tsuchimatsu, \& Purugganan,
2008).

{} A classic model for the study of mating system transitions is the
shift from heterostyly to homostyly (Barrett, 2019). Heterostyly,
occurring in at least 28 angiosperm families, is a floral polymorphism
that promotes outcrossing, whereby two (distyly) or three (tristyly)
floral morphs differing in the position of reproductive organs co-occur
in a population (Darwin, 1877). In distylous species, flowers of the
long-styled morph, henceforth L-morph (also known as pin), exhibit the
stigma above the anthers, whereas flowers of the short-styled morph,
henceforth S-morph (also known as thrum), show the reciprocal
arrangement (Ganders, 1979; Figure 1). The reciprocal position of male
and female sexual organs in flowers of heterostylous taxa simultaneously
promotes pollen export to individuals of the opposite floral morph and
reduces self-pollination (Keller, Thomson, \& Conti, 2014; Lloyd \&
Webb, 1992). Moreover, heterostyly is often, but not always, associated
with a heteromorphic self-incompatibility system that reduces self- and
intra-morph fertilization (Dulberger, 1992). Heterostyly is controlled
by a single Mendelian locus known as the heterostyly supergene or
S-locus (Lewis \& Jones, 1992). Recent genomic and functional studies in
\emph{Primula} indicated that the S-locus is entirely absent from the
L-morph and hemizygous in the S-morph, where it consists of five tightly
linked genes (Huu et al., 2016; Li et al., 2016). Two genes in the
S-locus, GLOBOSA2 (\emph{GLO?}\selectlanguage{english}) and CYP734A50 (\emph{CYP?}\selectlanguage{english}), have been
experimentally demonstrated to control anther and stigma positions,
respectively (Huu et al., 2016; 2020).

{} As common for complex traits, heterostyly has been repeatedly lost.
Conforming to predictions about mating system transitions, losses of
heterostyly are usually associated with the decline or absence of
specialized pollinators (P\selectlanguage{ngerman}érez-Barrales \& Arroyo, 2010; Yuan et al.,
2017) or the absence of plant mates, as in colonization after
long-distance dispersal (Ness, Wright, \& Barrett, 2010) and
recolonization after glacial retreat (Guggisberg, Mansion, Kelso, \&
Conti, 2006). The breakdown of heterostyly can occur either through the
loss of one or more of the heterostylous floral morphs (e.g.,
transitions from tristyly and distyly to monomorphism; Ness et al.,
2010; Pérez-Barrales, Pino, Albaladejo, \& Arroyo, 2009) or the invasion
of homostylous plants (Barrett, 2019). Homostylous flowers are
characterized by reduced distance between male and female sexual organs
(i.e., reduced herkogamy) and, typically, loss of SI (de Vos et al.,
2014)\selectlanguage{english}. Both long- and short-homostyles are known, the latter being rare
(Lewis \& Jones, 1992). Long- and short-homostyles, respectively, bear
both anthers and stigmas either higher (H-morph in Figure 1) or lower in
the corolla tube. Homostyly enables reproductive assurance through
autonomous self-pollination, as demonstrated in \emph{Primula} (Carlson,
Gisler, \& Kelso, 2008; de Vos, Keller, Isham, Kelso, \& Conti, 2012;
Yuan et al., 2017)\selectlanguage{english}, and increases selfing rates (Husband \& Barrett,
1993; Schoen, Johnston, L'Heureux, \& Marsolais, 1997; Zhong et al.,
2019)\selectlanguage{english}. Repeated shifts from heterostyly to homostyly have been reported
both between (de Vos, Hughes, Schneeweiss, Moore, \& Conti, 2014;
Kissling \& Barrett, 2013; Kohn, Graham, Morton, Doyle, \& Barrett,
1996; Ruiz-Martin et al., 2018) and within species (Ness et al., 2010;
Shao et al., 2019; Zhou et al., 2017). However, none of these studies
identified the molecular basis of such transitions, likely because the
necessary genomic resources have only recently become available.

{} Recent studies in \emph{Primula} revealed that the transition from
heterostyly to homostyly is achieved \emph{via} disruption of either
\selectlanguage{english}\emph{GLO?} or \emph{CYP?}. \selectlanguage{english}Silencing \emph{GLO?} \selectlanguage{english}caused a reduction of
anther height in S-individuals, leading to the expression of
short-homostyly in \emph{Primula forbesii} (Huu, Keller, Conti, Kappel,
\& Lenhard, 2020). However, heteromorphic SI was retained, \emph{}
preventing self-fertilization and possibly explaining, together with
the lower mutation rate of \emph{GLO?} \selectlanguage{english}than \emph{CYP?}, \selectlanguage{english}the low
frequency of short-homostyles observed in nature (Huu et al., 2020).
Similarly, a transposable element insertion likely compromising the
function of \emph{GLO?} \selectlanguage{english}was detected in one short-homostyle of a
cultivated \emph{P. vulgaris} variety (Li et al., 2016). Conversely,
silencing \emph{CYP?} \selectlanguage{english}in S-individuals of \emph{P. veris} increased
stigma height, causing the shift to long-homostyly (Huu et al., 2016).
Correspondingly, the expression of \emph{CYP?} \selectlanguage{english}was reduced in
long-homostyles compared to S-individuals of both \emph{P. vulgaris} and
\emph{P. oreodoxa} (Huu et al., 2016; Zhao, Luo, Yuan, Mei, \& Zhang,
2019). Furthermore, all five exons of \emph{CYP?} \selectlanguage{english}were likely lost in
long-homostyles of the ancestrally heterostylous \emph{P. forbesii},
whereas exon three of \emph{CYP?} \selectlanguage{english}was not detected by PCR assays in the
long-homostylous \emph{P. grandis, P. halleri,} and \emph{P. scotica}
(Huu et al., 2016). Moreover, two loss-of-function mutations in
\emph{CYP?} \selectlanguage{english}were reported in self-compatible long-homostyles of \emph{P.
vulgaris,} suggesting that \emph{CYP?} \selectlanguage{english}might additionally be involved in
the control of SI (Li et al., 2016), but this has not been
experimentally demonstrated \emph{} (Kappel, Huu, \& Lenhard, 2017). The
findings above demonstrate that any genetic changes disrupting the
function or expression of \emph{CYP?}, \selectlanguage{english}a gene present only in
S-individuals, cause long-homostyly. However, comprehensive analyses of
the different types of mutations associated with shifts to homostyly
based on extensive sampling of heterostyles and homostyles from multiple
natural populations have never been performed.

{} After homostyles originate, their reproductive advantages over
heterostyles should prompt changes in floral morph composition and
genetic diversity of populations. Previous theoretical models predicted
that homostyles should increase in frequency by replacing S-individuals
(Crosby, 1949; 1960). Concomitantly, occasional crosses between
homostyles and L-individuals and reduced fitness of homostyles,
presumably caused by homozygosity of genes at the S-locus, were expected
to prevent the fixation of homostyly (Crosby, 1949; 1960). Both
mechanisms should ultimately result in the co-occurrence of homostyles
with L-individuals, the latter remaining at low frequency. The
predictions above were supported by two surveys of the frequency of L-,
S- and homostylous individuals from natural populations of \emph{P.
vulgaris} performed at a 36-yr. interval (Crosby, 1949; Curtis \&
Curtis, 1985). Subsequent studies testing whether homostyles were
affected by inbreeding depression found no evidence of seed-mass
differences between heterostyles and homostyles (Piper, Charlesworth, \&
Charlesworth, 1986). However, inbreeding depression at subsequent
life-cycle stages has never been investigated. At the population genetic
level, higher selfing in homostyles should increase homozygosity, thus
decreasing diversity within and increasing differentiation among
populations (Hamrick \& Godt, 1996). These microevolutionary
consequences of mating system transitions have been quantified only in a
few homostylous species (Ness et al., 2010; Yuan et al., 2017; Zhou et
al., 2017; Zhong et al., 2019).

\emph{Primula vulgaris} Huds. (the common primrose) represents a
classical model for the study of heterostyly and intraspecific
transitions to homostyly (Bodmer, 1958; Cocker et al., 2018; Crosby,
1959; Piper et al., 1986). While distylous populations consisting
exclusively of L- and S-individuals are common throughout Eurasia, mixed
populations with varying frequencies of long-homostylous, S-, and
L-individuals have been reported and widely documented for decades only
in Somerset, England (Crosby, 1940; Curtis \& Curtis, 1985). The
occurrence of homostyles in multiple populations and the wide variation
in their frequency in this region prompted Crosby (1949) to suggest that
homostyles migrated from an initial center of origin to neighboring
populations. However, two different loss-of-function mutations in
\emph{CYP?} \selectlanguage{english}were recently detected in two long-homostyles from different
natural populations of \emph{P. vulgaris} in England (Li et al., 2016),
suggesting that homostyly might arise independently rather than
\emph{via} migration. Therefore, population genetic analyses aimed at
clarifying the role of migration in the spread of homostyly are needed.

{} This study investigates the origins of homostyly, its effects on the
mating system, and its population genetic consequences using the
well-characterized \emph{P. vulgaris} model. We thus combined surveys of
floral morph frequencies from multiple natural populations of the
ancestrally heterostylous \emph{P. vulgaris} with both DNA sequence
analyses of \emph{CYP?} \selectlanguage{english}and microsatellite data from multiple
individuals to ask the following questions: 1) Is an increase in the
frequency of homostyles associated with a reduction of S-individuals and
the maintenance of L-individuals at low frequencies, as predicted by
Crosby's (1949) model and found in earlier studies (Curtis \& Curtis,
1985)? 2) Is the transition to homostyly associated with a single or
multiple molecular changes in \emph{CYP?} \selectlanguage{english}and which ones? 3) Is the loss
of heterostyly accompanied by a shift from outcrossing to selfing? And,
if so, is there evidence of inbreeding depression as well as increased
inbreeding within and genetic differentiation among populations? Our
study generates new knowledge on the molecular basis of transitions from
outcrossing to selfing, its evolutionary consequences, and the putative
mechanisms precluding the fixation of homostyly within species.

\textbf{2 METHODS}

\emph{} \emph{2.1 Study species}

{} The ancestrally heterostylous \emph{Primula vulgaris} is a perennial,
rosette-forming diploid (2n = 22) blooming in early spring
(February--April) with pale-yellow fused corollas reaching up to 40 mm
in width formed by broad, overlapping, v-notched lobes (Richards, 2003).
\emph{Primula vulgaris} has dimorphic pollen grains and stigma papillae
(Mast, Kelso, \& Conti, 2006). Pollen grains are small in L-individuals
but large in S-individuals, whereas stigma papillae are long in
L-individuals but short in S-individuals (Richards, 2003). The species
reproduces mainly sexually; clonal reproduction, uncommon in the wild,
has been reported under cultivation (Boyd, Silvertown, \& Tucker,
1990). Pollination in natural populations of \emph{P. vulgaris} is
mostly performed by long-tongued bees, bumblebees, and bee-flies, but
butterflies, pollen-gathering bees and moths have been occasionally
reported (Boyd et al., 1990; Keller et al., 2020; Woodell, 1960).
Seeds, produced in capsules, are dispersed by ants and rodents (Valverde
\& Silvertown, 1995). \emph{Primula vulgaris} has a patchy distribution
with populations ranging from dozens to hundreds of individuals mostly
found in woodlands and hedgerows and less frequently in grasslands
(Jacquemyn, Endels, Brys, Hermy, \& Woodell, 2009). Previous studies
indicated that long-homostyles of \emph{P. vulgaris} are self-compatible
and produce seeds in the absence of pollinators (Piper et al., 1986).
Analyses of the mating system of homostyles estimated a range of
outcrossing rates from 0.05 to 0.8 (Bodmer, 1958; Crosby, 1959; Piper,
Charlesworth, \& Charlesworth, 1984). It has been suggested that
reduced pollinator service due to intensive pastoral and agricultural
activities could be the selective pressure favoring the increase of
self-fertilizing homostyles in Somerset, England (Piper et al., 1986).
Recently, long-homostyles have been reported in one population of an
agriculturally fragmented area in the Netherlands (Barmentlo, Meirmans,
Luijten, Triest, \& Oostermeijer, 2017).

\emph{2.2 Frequency variation of floral morphs in natural populations
of} Primula vulgaris

{} To test whether homostyles increase concomitantly with the decrease
of S-individuals, while L-individuals are maintained at low frequencies,
we quantified floral morph composition in 22 natural populations of the
Somerset region (England) in the spring of 2019 (Figure 2A). The
approximate location of these populations was based on Curtis and Curtis
(1985). The frequency of floral morphs was estimated by randomly
scoring 100 individuals whenever possible (Table 1). Flowering plants
were visually assigned to S-, L- and H-morphs (Figure 1). To determine
whether S- and L-morphs segregated at equal frequencies in both
dimorphic heterostylous (i.e., populations consisting of S- and
L-individuals) and trimorphic populations (i.e., populations consisting
of S-, L- and H-individuals), we tested whether morph ratios deviated
significantly from isoplethy (i.e., equal frequencies of S- and
L-morphs) using G-tests of goodness-of-fit as implemented in the R
package `DescTools' v0.99.38 (Signorell et al., 2020). We estimated a
ternary plot using the R package `ggtern' v3.1.0 (Hamilton \& Ferry,
2018) to assess whether the increase in the frequency of H-morphs
occurred at the cost of only one or both S- and L-morphs.

\emph{2.3 Sampling for genetic analyses}

In the spring of 2019, we randomly collected leaf tissue from 11 to 30
individuals in each of the 22 natural populations described above for a
total of 613 individuals. Leaf tissue was dried in silica gel. We
sampled the same number of each floral morph where possible; sampled
individuals were at least 1--2 meters apart from each other. We
additionally collected leaf tissue and 2--3 capsules with seeds from
nine individuals (four pins and five thrums) of one dimorphic
heterostylous population (D07) and seven individuals from the fully
homostylous population (M01). Collected seeds from each individual were
sowed in a single pot; seedlings from a single individual formed a
family. All pots were placed in water containing gibberellic acid (70
ppm) for 24 hours and later placed in growth chambers at 18 \selectlanguage{ngerman}°C with a
12:12 h light-dark period. The proportion of germinated seeds per family
was recorded after eight weeks, and leaf tissue from a total of ten
seedlings per family was collected whenever possible and dried in silica
gel. DNA of individuals from natural populations and progeny arrays was
extracted with a modified CTAB protocol (Doyle \& Doyle, 1987) and used
for subsequent genetic analyses.

\emph{2.4 Molecular basis of the transition to homostyly}

To determine whether the transition to homostyly is associated with a
single or multiple molecular changes in \selectlanguage{english}\emph{CYP?} and which ones, we
sequenced \emph{CYP?} \selectlanguage{english}from multiple S-individuals and homostyles of
\emph{P. vulgaris} from different populations.

{} \emph{Data generation} - We Sanger-sequenced all five exons of
\emph{CYP?} \selectlanguage{english}from 38 individuals (27 H-individuals from 10 populations
and 11 S-individuals from five populations, corresponding to 2--3
individuals per population). Additionally, we obtained \emph{CYP?}
\selectlanguage{english}sequences \emph{} from whole genome resequencing data of six
S-individuals from two populations (three individuals per population;
EMC, unpublished data). Finally, we complemented our data set with
previously published sequences from two H-individuals (\emph{CYP?
S}\selectlanguage{english}\textsuperscript{\emph{LH1}} and
\emph{S}\textsuperscript{\emph{LH2}}\emph{,} respectively, of Li et al.,
\emph{} 2016). In total, we thus analyzed \emph{CYP?} \selectlanguage{english}sequences from 46
individuals (29 H- and 17 S-morphs).

{} For Sanger sequencing\emph{,} we used previously designed primers
based on all five \emph{CYP?} \selectlanguage{english}exons (Huu et al., 2016). Reverse primers
for exons 1, 2, and 3 were newly designed to obtain longer sequences
(Supplementary material Table S1). Detailed PCR conditions for the
amplification of all \emph{CYP?} \selectlanguage{english}exons are included in Supplementary
Material. Sanger sequencing was performed in an ABI Prism 3130 genetic
analyzer (Applied Biosystems). Forward and reverse sequences were
visually inspected and aligned with MUSCLE, as implemented in MEGA X
(Kumar, Stecher, Li, Knyaz, \& Tamura, 2018). All exon sequences were
concatenated with Mesquite v3.61 (Maddison \& Maddison, 2019). Library
preparation and high throughput sequencing for the whole genome
resequencing data were performed by RAPiD GENOMICS (Gainesville,
Florida, USA) using paired-end sequencing (\textasciitilde{}150bp
sequence reads) in NovaSeq 6000 (Illumina). We used HybPiper v1.3.1
(Johnson et al., 2016) with default parameters (except for the
coverage-cutoff level for assemblies set to 4) to target and extract the
sequence of \emph{CYP?} \selectlanguage{english}from whole genome resequencing data. To identify
synonymous and nonsynonymous mutations, we used the Open Reading Frame
(ORF) from \emph{P. vulgaris} \emph{CYP?} \selectlanguage{english}sequence deposited in GenBank
(KT257665.1; Li et al., 2016).

{} \emph{Analyses} -- To determine the extent of \emph{CYP?} \selectlanguage{english}variation
in S-individuals and homostyles, we estimated nucleotide diversity at
synonymous (\selectlanguage{greek}π\selectlanguage{english}S) and nonsynonymous sites (\selectlanguage{greek}π\selectlanguage{english}NS) using DNAsp v6.12.03
(Rozas et al., 2017). To determine the relationships among \emph{CYP?}
\selectlanguage{english}sequences from H- and S-individuals, we first estimated a haplotype
network with the R package `pegas' v0.12 (Paradis, 2010). Such
analysis is appropriate because \emph{CYP?} \selectlanguage{english}is located in the
hemizygous, non-recombining S-locus. In addition to single nucleotide
substitutions, we scored each insertion or deletion (i.e., indel) as a
character following the ``simple indel contig'' guidelines by Simmons
and Ochoterena (2000). Briefly, each indel of any length was scored as
a presence/absence character in every individual. The resulting
presence/absence matrix was coded \emph{} as pseudo-nucleotides using
`A' as absence, `T' as presence and `-' as unknown, and was concatenated
to the sequence alignment of \emph{CYP?}. \selectlanguage{english}Secondly, we estimated the
\emph{CYP?} \selectlanguage{english}phylogeny using a partitioned Maximum Likelihood (ML)
analysis in RAxML v8.01 (Kozlov, Darriba, Flouri, Morel, \& Stamatakis,
2019) with a GTR-GAMMA substitution model for nucleotide substitutions
and a binary (BIN) model for indels. Branch support was estimated with
1000 standard bootstrap re-samplings. A publicly available \emph{CYP?}
\selectlanguage{english}sequence from \emph{P. veris} (KX589238; Huu et al., 2016) was used as
an outgroup to root the \emph{CYP?} \selectlanguage{english}tree.

\emph{2.5 Evolutionary consequences of transitions to homostyly}

We \emph{} tested whether outcrossing rates from progeny arrays are
higher in heterostyles than homostyles with microsatellites. We
estimated the strength of inbreeding depression at the seed-germination
stage. We tested whether increasing frequencies of homostyles are
associated with higher population-level selfing rates, inbreeding
coefficients within populations, and genetic differentiation among
populations.

\emph{2.5.1 Outcrossing rates from progeny arrays}

{} Outcrossing rates were estimated from nine heterostylous families
(stemming from the dimorphic heterostylous population D07) with a mean
number of 10 seedlings per family and seven homostylous families
(stemming from the fully homostylous population M01) with 7.86 \selectlanguage{ngerman}± 1.06
(mean ± SE) seedlings per family. Both maternal and seedling plants were
genotyped with the same microsatellites as in section 2.5.2. Multi-locus
(\emph{t}\textsubscript{\emph{m}}) and single-locus
(\emph{t}\textsubscript{\emph{s}}) outcrossing rates and their
associated standard errors were estimated with MLTR v.3.4 (Ritland,
2002)\selectlanguage{english}\selectlanguage{english}, using 10000 bootstrap resamplings of maternal families. We used
a one-sided Wilcoxon signed-rank test to assess whether homostylous
individuals have a lower outcrossing rate than distylous individuals.
Moreover, the germination rate of heterostylous and homostylous families
was used as a fitness proxy to estimate the strength of inbreeding
depression (\selectlanguage{greek}δ\selectlanguage{english}) as 1 -- (\emph{W}\textsubscript{\emph{HO}} /
\emph{W}\textsubscript{\emph{HE}}) (Goodwillie, Kalisz, \& Eckert,
2005) where \emph{W}\textsubscript{\emph{HO}} \textsubscript{} and
\emph{W}\textsubscript{\emph{HE}} are the mean germination rates of
homostyles and heterostyles, respectively.

\emph{2.5.2 Population-level estimates of selfing rates}

{} \emph{Data generation -} We genotyped all 613 individuals collected
from the 22 natural populations using 12 microsatellites previously
developed for \emph{P. vulgaris} (Triest, Sierens, \& Van Rossum,
2015). Briefly, the loci were amplified in two multiplex reactions with
a QIAGEN Multiplex PCR kit. Detailed conditions for the amplification of
microsatellites are included in Supplementary Material. Fragment size
estimations were performed using an ABI Prism 3130xl genetic analyzer
(Applied Biosystems) with GeneScan 500 LIZ (Applied Biosystems) as
internal standard. Alleles were scored using GeneMapper v4.1 (Applied
Biosystems). The presence of null alleles was estimated with
Micro-checker v2.2.3 (Van Oosterhout, Hutchinson, Wills, \& Shipley,
2004).

{} \emph{Analyses} -- Population-level selfing rates were inferred using
model-based Bayesian estimation (BES v0.1.3; Redelings et al., 2015)
from microsatellite data of 613 individuals from 22 natural populations
(see section 2.3). The main advantage of this method compared to
estimates generated by other programs such as RMES (David, Pujol, Viard,
\& Goudet, 2007) is that BES makes use of all loci regardless of
whether they are homozygous or heterozygous. Selfing rates and 95\% CIs
were inferred using 10000 bootstrap iterations. To test whether higher
frequencies of homostyles in natural populations are associated with
higher selfing rates, we used generalized additive models where the
frequency of homostyles was used as an independent variable and by
specifying a beta distribution for the estimates of selfing rates as
implemented in the R package `gamlss' \emph{} v5.1.6 (Rigby,
Stasinopoulos, \& Lane, 2005).

\emph{2.5.3 Population genetic analyses}

{} Allelic richness (\emph{Ar}), observed
(\emph{H}\textsubscript{\emph{o}}) and expected heterozygosity
(\emph{H}\textsubscript{\emph{e}}), and inbreeding coefficient
(\emph{F}\textsubscript{\emph{is}}) were calculated with the R package
`diveRsity' v9.90 (Keenan, Mcginnity, Cross, Crozier, \& Prod\selectlanguage{ngerman}öhl, 2013)\selectlanguage{english}
from microsatellite data of 613 individuals from 22 populations (see
section 2.3)\selectlanguage{english}. Each locus and population were tested for Hardy-Weinberg
equilibrium using the R package `genepop' v1.1.7 (Rousset, 2008).\selectlanguage{english} To
assess whether an increase in the frequency of homostyly is associated
with an increase in \emph{F}\textsubscript{\emph{is}}\emph{,} we used a
standard linear model with homostyle frequency as the independent
variable. Given that population size can influence
\emph{F}\textsubscript{\emph{is}}\emph{,} we included census size as
covariable with no interaction among independent variables.

{} Pairwise genetic differentiation measured as
\emph{F}\textsubscript{\emph{st}} (Weir \& Cockerham, 1984)\selectlanguage{english}\selectlanguage{english} and its
normalized counterpart \emph{G'}\textsubscript{\emph{st}} (Meirmans \&
Hedrick, 2011)\selectlanguage{english} were estimated with the R package `diveRsity' v9.90
(Keenan et al., 2013). To discover whether higher homostyle frequencies
increase genetic differentiation, we tested whether pairwise differences
in homostyle frequencies between populations are positively correlated
with pairwise \emph{G'}\textsubscript{\emph{st}} \emph{} values. To test
the significance of association between these two matrices, we used a
Mantel test as implemented in the R package `adegenet' v2.1.3 (Jombart,
2008)\selectlanguage{english}.

\emph{2.6 Genetic structure and gene flow}

{} To estimate genetic structure among \emph{P. vulgaris} populations,
we used InStruct v3.2.09 (Gao, Williamson, \& Bustamante, 2007)\selectlanguage{english}, which
extends the Bayesian model implemented in Structure (Pritchard,
Stephens, \& Donnelly, 2000) \selectlanguage{english}to situations of partial
self-fertilization and recurrent inbreeding. We ran the analyses with
100000 burn-in iterations and 20 independent MCMC chains per run from 2
to 22 \emph{k}, as suggested by Gilbert et al. (2012)\selectlanguage{english}. The most likely
number of clusters (\emph{k}) was detected with the \selectlanguage{greek}\emph{Δk}\selectlanguage{english} method
(Evanno, Regnaut, \& Goudet, 2005) implemented in CLUMPAK (Kopelman,
Mayzel, Jakobsson, Rosenberg, \& Mayrose, 2015). Additionally, a
Discriminant Analysis of Principal Components (DAPC) was performed with
the R package `adegenet' v2.1.1 (Jombart, 2008). To determine whether
there is Isolation By Distance (\emph{IBD}), we performed a Mantel test
implemented in `adegenet' \emph{} using both
\emph{F}\textsubscript{\emph{st}} and \emph{G'}\textsubscript{\emph{st}}
values of genetic differentiation.

{} Finally, to determine whether patterns of genetic structure are
explained by gene flow, we estimated the size-corrected effective number
of migrants per generation (\emph{Nm}) using the frequency of private
alleles with the R package `genepop' v1.1.7 (Rousset, 2008). Moreover,
we used Approximate Bayesian Computations (ABC), as implemented in
popABC v1.0 (Lopes, Balding, \& Beaumont, 2009), to compare models of
genetic divergence with and without gene flow using the coalescent
framework developed by Nielsen and Wakely (2001) and Hey and Nielsen
(2004). Specifically, we compared a scenario of isolation without gene
flow (i.e., proportion of migrants per generation, \emph{m} = 0) with
seven scenarios of increasing migration rates (\emph{m} = 0.001, 0.01,
0.1, 0.2, 0.3, 0.4 and 0.5). Priors used for this analysis are specified
in Supplementary Material (Table S2). Each model was simulated with
100000 iterations and tolerance for the rejection step was set to 0.01.
Model selection was based on the posterior probability estimated by
categorical regression (Beaumont, Zhang, \& Balding, 2002) using custom
R scripts included in the popABC documentation. A major assumption of
the models implemented in popABC is that there is no genetic structure
within populations. Given that H-individuals may cause within-population
structure due to selfing, we decided to use only the ten dimorphic
heterostylous populations for popABC analyses because no genetic
substructure is expected in populations consisting entirely of
outcrossing individuals. No differences in pollen and seed dispersal are
known for H-, S- and L-individuals, thus excluding populations with
H-individuals from popABC analyses should not affect estimates of
migration between populations.

\textbf{3 RESULTS}

\emph{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 \emph{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 (\emph{G}
= 7.13, \emph{P} = 0.007 and \emph{G} = 6.11, \emph{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 (\emph{G} = 0.13, \emph{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).

\emph{3.2 Molecular basis of the transition to homostyly}

{} We sequenced \textasciitilde{}97\% of the 1,587 bp comprising the
five exons of \emph{CYP?} \selectlanguage{english}from 17 S- and 27 H-individuals of \emph{P.
vulgaris} from seven and ten populations, respectively. In total, we
detected seven \emph{CYP?} \selectlanguage{english}alleles (\emph{CYP?-1}\selectlanguage{english} to \emph{-7}; Figure
3A). All 17 S-individuals shared the same allele (\emph{CYP?-1}\selectlanguage{english}),
identical to the functional copy of \emph{CYP?} \selectlanguage{english}previously reported in
\emph{P. vulgaris} (GenBank: KT257665.1; Li et al., 2016). Unexpectedly,
the putatively functional \emph{CYP?-1}\selectlanguage{english} allele also occurred in six
H-individuals. In the 21 remaining H-individuals, we detected six
\emph{CYP?} \selectlanguage{english}alleles (\emph{CYP?-2}\selectlanguage{english} to \emph{-7}), comprising altogether
a single synonymous mutation, four nonsynonymous mutations, and two
deletions. Specifically, \emph{CYP?-4}\selectlanguage{english} contained two private mutations
(one synonymous and one nonsynonymous) and five \emph{CYP?} \selectlanguage{english}alleles
contained a single private mutation \emph{} (three different
nonsynonymous mutations in \emph{CYP?-2, -3, -7}\selectlanguage{english} and two different
deletions \emph{} in \emph{CYP?-5}\selectlanguage{english} and \emph{-6})\emph{.} The
nonsynonymous mutation in \emph{CYP?-2}\selectlanguage{english} introduced a premature stop
codon in exon 2, whereas a 31 bp deletion in exon 5 of \emph{CYP?-5}\selectlanguage{english} and
an 8 bp deletion in exon 1 of \emph{CYP?-6}\selectlanguage{english} introduced a frameshift
leading to an early stop codon\emph{.} Moreover, a mutation in exon 5 of
\emph{CYP?-3}\selectlanguage{english} changed a non-polar (phenylalanine) to a polar amino acid
(serine), whereas a mutation in exon 3 of \emph{CYP?-4}\selectlanguage{english} changed an
amidic (asparagine) to a hydroxylic (serine) amino acid. Finally, the
arginine to histidine mutation in \emph{CYP?-7}\selectlanguage{english} caused no change in
amino acid side-chain polarity. Overall, nucleotide diversity in
homostyles was \selectlanguage{greek}π\selectlanguage{english}S = 0.0005 and \selectlanguage{greek}π\selectlanguage{english}NS = 0.001, while S-morph individuals
showed no genetic variation in \emph{CYP?} \selectlanguage{english}(both \selectlanguage{greek}π\selectlanguage{english}S and \selectlanguage{greek}π\selectlanguage{english}NS were zero).

{} In the haplotype network, \emph{CYP?-1}\selectlanguage{english}, the known functional allele,
occupied a central node and was connected by a single mutational step to
all other alleles except for \emph{CYP?-5,}\selectlanguage{english} which differed by two
mutational steps (Figure 3B)\emph{.} Additionally, our phylogenetic
analysis recovered lack of structure in the backbone of the inferred
tree but well-supported clades where individuals sharing the same
\emph{CYP?} \selectlanguage{english}allele are grouped together (Figure 4). Specifically,
\emph{CYP?-2}\selectlanguage{english} was shared among homostyles from different populations
(T03, T04, T07, T10, and D*11; Figure 4); \emph{CYP?-3}\selectlanguage{english} was unique to
population T05, and \emph{CYP?-5}\selectlanguage{english} was unique to population T09. Finally,
two \emph{CYP?} \selectlanguage{english}alleles were detected in each of three populations:
\emph{CYP?-4}\selectlanguage{english} and \emph{CYP?-6}\selectlanguage{english} were found in population T06,
\emph{CYP?-2}\selectlanguage{english} and \emph{CYP?-4}\selectlanguage{english} in population T10, and \emph{CYP?-2}\selectlanguage{english} and
\emph{CYP?-7}\selectlanguage{english} in population D*11 (Figure 4).

\emph{3.3 Evolutionary consequences of transitions to homostyly}

\emph{Mating system in homostyles -} Outcrossing rates estimated from
progeny arrays were more variable and significantly lower in homostylous
than heterostylous families (single locus
{[}\emph{t}\textsubscript{s}{]}: mean \selectlanguage{ngerman}± SE; 0.24 ± 0.10
{[}homostylous{]} and 0.83 ± 0.05 {[}heterostylous{]}; multi-locus
{[}\emph{t}\textsubscript{\emph{m}}{]}: 0.14 ± 0.06 and 0.98 ± 0.02,
respectively; Wilcoxon signed-rank test \emph{W} = 63, \emph{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 (\emph{pseudo-R²} = 0.447, \emph{P} \textless{} 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; \emph{W} =
308, \emph{P} = 0.04) resulting in an estimated inbreeding depression
(\selectlanguage{greek}δ\selectlanguage{english}) of 0.58.

\emph{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 (\emph{Ar})
varied from 1.58 to 3.78 (mean \selectlanguage{ngerman}± SE, 3.19 ± 0.09), expected
heterozygosity (\emph{H}\textsubscript{\emph{e}}) from 0.14 to 0.53
(0.42 ± 0.01), and observed heterozygosity
(\emph{H}\textsubscript{\emph{o}}) from 0.11 to 0.46 (0.37 ± 0.02).
Inbreeding coefficient (\emph{F}\textsubscript{\emph{is}}) \emph{}
ranged \emph{} 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 (\emph{R²} = 0.405, \emph{P} = 0.03), but not with population
size (\emph{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
\emph{F}\textsubscript{\emph{st}} and
\emph{G'}\textsubscript{\emph{st}}, respectively (Table S4), and were
strongly correlated with each other (\emph{pseudo-R²} = 0.974, \emph{P}
= 0.001). Pairwise population values of
\emph{F}\textsubscript{\emph{st}} and \emph{G'}\textsubscript{\emph{st}}
are shown in Supplementary Material Table S5. As expected, pairwise
homostyle frequencies showed a significant positive relationship with
pairwise \emph{G'}\textsubscript{\emph{st}} (Mantel test: \emph{R²} =
0.481, \emph{P} \textless{} 0.001).

\emph{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, \emph{P} =
0.457). Together, these results suggest that genetic differentiation
among populations is moderate, corroborating
\emph{F}\textsubscript{\emph{st}} \emph{} results presented above. Using
the frequency of private alleles, we estimated that the effective number
of migrants per generation (\emph{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 (\emph{m}), the one with \emph{m} = 0.2 had
the highest posterior probability, indicating moderate gene flow (Figure
S2).

\textbf{4 DISCUSSION}

{} The present study of shifts from heterostyly to homostyly provides
novel insights into the molecular basis of intraspecific mating system
transitions and their evolutionary consequences. Previous studies
suggesting independent origins of homostyly used molecular markers not
associated with the genetic control of heterostyly (Ness et al., 2010;
Shao et al., 2019; Yuan et al., 2017; Zhou et al., 2017). By leveraging
new knowledge on the molecular basis of heterostyly and extensive
sampling of populations with different frequencies of homostyles, this
is the first study to demonstrate that heterostyly has been repeatedly
lost likely \emph{via} multiple, presumably loss-of-function mutations
of the S-locus gene controlling stigma position (\selectlanguage{english}\emph{CYP?}).
Additionally, we show that the transition to homostyly promotes a change
of mating system from outcrossing to selfing, with ensuing population
genetic consequences. Furthermore, we discuss the potential roles of
gene flow in the spreading of homostyly. Finally, we propose that
inbreeding depression might explain why homostyly, despite its apparent
advantages, has not become fixed in \emph{P. vulgaris}. Thus, the
present contribution provides a benchmark for similar studies on the
loss of heterostyly across angiosperms.

\emph{\textbf{4.1 Homostyles replace short-styled individuals within
populations}}

The variation in the frequency of floral morphs among populations
provides an opportunity to elucidate how homostylous phenotypes
originate and spread. According to previous evolutionary models, the
origin of homostyles from S-individuals should cause an initial
reduction of the latter, triggering an increase in the number of
homostyles (Crosby, 1949; 1959). When homostyles first arise, their
pollen can be used for both self- and cross-fertilization of
L-individuals, but not for fertilization of S-individuals, due to lack
of spatial reciprocity between high anthers of homostyles and low stigma
of S-individuals and incompatibility between homostyle pollen and
S-morph stigmas (Crosby, 1949; Figure 1). Furthermore, homostyles can be
fertilized by compatible pollen of S-morphs, but this cross is thought
to occur rarely in nature due to stigma clogging with self-pollen in
homostyles (Crosby, 1959; but see Bodmer, 1959). Hence, homostyles have
a reproductive advantage over L- and S-morphs stemming from their
ability to fertilize both themselves and the L-morphs (Richards, 2003).
Finally, self-fertilization of homostyles (Figures S3A and B) and
occasional crosses from homostyles to L-individuals produce homostylous
and L-progeny, but no S-progeny (Figures S3C and D).

{} Congruently with Crosby's predictions (1949) and findings from
previous studies (Curtis \& Curtis, 1985), our results demonstrate
that, as homostyles increase, more S- than L-individuals are lost within
populations (Table 1 and Figures 2A and B). After 72 and 36 years since
Crosby's (1949) and Curtis and Curtis (1985) observations, respectively,
we also found that homostyly has not become fixed in Somerset
populations. Furthermore, the fact that L-individuals are maintained at
low frequencies, even in the absence of S-individuals, suggests that
they are occasionally fertilized by homostyles. Experiments aimed at
quantifying pollen delivery and receipt among the three floral morphs
would help to clarify the selective regimes favoring the increase of
homostyles. Overall, our results are consistent with the hypothesis that
long-homostyles originate from S-individuals and explain why the latter
occur at lower frequencies than L-individuals across populations with
all three morphs.

\emph{\textbf{4.2 Multiple losses of heterostyly within species are
associated with different molecular changes in the CYP? gene}}\selectlanguage{english}

{} Recent advancements in genomics have enabled the identification of
loci responsible for repeated losses of phenotypic traits (Albalat \&
Ca\selectlanguage{ngerman}ñestro, 2016; Sharma et al., 2018). For instance, multiple
loss-of-function mutations were detected in \emph{SCR} or
\emph{SCR}-like genes corresponding to independent losses of
self-incompatibility in \emph{Arabidopsis thaliana} and
\emph{Laevenworthia alabamica} (Chantha et al., 2013; Shimizu et al.,
2008)\selectlanguage{english}. In the case of homostyly, recent experimental studies
established that \selectlanguage{english}\emph{CYP?} disruption causes a shift to homostyly (Huu
et al., 2016; Kappel et al., 2017). Accordingly, all S-individuals of
\emph{P. vulgaris} analyzed in the present study share the same
functional \emph{CYP?} \selectlanguage{english}allele (\emph{CYP?-1}\selectlanguage{english}), whereas private alleles
with nonsynonymous mutations or indels (\emph{CYP?-2}\selectlanguage{english} to -\emph{9}) are
found exclusively in homostyles (Figures 3A, 3B and 4). These mutations
cause either changes in the polarity of the amino acid side chain
(\emph{CYP?-3}\selectlanguage{english} and \emph{-4}) or premature stop codons (\emph{CYP?-2}\selectlanguage{english},
\emph{-5} and \emph{-6}), both types of mutation presumably inducing
loss or reduction of function (Zhang, 2000). In one case
(\emph{CYP?-7}\selectlanguage{english}), a nonsynonymous mutation did not change amino acid
side-chain polarity; thus, this mutation may or may not affect protein
function (Figure 3A). Notably, the alleles found in our study differ
from the ones previously reported in two long-homostyles of \emph{P.
vulgaris} (\emph{CYP?-8}\selectlanguage{english} and -\emph{9}, originally named \emph{CYP?
S}\selectlanguage{english}\textsuperscript{\emph{LH1}} \textsuperscript{} and
\emph{S}\textsuperscript{\emph{LH2}}, respectively; Li et al., 2016).
Thus, a total of eight different putative loss-of-function mutations in
\emph{CYP?} \selectlanguage{english}have been identified so far in homostylous individuals of
\emph{P. vulgaris}. Future work should focus on elucidating how the
found nonsynonymous mutations affect \emph{CYP?} \selectlanguage{english}function\emph{.}
Finally, the groupings of the haplotype network were congruent with
those of the phylogeny (Figures 3B and 4), as expected given the lack of
recombination in the S-locus, hence \emph{CYP?.} \selectlanguage{english}Overall, our results
show that a diversity of changes in the same gene is associated with
multiple, independent origins of homostyly within species.

{} Loss of phenotypic traits can be caused by mutations not only in
coding regions of specific genes but also in their promoter regions or
by structural rearrangements. For instance, some individuals of
self-compatible \emph{A. thaliana} were characterized by mutations in
\emph{cis}-regulatory regions, inversions, or splicing variants, causing
loss of function in \emph{SCR} or \emph{SRK} genes controlling
self-incompatibility (Dwyer et al., 2013; Shimizu et al., 2008; Shimizu
\& Tsuchimatsu, 2015). In \emph{P. vulgaris}, a surprising result is
that six out of 27 homostyles have the same \emph{CYP?} \selectlanguage{english}allele as the
S-individuals (\emph{CYP?}-\selectlanguage{english}1; Figures 3A, 3B and 4). These six
homostyles occur in two geographically close populations (M01 and T08;
Figures 2A). The occurrence of a \emph{CYP?} \selectlanguage{english}allele with no mutations in
homostyles enables us to suggest for the first time that additional
mechanisms might be responsible for the inactivation of \emph{CYP?},
\selectlanguage{english}including mutations in its promoter region or structural rearrangements
that cannot be detected \emph{via} exon sequencing. \emph{} In general,
this result underscores that mutations outside genes controlling
specific traits might represent essential sources for repeated trait
losses.

The interpretation of \emph{CYP?} \selectlanguage{english}sequencing results presented above is
compatible with multiple origins of homostyly in \emph{P. vulgaris}, but
alternative scenarios involving a single origin should be considered.
For example, as discussed above, the initial homostylous phenotype could
have arisen \emph{via} mutations in the promoter region, structural
rearrangements, or changes in the expression of \emph{CYP?} \selectlanguage{english}but not in
the \emph{CYP?} \selectlanguage{english}coding region \emph{per se}, as in the six homostyles
carrying \emph{CYP?-1,}\selectlanguage{english} the haplotype of S-morphs\emph{.} Such homostyle
could have originated once in a large ancestral population that later
became fragmented and/or migrated to neighboring populations, gaining
nonsynonymous mutations and indels in \emph{CYP?} \selectlanguage{english}after the gene had
already been rendered functionless. However, in the scenario of a single
origin of homostyly triggered by a mutation outside
\emph{CYP}\textsuperscript{T}, we would expect homostyles carrying
\emph{CYP?-1}\selectlanguage{english} to \emph{} co-occur in the same population with other
homostyles carrying mutated \emph{CYP?} \selectlanguage{english}alleles (i.e., \emph{CYP?-2 to
-9}\selectlanguage{english}), but our findings do not match this prediction (Figure 4). Thus,
the most parsimonious interpretation of the haplotype network and
phylogeny inferred from \emph{CYP?} \selectlanguage{english}sequences (Figures 3 and 4) is
congruent with multiple origins of homostyly \emph{via} independent
mutations at \emph{CYP?}. \selectlanguage{english}Furthermore, \emph{CYP?-9}\selectlanguage{english}
(\emph{S}\textsuperscript{\emph{LH2}} from Li et al., \emph{} 2016) was
found exclusively in a homostyle from Chiltern Hills located
approximately 200 km from Somerset. Given the limited dispersal
abilities reported for \emph{P. vulgaris} (Cahalan \& Gliddon, 1985),
this finding further supports independent origins of homostyly. Future
transcriptomic or RT-PCR analyses of homostyles with different
\emph{CYP}? \selectlanguage{english}alleles aimed at testing the expression of both mutated and
unmutated \emph{CYP}? \selectlanguage{english}should further elucidate the question of single
vs. multiple origins of homostyly in \emph{P. vulgaris}. Finally, our
microsatellite data do not enable us to determine whether mutations in
\emph{CYP}? \selectlanguage{english}accumulated before or after the fragmentation of a
hypothetical ancestral population. Whole genome sequences from multiple
heterostyles and homostyles aimed at improving phylogenetic resolution
and demographic analyses are necessary to elucidate the potential role
of fragmentation in the evolution of \emph{P. vulgaris} homostyles.

The multiple \emph{CYP?} \selectlanguage{english}mutations found in the homostyles of a small
geographic region in Somerset could be viewed as a puzzling result.
However, loss of gene function is usually followed by degeneration
through accumulation of mutations, or even gene loss (Sharma et al.,
2018). Additionally, lack of recombination and hemizygosity of the
S-locus should both decrease the efficacy of purifying selection on its
genes, including \emph{CYP?,} \selectlanguage{english}possibly contributing to accumulation of
mutations (Becher, Jackson, \& Charlesworth, 2020; Gossmann, Woolfit, \&
Eyre-Walker, 2011)\emph{.} Furthermore, a higher mutation rate was
reported for \emph{CYP?} \selectlanguage{english}than \emph{GLO?} \selectlanguage{english}(Huu et al., 2020). Thus, the
multiple \emph{CYP?} \selectlanguage{english}mutations of homostyles detected here are
compatible with both theoretical expectations and previous empirical
evidence. Future sequences of other S-locus genes, their paralogs, and
additional genes outside the S-locus from both \emph{P. vulgaris} and
other primroses would provide a broader context to compare selection on
\emph{CYP?} \selectlanguage{english}vs. other genes in both homostyles and heterostyles.

\textbf{4.3 Consequences of shifts to homostyly on mating system and
population genetics}

{} Changes in floral morphology can have profound effects on the mating
system (Opedal, 2018). Specifically, the reduction of anther and stigma
separation and loss of self-incompatibility of homostylous morphs should
increase selfing compared to heterostyles. Evidence of this mating
system transition in homostyles was reported for \emph{Eichhornia}
(Husband \& Barrett, 1993), \emph{Turnera} (Belaoussoff \& Shore,
1995), \emph{Amsinckia} (Schoen et al., 1997), and some \emph{Primula}
species (Yuan et al., 2017; Zhou et al., 2017; Zhong et al., 2019). In
\emph{P. vulgaris}, previous estimates of outcrossing rates in
homostyles \emph{} ranged broadly from 0.05 to 0.80 \emph{} (Bodmer,
1958, 1984; Crosby, 1958, Piper et al., 1986).

Our results from progeny arrays indicate that outcrossing rates in
homostyles are significantly lower than in heterostyles
(\emph{t}\textsubscript{\emph{m}} = 0.14 vs. 1.0, respectively), as
expected. Such outcrossing rates were estimated based exclusively on
germinated seeds. However, ungerminated seeds were likely those with
higher homozygosity, hence higher inbreeding depression, thus their
inclusion would have further decreased our estimates of outcrossing
rates. Also as expected, the frequency of homostyles was associated with
an increase in population-level selfing rates (Figure 5A). Nevertheless,
the fully homostylous population (M01) had lower selfing rates than
expected, possibly caused by increased error rates due to its small
sample size (Redelings et al., 2015). The estimated outcrossing rates of
homostyles firmly place them in the category of selfers (i.e., \emph{}
\emph{t}\textsubscript{\emph{m}} [?] 0.2; \emph{sensu} Goodwillie et al.,
2005; Schemske \& Lande, 1985). Overall, our findings clarify previous
conflicting results about the effects of homostyly on the mating system
of \emph{P. vulgaris}, confirming that transitions to homostyly increase
selfing.

{} Mating systems shape how genetic diversity is partitioned within and
among populations (Barrett, 2010; Wright et al., 2013). In concordance
with theoretical expectations, we found that intrapopulation homostyle
frequency significantly increased inbreeding coefficients (Figure 5B).
These results conform with the reduction of genetic diversity associated
with transitions from heterostyly to homostyly documented in other
systems (Husband \& Barrett, 1993; Ness et al., 2010; Yuan et al., 2017;
Zhong et al., 2019; Zhou et al., 2017). Moreover, our population genetic
results indicate that increased frequency of homostyles in populations
increases homozygosity, thus lowering intrapopulation genetic diversity
and elevating inter-population genetic differentiation. Altogether, our
findings confirm the population genetic consequences predicted for the
transition to selfing.

\textbf{4.4 Does homostyly spread among populations and why has it not
become fixed in} \emph{\textbf{Primula vulgaris?}}

{} Gene flow plays a central role in spreading advantageous mutations
(Morjan \& Rieseberg, 2004; Ralph \& Coop, 2010) and could favor the
migration of homostylous alleles among populations. In \emph{P.
vulgaris,} it has been proposed that homostyles could have migrated from
initial places of origin to neighboring populations through pollen or
seed, the former being more likely (Crosby, 1949; 1960). Previous
studies estimated that dispersal is restricted to a maximum of a few
hundred meters from the parental plant in \emph{P. vulgaris} (Cahalan \&
Gliddon, 1985). However, occasional pollen flow over 1--3 kilometers has
been reported (Van Geert, Van Rossum, \& Triest, 2010). Indeed, our
population genetic results support a model of moderate gene flow among
populations (\emph{Nm} = 1.96; Figure S2) consistent with moderate
genetic differentiation among populations
(\emph{F}\textsubscript{\emph{st}} = 0.083; Figures 6A and B), the
latter being lower than previous estimates for \emph{P. vulgaris}
(\emph{F}\textsubscript{\emph{st}} = 0.086--0.508; van Geert et al.,
2015). Furthermore, we did not detect isolation by distance (Mantel
test, \emph{P} = 0.457), further favoring the hypothesis of gene flow.
Moreover, the geographic distribution of \emph{CYP?} \selectlanguage{english}alleles in our
study indicates that 12 homostylous individuals from five populations
(D*11, T03, T04, T07 and T10) separated by 3-18 Km shared the same
\emph{CYP?} \selectlanguage{english}allele (\emph{CYP?}-\selectlanguage{english}2; Figure 4), suggesting that some
homostyles could have migrated between neighboring populations. To
summarize, Crosby (1949; 1959) proposed a single origin of homostyles
followed by their spreading \emph{via} migration. While supporting
multiple origins of homostyly, our results imply that it can also spread
\emph{via} gene flow, providing partial corroboration of Crosby's
hypothesis.

{} A crucial question remains as to why homostyly has not become fixed
in \emph{P. vulgaris} despite automatic selfing advantage (Fisher, 1941;
1949) and reproductive assurance (Piper et al., 1986). In fact, the
opposite has been found, since previous studies recorded the loss or
decrease of homostyles in revisited populations (Boyd et al., 1990;
Curtis \& Curtis, 1985). Crosby (1949) proposed that selfed progeny from
homostyles homozygous for the S-locus should have reduced fitness,
likely caused by potentially deleterious effects of dominant
homozygosity in S-locus genes. Such negative effects were suggested by
crossing experiments with S-morphs of \emph{Primula sinensis} (Mather \&
de Winton, 1941) and \emph{P. oreodoxa} (Yuan et al., 2018). However,
whether homostyles with one or two copies of the S-locus have
differences in fitness has not been experimentally tested.

{} An additional explanation for the failure of homostyly to completely
replace heterostyly is that the negative effects of inbreeding
depression (Goodwillie et al., 2005; Richards 1984) can outweigh the
advantages of homostyly (Boyd et al., 1990; Piper et al., 1984).
Theoretical models suggest that inbreeding depression values above 0.5
should prevent the invasion of selfing individuals in an outcrossing
population (Lande \& Schemske, 1985). Two recent studies investigated
the relationship between inbreeding depression and the spread of
self-compatibility, finding contrasting results. In experimental
populations of \emph{Linaria cavanillesi}, self-compatible individuals
with no inbreeding depression displaced self-incompatible individuals in
just three generations (Voillemot \& Pannell, 2017; Voillemot,
Encinas-Viso, \& Pannell, 2019). In contrast, inbreeding depression
values of 0.54 prevented the spread of self-compatible plants in
experimental patches of self-incompatible plants in \emph{Laevenworthia
alabamica} (Layman, Fernando, Herlihy, \& Busch, \emph{} 2017). Our
estimate of inbreeding depression at seed germination stage (0.58) in
\emph{P. vulgaris} is closer to the latter study and above the threshold
value of 0.5 predicted to prevent the invasion of selfing (Lande \&
Schemske, 1985). The effect of inbreeding depression in \emph{P.
vulgaris} might be even higher than in our current estimates, for it
could also act at later stages of the life cycle. Altogether, our
results imply that inbreeding depression in homostyles is sufficiently
high to prevent the fixation of homostyly within populations.

{} While theory proposes that increased selfing should purge inbreeding
depression over time by eliminating recessive deleterious alleles from
populations (Schemske \& Lande, 1985), transitions to homostyly in
\emph{P. vulgaris} may be too recent for the purging process to have
occurred. Moreover, recessive deleterious alleles could have been
re-introduced into populations \emph{via} gene flow, thus slowing the
decrease of inbreeding depression. Our study estimated moderate levels
of gene flow among populations of \emph{P. vulgaris} (\emph{Nm} = 1.96;
Supplementary Figure S2), implying that re-introduction of genetic load
can occur. In conclusion, our findings suggest that inbreeding
depression contributes to maintaining heterostylous morphs within
populations despite the reproductive advantages of homostyly.

\textbf{ACKNOWLEDGMENTS}

We thank Natural England for permits and landowners for access to
populations; Lindsey and Derek Cullen for accommodation and assistance
during fieldwork; Markus Meierhofer and Rayko Jonas for help in the
greenhouse; Francesco Rivetti for assistance in the laboratory; Members
of the Conti group for fruitful discussions. We thank Kathryn Hodgins
and three anonymous reviewers for their comments that greatly improved
this manuscript.

\textbf{AUTHOR CONTRIBUTIONS}

EMC, BK, JdV, PS and EC designed the research. EMC collected the data.
EMC, RS and BK performed laboratory work. EMC, RS, EL-B and BK analyzed
the data. EMC and EC wrote the paper and all authors provided feedback
on the manuscript.

\textbf{FUNDING}

This research was supported by the Graduate Campus office at the
University of Zurich through a GRC-Travel Grant to EMC, by the Swiss
Government Excellence Scholarship grant no. 2018.0475 to EMC, and by the
Swiss National Science Foundation grant nos. 3100-061674.00/1 and
31003A\_175556 to EC. JdV is supported in part by Swiss National Science
Foundation grant 310030\_185251.

\textbf{CONFLICT OF INTEREST STATEMENT}

The authors declare that the research was conducted in the absence of
any commercial or financial relationships that could be construed as a
potential conflict of interest.

\textbf{DATA AVAILABILITY STATEMENT}

-Microsatellite genotypes from the 613 individuals from the 22 natural
populations analyzed in this study, as well as DNA alignment
for\emph{CYP\textsuperscript{T}} haplotype network and phylogeny will be
made available in Dryad.

-Scripts for analyses and figures in R are being curated and will be
made available in Github.

\textbf{REFERENCES}

Albalat, R. \& Ca\selectlanguage{ngerman}ñestro, C. (2016). Evolution by gene loss. \emph{Nature
Reviews Genetics} 17, 379--391

Barmentlo, S. H., Meirmans, P. G., Luijten, S. H., Triest, L. \&
Oostermeijer, J. G. B. (2017). Outbreeding depression and breeding
system evolution in small, remnant populations of Primula vulgaris:
consequences for genetic rescue. \emph{Conservation Genetics}
\textbf{19}, 545--554.

Barrett, S. C. H. (1992). \emph{Evolution and function of heterostyly}.
Berlin, Germany: Springer-Verlag.

Barrett, S. C. H. (2010). Understanding plant reproductive diversity.
\emph{Philosophical Transactions of the Royal Society B: Biological
Sciences,} \emph{365}(1537), 99--109. doi:10.1098/rstb.2009.0199

Barrett, S. C. H. (2019). `A most complex marriage arrangement': Recent
advances on heterostyly and unresolved questions. \emph{New
Phytologist,} \emph{224}(3), 1051--1067. doi:10.1111/nph.16026

Barrett, S. C. H., Ness, R. W., \& Vallejo-Marín, M. (2009).
Evolutionary pathways to self-fertilization in a tristylous plant
species. \emph{New Phytologist,} \emph{183}(3), 546--556.
doi:10.1111/j.1469-8137.2009.02937.x

Beaumont, M. A., Zhang, W., \& Balding, D. J. (2002). Approximate
Bayesian Computation in population genetics. \emph{Genetics,}
\emph{162}(162), 2025--2035. doi:10.1111/j.1937-2817.2010.tb01236.x

Becher, H., Jackson, B. C. \& Charlesworth, B. (2020). Patterns of
genetic variability in genomic regions with low rates of recombination.
\emph{Curr. Biol.} 30, 94-100.e3

Belaoussoff, S., \& Shore, J. S. (1995). Floral correlates and fitness
consequences of mating-system variation in \emph{Turnera ulmifolia}.
\emph{Evolution,} \emph{49}(3), 545--556. doi:10.2307/2410278

Bodmer, W. F. (1958). Natural crossing between homostyle plants of
\emph{Primula vulgaris}. \emph{Heredity,} \emph{12}(3), 363--370.

Bodmer, W. F. (1960). The genetics of homostyly in populations of
\emph{Primula vulgaris}. \emph{Philosophical Transactions of the Royal
Society B: Biological Sciences,} \emph{242}(696), 517--549.

Bodmer, W. F. (1984). Sex and generations of primroses. \emph{Nature,}
\emph{310}, 731. doi:10.1038/310731a0

Boyd, M., Silvertown, J., \& Tucker, C. (1990). Population ecology of
heterostyle and homostyle \emph{Primula vulgaris}: Growth, survival and
reproduction in field populations. \emph{Journal of Ecology,}
\emph{78}(3), 799--813.
doi:\href{https://doi.org/10.2307/2260900\#_blank}{10.2307/2260900}

Busch, J. W., \& Delph, L. F. (2012). The relative importance of
reproductive assurance and automatic selection as hypotheses for the
evolution of self-fertilization. \emph{Annals of Botany,} \emph{109}(3),
553--562. doi:10.1093/aob/mcr219

Cahalan, C. M., \& Gliddon, C. (1985). Genetic neighbourhood sizes in
\emph{Primula vulgaris}. \emph{Heredity,} \emph{54}, 65--70.
doi:10.1038/hdy.1985.9

Carlson, M. L., Gisler, S. D., \& Kelso, S. (2008). The role of
reproductive assurance in the arctic: A comparative study of a
homostylous and distylous species pair. \emph{Arctic, Antarctic, and
Alpine Research.} \emph{40}(1), 39--47.
doi:10.1657/1523-0430(06-080){[}CARLSON{]}2.0.CO;2

Chantha, S. C., Herman, A. C., Platts, A. E., Vekemans, X., \& Schoen,
D. J. (2013). Secondary evolution of a self-incompatibility locus in the
Brassicaceae genus \emph{Leavenworthia}. \emph{PLoS Biology.}
\emph{11}(5), e1001560. doi:10.1371/journal.pbio.1001560

Charlesworth, B., \& Charlesworth, D. (1979). The maintenance and
breakdown of distyly. \emph{American Naturalist,} \emph{114}(4),
499--513. doi:\href{https://doi.org/10.1086/283497}{10.1086/283497}

Cocker, J. M., Webster, M. A., Li, J., Wright, J., Kaithakottil, G.,
Swarbreck, D., \& Gilmartin, P. M. (2015). Oakleaf: An S locus-linked
mutation of \emph{Primula vulgaris} that affects leaf and flower
development. \emph{New Phytologist,} \emph{208}(1), 149--161.
doi:10.1111/nph.13370

Cocker, J. M., Wright, J., Li, J., Swarbreck, D., Dyer, S., Caccamo, M.,
\& Gilmartin, P. M. (2018). \emph{Primula vulgaris} (primrose) genome
assembly, annotation and gene expression, with comparative genomics on
the heterostyly supergene. \emph{Scientific Reports,} \emph{8}(1),
17942. doi:10.1038/s41598-018-36304-4

Costa, J., Torices, R., \& Barrett, S. C. H. (2019). Evolutionary
history of the buildup and breakdown of the heterostylous syndrome in
Plumbaginaceae. \emph{New Phytologist}, \emph{224}(3), 1278-1289.
doi:10.1111/nph.15768

Crosby, J. L. (1940). High proportions of homostyle plants in
populations of \emph{Primula vulgaris}. \emph{Nature,} \emph{145},
672--673. doi:\href{https://doi.org/10.1038/145672c0}{10.1038/145672c0}

Crosby, J. L. (1949). Selection of an unfavourable gene-complex.
\emph{Evolution,} \emph{3}(3), 212--230.
doi:\href{https://doi.org/10.1111/j.1558-5646.1949.tb00022.x}{10.1111/j.1558-5646.1949.tb00022.x}

Crosby, J. L. (1959). Outcrossing on homostyle primroses.
\emph{Heredity,} \emph{13}, 127--131.
doi:\href{https://doi.org/10.1038/hdy.1959.9}{10.1038/hdy.1959.9}

Crosby, J. L. (1960). The use of electronic computation in the study of
random fluctuations in rapidly evolving populations. \emph{Philosophical
Transactions of the Royal Society B: Biological Sciences,}
\emph{242}(697), 551--573.
doi:\href{https://doi.org/10.1098/rstb.1960.0002}{10.1098/rstb.1960.0002}

Cutter, A. D. (2019). Reproductive transitions in plants and animals:
selfing syndrome, sexual selection and speciation. \emph{New
Phytologist} 224, 1080--1094.

Curtis, J., \& Curtis, C. F. (1985). Homostyle primroses re-visited. I.
Variation in time and space. \emph{Heredity,} \emph{54}, 227--234.
doi:10.1038/hdy.1985.30

Darwin, C. (1876). \emph{The effects of cross and self fertilization in
the vegetable kingdom}. London, UK: John Murray.

Darwin, C. (1877). \emph{The different forms of flowers on plants of the
same species.} London, UK: John Murray.

David, P., Pujol, B., Viard, F., Castella, V., \& Goudet, J. (2007).
Reliable selfing rate estimates from imperfect population genetic data.
\emph{Molecular Ecology,} \emph{16}(12), 2474--2487.
doi:10.1111/j.1365-294X.2007.03330.x

de Vos, J. M., Wüest, R. O., \& Conti, E. (2014). Small and ugly?
Phylogenetic analyses of the ``selfing syndrome'' reveal complex
evolutionary fates of monomorphic primrose flowers. \emph{Evolution,}
\emph{68}(4), 1042--1057. doi:10.1111/evo.12331

de Vos, J. M., Hughes, C. E., Schneeweiss, G. M., Moore, B. R., \&
Conti, E. (2014). Heterostyly accelerates diversification \emph{via}
reduced extinction in Darwin's primroses. \emph{Philosophical
Transactions of the Royal Society B: Biological Sciences,}
\emph{281}(1784), 20140075. https://doi.org/10.1098/rspb.2014.0075

de Vos, J. M., Keller, B., Isham, S. T., Kelso, S., \& Conti, E. (2012).
Reproductive implications of herkogamy in homostylous primroses:
Variation during anthesis and reproductive assurance in alpine
environments. \emph{Functional Ecology,} \emph{26}(4), 854--865.
doi:10.1111/j.1365-2435.2012.02016.x

de Vos, J. M., Keller, B., Zhang, L. R., Nowak, M. D., \& Conti, E.
(2018). Mixed mating in homostylous species: Genetic and experimental
evidence from an alpine plant with variable herkogamy, \emph{Primula
halleri}. \emph{International Journal of Plant Sciences,} \emph{179}(2),
87--99. doi:10.1086/695527

Dieringer, D., \& Schlötterer, C. (2003). Microsatellite Analyser (MSA):
A platform independent analysis tool for large microsatellite data sets.
\emph{Molecular Ecology Notes,} \emph{3}(1), 167-169.
doi:\href{https://doi.org/10.1046/j.1471-8286.2003.00351.x}{10.1046/j.1471-8286.2003.00351.x}

Doyle, J., \& Doyle, J. (1987). A rapid DNA isolation procedure for
small quantities of fresh leaf tissue. \emph{Phytochemical Bulletin},
\emph{19}, 11--15.

\protect\hypertarget{__DdeLink__3968_383804172}{}{} Dulberger, R.
(1992). Floral polymorphism and their functional significance in the
heterostylous syndrome. In S. C. H. Barrett (Ed.), \emph{Evolution and
Function of Heterostyly} (pp. 41--84). Berlin, Germany: Springer-Verlag.

Dwyer, K. G., Berger, M. T., Ahmed, R., Hritzo, M. K., Mcculloch, A. A.,
Price, M. J., \ldots{} Nasrallah, M. E. (2013). Molecular
characterization and evolution of self-incompatibility genes in
\emph{Arabidopsis thaliana}: The case of the \emph{Sc} haplotype.
\emph{Genetics,} \emph{193}(3), 985--994.
doi:10.1534/genetics.112.146787

Evanno, G., Regnaut, S., \& Goudet, J. (2005). Detecting the number of
clusters of individuals using the software STRUCTURE: A simulation
study. \emph{Molecular Ecology,} \emph{14}(8), 2611--2620.
doi:10.1111/j.1365-294X.2005.02553.x

Felsenstein, J. (2005). PHYLIP (Phylogeny Inference Package) version
3.6. Distributed by the author. Department of Genome Sciences,
University of Washington, Seattle.

Fisher, R. A. (1949). A theoretical system of selection for homostyle
\emph{Primula}. \selectlanguage{polish}\emph{Sankhyā:} \selectlanguage{english}\emph{The Indian Journal of Statistics,}
\emph{9}(4), 325--342. doi:10.1002/0471667196.ess1229.pub2

Foxe, J. P., Stift, M., Tedder, A., Haudry, A., Wright, S. .I, \& Mable,
B. K. (2010). Reconstructing origins of loss of self-incompatibility and
selfing in North American Arabidopsis lyrata: A population genetic
context. \emph{Evolution,} 64, 3495--3510

Ganders, F. R. (1979). The biology of heterostyly. \emph{New Zealand
Journal of Botany,} \emph{17}(4), 607--635.
doi:10.1080/0028825X.1979.10432574

Gao, H., Williamson, S., \& Bustamante, C. D. (2007). A Markov Chain
Monte Carlo approach for joint inference of population structure and
inbreeding rates from multilocus genotype data. \emph{Genetics,}
\emph{176}(3), 1635--1651. doi:10.1534/genetics.107.072371

Gilbert, K. J., Andrew, R. L., Bock, D. G., Franklin, M. T., Kane, N.
C., Moore, J. -S., \ldots{} Vines T. H. (2012). Recommendations for
utilizing and reporting population genetic analyses: The reproducibility
of genetic clustering using the program STRUCTURE. \emph{Molecular
Ecology,} \emph{21}(20), 4925--4930.
doi:10.1111/j.1365-294X.2012.05754.x

Goodwillie, C. (1999). Multiple origins of self-compatibility in
Linanthus section Leptosiphon (Polemoniaceae): Phylogenetic evidence
from internal-transcribed-spacer sequence data. \emph{Evolution,} 53,
1387--1395

Goodwillie, C., Kalisz, S., \& Eckert, C. G. (2005). The evolutionary
enigma of mixed mating systems in plants: Occurrence, theoretical
explanations, and empirical evidence. \emph{Annual Review of Ecology,
Evolution and Systematics,} \emph{36}(1), 47--79.
doi:10.1146/annurev.ecolsys.36.091704.175539

Gossmann, T. I., Woolfit, M. \& Eyre-Walker, A. (2011). Quantifying the
variation in the effective population size within a genome.
\emph{Genetics} 189, 1389--1402

Guggisberg, A., Mansion, G., Kelso, S. \& Conti, E. (2006). Evolution of
biogeographic patterns, ploidy levels, and breeding systems in a
diploid-polyploid species complex of Primula. \emph{New Phytologist}
171, 617--632

Hamilton N.E., \& Ferry M. (2018). ggtern: Ternary Diagrams using
ggplot2. \emph{Journal of Statistical Software}, \emph{87}(3), 1--17.
doi:10.18637/jss.v087.c03

Hamrick, J. L., \& Godt, M. J. W. (1996). Effects of life history traits
on genetic diversity in plant species. \emph{Philosophical Transactions
of the Royal Society B: Biological Sciences,} \emph{351}(1345),
1291--1298. doi:10.1098/rstb.1996.0112

Hey, J., \& Nielsen, R. (2004). Multilocus methods for estimating
population sizes, migration rates and divergence time, with applications
to the divergence of \emph{Drosophila pseudoobscura} and \emph{D.
persimilis}. \emph{Genetics,} \emph{167}(2), 747--760.
doi:10.1534/genetics.103.024182

Husband, B. C., \& Barrett, S. C. H. (1993). Multiple origins of
self-fertilization in tristylous \emph{Eichhornia paniculata}
(Pontederiaceae): Inferences from style morph and isozyme variation.
\emph{Journal of Evolutionary Biology,} \emph{6}(4), 591--608.
doi:\href{https://doi.org/10.1046/j.1420-9101.1993.6040591.x}{10.1046/j.1420-9101.1993.6040591.x}

Huu, C. N., Kappel, C., Keller, B., Sicard, A., Takebayashi, Y.,
Breuninger, H., \ldots{} Lenhard, M. (2016). Presence versus absence of
\emph{CYP734A50} underlies the style-length dimorphism in primroses.
\emph{Elife,} \emph{5}, 1--15. doi:10.7554/eLife.17956

Huu, C. N., Keller, B., Conti, E., Kappel, C., \& Lenhard, M. (2020).
Supergene evolution via stepwise duplications and neofunctionalization
of a floral-organ identity gene. \emph{Proceedings of the National
Academy of Science,} \emph{117}(37), 23148--23157.
doi:10.1073/pnas.2006296117

Jacquemyn, H., Endels, P., Brys, R., Hermy, M., \& Woodell, S. R. J.
(2009). Biological flora of the British Isles: \emph{Primula vulgaris}
Huds. (\emph{P. acaulis} (L.) Hill). \emph{Journal of Ecology,}
\emph{97}(4), 812--833. doi:10.1111/j.1365-2745.2009.01513.x

Johnson, M. G., Gardner, E. M., Liu, Y., Medina, R., Goffinet, B., Shaw,
A. J., \ldots{} Wickett, N. J. (2016). HybPiper: Extracting coding
sequence and introns for phylogenetics from high-throughput sequencing
reads using target enrichment. \emph{Applications in Plant Sciences,}
\emph{4}(7), 1600016. doi:10.3732/apps.1600016

Jombart, T. (2008). Adegenet: A R package for the multivariate analysis
of genetic markers. \emph{Bioinformatics,} \emph{24}, 1403--1405.
doi:10.1093/bioinformatics/btn129

Kappel, C., Huu, C. N., \& Lenhard, M. (2017). A short story gets
longer: Recent insights into the molecular basis of heterostyly.
\emph{Journal of Experimental Botany,} \emph{68}(21-22), 5719--5730.
doi:10.1093/jxb/erx387

Keenan, K., Mcginnity, P., Cross, T. F., Crozier, W. W., \& Prod\selectlanguage{ngerman}öhl, P.
A. (2013). DiveRsity: An R package for the estimation and exploration of
population genetics parameters and their associated errors.
\emph{Methods in Ecology and Evolution,} \emph{4}(8), 782--788.
doi:10.1111/2041-210X.12067

Keller, B., Ganz, R., Mora-Carrera, E., Nowak, M. D., Theodoridis, S.,
Koutroumpa, K., \ldots{} Conti, E. (2020). Asymmetries of reproductive
isolation are reflected in directionalities of hybridization:
Integrative evidence on the complexity of species boundaries. \emph{New
Phytologist.} doi:10.1111/nph.16849

Keller, B., Thomson, J. D., \& Conti, E. (2014). Heterostyly promotes
disassortative pollination and reduces sexual interference in Darwin's
primroses: Evidence from experimental studies. \emph{Functional
Ecology,} \emph{28}(6), 1413--1425. doi:10.1111/1365-2435.12274

Kissling, J., \& Barrett, S. C. H. (2013). Variation and evolution of
herkogamy in \emph{Exochaenium} (Gentianaceae): Implications for the
evolution of distyly. \emph{Annals of Botany,} \emph{112}(1), 95--102.
doi:10.1093/aob/mct097

Kohn, J. R., Graham, S. W., Morton, B., Doyle, J. J., \& Barrett, S. C.
H. (1996). Reconstruction of the evolution of reproductive characters in
Pontederiaceae using phylogenetic evidence from chloroplast DNA
restriction-site variation. \emph{Evolution,} \emph{50}(4), 1454-1469.
doi:\href{https://doi.org/10.1111/j.1558-5646.1996.tb03919.x}{1558-5646.1996.tb03919.x}

Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., \&
Mayrose, I. (2015). Clumpak: A program for identifying clustering modes
and packaging population structure inferences across K. \emph{Molecular
Ecology Resources,} \emph{15}(5), 1179--1191.
doi:10.1111/1755-0998.12387

Kozlov, A. M., Darriba, D., Flouri, T., Morel, B., \& Stamatakis, A.
(2019). RAxML-NG: A fast, scalable and user-friendly tool for maximum
likelihood phylogenetic inference. \emph{Bioinformatics,} \emph{35}(21),
4453--4455. doi:10.1093/bioinformatics/btz305

Kumar, S., Stecher, G., Li, M., Knyaz, C., \& Tamura, K. (2018). MEGA X:
Molecular evolutionary genetics analysis across computing platforms.
\emph{Molecular Biology and Evolution,} \emph{35}(6), 1547--1549.
doi:10.1093/molbev/msy096

Lande, R., \& Schemske, D. W. (1985). The evolution of
self-fertilization and inbreeding depression in plants. I. Genetic
Models. \emph{Evolution,} \emph{39}(1), 24--40. doi:10.2307/2408514

Layman, N. C., Fernando, M. T. R., Herlihy, C. R., \& Busch, J. W.
(2017). Costs of selfing prevent the spread of a self-compatibility
mutation that causes reproductive assurance. \emph{Evolution,}
\emph{71}(4), 884--897. doi:10.1111/evo.13167

Lewis, D., \& Jones, D. A. (1992). The genetics of heterostyly. In S. C.
H. Barrett (Ed.), \emph{Evolution and Function of Heterostyly} (pp.
129--151). Berlin, Germany: Springer-Verlag.

Li, J., Cocker, J. M., Wright, J., Webster, M. A., McMullan, M., Dyer,
S., \ldots{} Gilmartin, P. M. (2016). Genetic architecture and evolution
of the S locus supergene in \emph{Primula vulgaris}. \emph{Nature
Plants,} \emph{2}(12), 16188. doi:10.1038/nplants.2016.188

Li, J., Webster, M. A., Wright, J., Cocker, J. M., Smith, M. C.,
Badakshi, F., \ldots{} Gilmartin, P. M. (2015). Integration of genetic
and physical maps of the \emph{Primula vulgaris} S locus and
localization by chromosome in situ hybridization. \emph{New
Phytologist,} \emph{208}(1), 137--148. doi:10.1111/nph.13373

Lloyd, D. G., \& Webb, C. J. (1992). The selection of heterostyly. In S.
C. H. Barrett (Ed.), \emph{Evolution and Function of Heterostyly} (pp.
179--207). Berlin, Germany: Springer-Verlag.

Lopes, J. S., Balding, D., \& Beaumont, M. A. (2009). popABC: A program
to infer historical demographic parameters. \emph{Bioinformatics,}
\emph{25}(20), 2747--2749. doi:10.1093/bioinformatics/btp487

Maddison, W. P. and D.R. Maddison. (2019). Mesquite: a modular system
for evolutionary analysis. Version 3.61. Retrieved from
http://www.mesquiteproject.org

Mast, A. R., Kelso, S., \& Conti, E. (2006). Are any primroses
(\emph{Primula}) primitively monomorphic? \emph{New Phytologist,}
\emph{171}(3), 605--616.
doi:\href{https://doi.org/10.1111/j.1469-8137.2006.01700.x}{10.1111/j.1469-8137.2006.01700.x}

Mather, K. \& Winton, D. DE. (1941). Adaptation and counter-adaptation
of the breeding system in Primula. \emph{Annals of Botany} 5, 297--311.

Meirmans, P. G., \& Hedrick, P. W. (2011). Assessing population
structure: FST and related measures. \emph{Molecular Ecology Resources,}
\emph{11}(1), 5--18. doi:10.1111/j.1755-0998.2010.02927.x

Morjan, C. L., and Rieseberg, L. H. (2004). How species evolve
collectively: Implications of gene flow and selection for the spread of
advantageous alleles. \emph{Molecular Ecology,} \emph{13}(6),
1341--1356. doi:10.1111/j.1365-294X.2004.02164.x

Ness, R. W., Wright, S. I., \& Barrett, S. C. H. (2010). Mating-system
variation, demographic history and patterns of nucleotide diversity in
the tristylous plant \emph{Eichhornia paniculata.} \emph{Genetics,}
184(2), 381--392. doi:10.1534/genetics.109.110130

Nielsen, R., \& Wakely, J. (2001). Distinguishing migration from
isolation: A Markov Chain Monte Carlo approach. \emph{Genetics,}
\emph{158}(2), 885--896.

Nowak, M. D., Russo, G., Schlapbach, R., Huu, C. N., Lenhard, M., \&
Conti, E. (2015). The draft genome of \emph{Primula veris} yields
insights into the molecular basis of heterostyly. \emph{Genome Biology,}
\emph{16}(1), 12. doi:10.1186/s13059-014-0567-z

Opedal, Ø. H. (2018). Herkogamy, a principal functional trait of plant
reproductive biology. \emph{International Journal of Plant Sciences,}
\emph{179}(9), 677--687. doi:10.1086/700314

Paradis, E. (2010). Pegas: An R package for population genetics with an
integrated-modular approach. \emph{Bioinformatics,} \emph{26}(3),
419--420. doi:10.1093/bioinformatics/btp696

Pérez-Barrales, R. \& Arroyo, J. (2010). Pollinator shifts and the loss
of style polymorphism in Narcissus papyraceus (Amaryllidaceae).
\emph{Journal of Evolutionary Biology} 23, 1117--1128

Pérez-Barrales, R., Pino, R., Albaladejo, R. G. \& Arroyo, J. (2009).
Geographic variation of flower traits in Narcissus papyraceus
(Amaryllidaceae): Do pollinators matter? \emph{Journal of Biogeography}
36, 1411--1422

Piper, J. G., \& Charlesworth, B. (1986). The evolution of distyly in
\emph{Primula vulgaris}. \emph{Biological Journal of the Linnean
Society,} \emph{29}(2), 123--137. doi:10.1111/j.1095-8312.1986.tb01827.x

Piper, J. G., Charlesworth, B., \& Charlesworth, D. (1984). A high rate
of self-fertilization and increased seed fertility on homostyle
primroses. \emph{Nature,} \emph{310}(5972), 50-51. doi:10.1038/309126a0

Piper, J. G., Charlesworth, B., \& Charlesworth, D. (1986). Breeding
system evolution in \emph{Primula vulgaris} and the role of reproductive
assurance. \emph{Heredity,} \emph{56}(2), 207--217.
doi:10.1038/hdy.1986.33

Pritchard, J. K., Stephens, M., \& Donnelly, P. (2000). Inference of
population structure using multilocus genotype data. \emph{Genetics,}
\emph{155}(2), 945--959. doi:10.1111/j.1471-8286.2007.01758.x

Ralph, P., \& Coop, G. (2010). Parallel adaptation: One or many waves of
advance of an advantageous allele? \emph{Genetics,} \emph{186}(2),
647--668. doi:10.1534/genetics.110.119594

Redelings, B. D., Kumagai, S., Tatarenkov, A., Wang, L., Sakai, A. K.,
Weller, S. G., \ldots{} Uyenoyama, M. K. (2015). A Bayesian approach to
inferring rates of selfing and locus-specific mutation. \emph{Genetics,}
\emph{201}(3), 1171--1188. doi:10.1534/genetics.115.179093

Richards, A. J. (2003). \emph{Primula} (2\textsuperscript{nd} ed.).
Portland, USA: Timber Press.

Rigby, R. A., Stasinopoulos, D. M., \& Lane, P. W. (2005). Generalized
additive models for location, scale and shape. \emph{Applied
Statistics,} \emph{54}(3), 507--554.
doi:10.1111/j.1467-9876.2005.00510.x

Ritland, K. (2002). Extensions of models for the estimation of mating.
\emph{Heredity,} \emph{88}(4), 221--228. doi:10.1038/sj/hdy/6800029

Rousset, F. (2008). GENEPOP'007: A complete re-implementation of the
GENEPOP software for Windows and Linux. \emph{Molecular Ecology
Resources,} \emph{8}(1), 103--106. doi:10.1111/j.1471-8286.2007.01931.x

Ruiz-Martín, J., Santos-Gally, R., Escudero, M., Midgley, J. J.,
Pérez-Barrales, R., \& Arroyo, J. (2018). Style polymorphism in Linum
(Linaceae): a case of Mediterranean parallel evolution? \emph{Plant
Biology} 20, 100--111

Schemske, D. W., \& Lande, R. (1985). The evolution of
self-fertilization and inbreeding depression in plants. II. Empirical
observations. \emph{Evolution,} \emph{39}(1), 41--52.
doi:10.1111/j.1558-5646.1985.tb04078.x

Schoen, D. J., Johnston, M. O., L'Heureux, A.-M., \& Marsolais J. V.
(1997). Evolutionary history of the mating system in \emph{Amsinckia}
(Boraginaceae). \emph{Evolution,} \emph{51}(4), 1090--1099.
doi:\href{https://doi.org/10.1111/j.1558-5646.1997.tb03956.x}{10.1111/j.1558-5646.1997.tb03956.x}

Shao, J. W., Wang, H. F., Fang, S. P., Conti, E., Chen, Y. J., \& Zhu,
H. M. (2019). Intraspecific variation of self-incompatibility in the
distylous plant \emph{Primula merrilliana}. \emph{AoB Plants,}
\emph{11}(3). doi:10.1093/aobpla/plz030

Sharma, V., Hecker, N., Roscito, J. G., Foerster, L., Bjoern, E. L., \&
Hiller, M. (2018). A genomics approach reveals insights into the
importance of gene losses for mammalian adaptations. \emph{Nature
Communications} 9, 1--9

Shimizu, K. K., Shimizu-Inatsugi, R., Tsuchimatsu, T., \& Purugganan, M.
D. (2008). Independent origins of self-compatibility in
\emph{Arabidopsis thaliana}. \emph{Molecular Ecology,} \emph{17}(2),
704--714. doi:10.1111/j.1365-294X.2007.03605.x

Shimizu, K. K., \& Tsuchimatsu, T. (2015). Evolution of selfing:
Recurrent patterns in molecular adaptation. \emph{Annual Review of
Ecology, Evolution and Systematics,} \emph{46}(1), 593--622.
doi:10.1146/annurev-ecolsys-112414-054249

Signorell, A., Aho, K., Alfons, A., Anderegg, N., Aragon, A.,
Arachchige, A., \ldots{} Zeileis (2020). DescTools: Tools for
Descriptive Statistics. R package version 0.99.39. Retrieved from
https://cran.r-project.org/package=DescTools.

Simmons, M. P., \& Ochoterena, H. (2000). Gaps as characters in
sequence-based phylogenetic analyses. \emph{Systematic Biology,}
\emph{49}(2), 369--381.
doi:\href{https://doi.org/10.1093/sysbio/49.2.369}{10.1093/sysbio/49.2.369}

Triest, L., Sierens, T., \& Van Rossum, F. (2015). Multiplexing 15
microsatellite loci for European primrose (\emph{Primula vulgari}s).
\emph{Conservation Genetic Resources,} \emph{7}(1), 279--281.
doi:10.1007/s12686-014-0357-7

Valverde, T., \& Silvertown, J. (1995). Spatial variation in the seed
ecology of a woodland herb (\emph{Primula vulgaris}) in relation to
light environment. \emph{Functional Ecology,} \emph{9}(6), 942--950.
doi:10.2307/2389993

Van Geert, A., Van Rossum, F., \& Triest, L. (2010). Do linear landscape
elements in farmland act as biological corridors for pollen dispersal?
\emph{Journal of Ecology,} \emph{98}(1), 178--187.
doi:10.1111/j.1365-2745.2009.01600.x

Van Oosterhout, C., Hutchinson, W. F., Wills, D. P. M., \& Shipley, P.
(2004). MICRO-CHECKER: Software for identifying and correcting
genotyping errors in microsatellite data. \emph{Molecular Ecology
Notes,} \emph{4}(3), 535--538. doi:10.1111/j.1471-8286.2004.00684.x

Vekemans, X., Poux, C., Goubet, P. M. \& Castric, V. (2014). The
evolution of selfing from outcrossing ancestors in Brassicaceae: What
have we learned from variation at the S-locus? \emph{Journal of
Evolutionary Biology,} 27, 1372--1385

Voillemot, M., Encinas-Viso, F., \& Pannell, J. R. (2019). Rapid loss of
self-incompatibility in experimental populations of the perennial
outcrossing plant \emph{Linaria cavanillesii}. \emph{Evolution,}
\emph{73}(5), 913--926. doi:10.1111/evo.13721

Voillemot, M., and Pannell, J. R. (2017). Inbreeding depression is high
in a self-incompatible perennial herb population but absent in a
self-compatible population showing mixed mating. \emph{Ecology and
Evolution,} \emph{7}(20), 8535--8544. doi:10.1002/ece3.3354

Weir, B. S., \& Cockerham, C. C. (1984). Estimating \emph{F}-statistics
for the analysis of population structure. \emph{Evolution,}
\emph{38}(6), 1358--1370.
doi:\href{https://doi.org/10.1111/j.1558-5646.1984.tb05657.x}{10.1111/j.1558-5646.1984.tb05657.x}

Woodell, S. R. J. (1960). What pollinates primulas? \emph{New
Scientist}, \emph{8}, 568--571.

Wright, S. I., Kalisz, S., \& Slotte, T. (2013). Evolutionary
consequences of self-fertilization in plants. \emph{Philosophical
Transactions of the Royal Society B: Biological Sciences,}
\emph{280}(1760), 20130133. doi:10.1098/rspb.2013.0133

Yuan, S., Barrett, S. C. H., Duan, T., Qian, X., Shi, M., and Zhang, D.
(2017). Ecological correlates and genetic consequences of evolutionary
transitions from distyly to homostyly. \emph{Annals of Botany,}
\emph{120}(5), 775--789. doi:10.1093/aob/mcx09

Yuan, S., Barrett, S.C.H., Li, C., Li, X., Xie, K., and Zhang D. (2018).
Genetics of distyly and homostyly in a self-compatible Primula.
\emph{Heredity,} 1--10 doi:10.1038/s41437-018-0081-2

Zhang, J. (2000). Rates of conservative and radical nonsynonymous
mucleotide substitutions in Mammalian nuclear genes. \emph{Journal of
Molecular Evolution}, \emph{50}(1), 56-68. doi:10.1007/s002399910007

Zhao, Z., Luo, Z., Yuan, S., Mei, L., \& Zhang, D. (2019). Global
transcriptome and gene co-expression network analyses on the development
of distyly in \emph{Primula oreodoxa}. \emph{Heredity,} \emph{123}(6),
784--794. doi:10.1038/s41437-019-0250-y

Zhong, L., Barrett, S. C. H., Wang, X., Wu, Z., Sun, H., Li, D., \&
Zhou, W. (2019). Phylogenomic analysis reveals multiple evolutionary
origins of selfing from outcrossing in a lineage of heterostylous
plants. \emph{New Phytologist,} \emph{224}(3), 1290--1303.
doi:10.1111/nph.15905

\protect\hypertarget{__UnoMark__72215_2042567327}{}{} Zhou, W., Barrett,
S. C. H., Li, H. D., Wu, Z. K., Wang, X. J., Wang, H., \& Li, D. Z.
(2017). Phylogeographic insights on the evolutionary breakdown of
heterostyly. \emph{New Phytologist,} \emph{214}(3), 1368--1380.
doi:10.1111/nph.14453\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-1/Fig-1}
\caption{{\textbf{Figure 1.-} Heterostyly and homostyly in \emph{Primula
vulgaris}. Heterostylous floral morphs can be either short-styled
(S-morphs) or long-styled (L-morphs). They are characterized by
reciprocal placement of male (anthers) and female (stigmas) sexual
organs (i.e., reciprocal herkogamy), promoting outcrossing.
Additionally, a physiological self-incompatibility mechanism prevents
self- and intramorph fertilization. Homostylous floral morphs (H-morphs)
are characterized by reduced or no herkogamy and are self-compatible,
favoring selfing. Most homostyles have both stigma and anthers at the
mouth of the corolla tube (i.e., long-homostyly); short homostyly is
rare (not shown). Long-homostyles stem from S-morphs and retain their
compatibility type, thus can fertilize only L-morphs. Conversely,
S-morphs cannot fertilize long homostyles due to stigma clogging in the
latter. Red arrows indicate compatible, likely pollinations between
floral morphs. Anthers are represented in yellow and stigmas in green.
{\label{277096}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-2-D/Fig-2-D}
\caption{{\textbf{Figure 2.-}~\textbf{A)} Geographic distribution of 22
populations of~\emph{Primula vulgaris} from Somerset, England. Pie
charts represent intra-population frequencies of long-styled (L; white),
short-styled (S; black) and homostylous (H; grey) floral morphs;
dimorphic (D) populations consist of L- and S-individuals, except for
population D\textsuperscript{*}11, which consists of L- and
H-individuals; trimorphic (T) populations consist of L-, S- and
H-individuals; the single monomorphic (M) population consists of
H-individuals.~\textbf{B)} Ternary plot representing the frequencies of
S-, L- and H-individuals sampled from the same 22 populations. Each
vertex of the triangle represents complete population monomorphy for
each floral morph; each side of the triangle represents different levels
of population dimorphism (bottom: H-morph absent, only L- and S-morphs;
left side: S-morph absent, only L- and H-morphs; right side: L-morph
absent, only S- and H-morphs); each point inside the triangle represents
trimorphic populations; arrow pointing to the base of the triangle
represents equal frequencies of S- and L- morphs (i.e., isoplethy); the
black dashed line represents the trajectory from isoplethy in distylous
populations to complete homostyly observed in the sampled populations;
the leftward shift of the line is caused by the greater reduction of S-
than L-individuals; the grey dashed line represents the trajectory from
isoplethy in distylous populations to complete homostyly with equal
reduction of S- vs. L-individuals.
{\label{599874}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-3/Fig-3}
\caption{{\textbf{Figure 3.-}~~Variation of \emph{CYP?} detected in 44 individuals
of \emph{Primula vulgaris} (17 short-styled and 27 homostylous
individuals) from ten natural populations sampled for this study\emph{,}
plus two \emph{CYP?} alleles (indicated with *) from two homostylous
individuals previously reported as \emph{CYP?
S}\textsuperscript{\emph{LH1}} and \emph{S}\textsuperscript{\emph{LH2}}
by Li et al. (2016); \emph{CYP?} is a 1587 bp long hemizygous gene that
comprises five exons (base-pair lengths of each exon are indicated in
parentheses). \textbf{A)} Graphical representation of the nine
\emph{CYP?} alleles with different types of mutations in
\emph{CYP?2--9}: vertical lines, white triangles, and black triangles
represent point mutations, deletions, and insertions, respectively, with
positions of each mutation reported above each exon; one synonymous, and
five nonsynonymous mutations, two deletions, and one insertion were
found. \textbf{B)} Haplotype network of the nine \emph{CYP?} alleles:
each circle (i.e., node in the network) represents a different
\emph{CYP?} allele (\emph{CYP?1--9}, listed next to each circle), with
circle size proportional to the number of homostylous (grey) and
short-styled (black) individuals carrying that allele (reported inside
the circle); nodes are connected by solid lines representing the lowest
number of mutational steps between alleles, in this case always a single
mutational step indicated by a single dash across the line, except for
\emph{CYP?-1} and \emph{CYP?-5} connected by two mutational steps; grey
dashed lines represent alternative connections among haplotypes with two
mutational steps between alleles, indicated by a double dash across the
line.
{\label{488201}}%
}\selectlanguage{english}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Screenshot-from-2021-04-28-14-33-57/Screenshot-from-2021-04-28-14-33-57}
\caption{{\textbf{Figure 4.-} Maximum likelihood tree of~\emph{CYP?}sequences from
44 individuals of~\emph{Primula vulgaris} (17 short-styled and 27
homostylous, boldfaced) sampled for this study, plus~\emph{CYP?}
sequences from two homostylous individuals (indicated with *) previously
reported as~\emph{CYP? S}\textsuperscript{\emph{LH1}}
and~\emph{S}\textsuperscript{\emph{LH2}} (Li et al., 2016), and
a~\emph{CYP?~}sequence from an S-individual of~\emph{Primula veris}
(GenBank: KX589238) as outgroup. Each accession at the tips is labeled
from left to right with population number, individual number, and floral
morph type (S and H). Accessions from H-individuals are boldfaced;
accessions from the same population have the same color (see also Fig.
2). Bootstrap support values [?] 80 are indicated below to the branches.
{\label{288412}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-5/Fig-5}
\caption{\selectlanguage{ngerman}{\textbf{Figure 5.-} Correlations between population genetic parameters
(y axis) and frequency of homostyles (x axis) in 22 natural populations
of \emph{Primula vulgaris} inferred from 12 microsatellites: A)
Population-level estimates ± 95 \% CI of selfing rates based on
multilocus linkage disequilibrium (See Methods) and B) Inbreeding
coefficient ± 95 \% CI.
{\label{692091}}%
}}
\end{center}
\end{figure}\selectlanguage{ngerman}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig-7/Fig-7}
\caption{{\textbf{Figure 6.-} Patterns of genetic differentiation for the 22
natural populations of~\emph{Primula vulgaris} . A) Results from
InStruct analysis of microsatellites assigning all populations to two
genetic clusters (k=2): Optimal number of clusters was selected
according to Evanno et al. (2009); B) Plot of the two first PCs from
Discriminant Analysis of Principal Components (DAPC) separating
population T08 (blue circles) from the rest of the populations
{\label{692135}}%
}}
\end{center}
\end{figure}

\par\null

\textbf{Table 1.} Summary data for 22 natural populations of Primula
vulgaris sampled in Somerset, England, with floral morph frequencies,
levels of inbreeding (Fis), and estimates of population-level selfing
rates (s(BES)) inferred using 12 microsatellites. Dimorphic (D)
populations consist of long-styled (L) and short-styled (S) morphs,
except for population D\textsuperscript{*}11, which consists of LS- and
homostylous (HO) morphs; trimorphic (T) populations consist of LS-, SS-
and HO-morphs; the single monomorphic (M) population consists of
HO-morphs.

\par\null\par\null\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Table-1/Table-1}

\end{center}
\end{figure}

Fis varies from -1 (heterozygotes excess) to + 1 (homozygotes excess due
to inbreeding): values closer to 0 indicate random outcrossing under
HWE; population level estimates of selfing rates vary from 0 (maximum
outcrossing) to 1 (maximum self-fertilization).

\selectlanguage{english}
\FloatBarrier
\end{document}

