How rapidly do self‐compatible populations evolve selfing? Mating system estimation within recently evolved self‐compatible populations of Azorean Tolpis succulenta (Asteraceae)

Abstract Genome‐wide genotyping and Bayesian inference method (BORICE) were employed to estimate outcrossing rates and paternity in two small plant populations of Tolpis succulenta (Asteraceae) on Graciosa island in the Azores. These two known extant populations of T. succulenta on Graciosa have recently evolved self‐compatibility. Despite the expectation that selfing would occur at an appreciable rate (self‐incompatible populations of the same species show low but nonzero selfing), high outcrossing was found in progeny arrays from maternal plants in both populations. This is inconsistent with an immediate transition to high selfing following the breakdown of a genetic incompatibility system. This finding is surprising given the small population sizes and the recent colonization of an island from self‐incompatible colonists of T. succulenta from another island in the Azores, and a potential paucity of pollinators, all factors selecting for selfing through reproductive assurance. The self‐compatible lineage(s) likely have high inbreeding depression (ID) that effectively halts the evolution of increased selfing, but this remains to be determined. Like their progeny, all maternal plants in both populations are fully outbred, which is consistent with but not proof of high ID. High multiple paternity was found in both populations, which may be due in part to the abundant pollinators observed during the flowering season.


| INTRODUC TI ON
The evolution of self-fertilization (selfing) is one of the most common trends in flowering plants (Barrett, 2013;Grossenbacher et al., 2017;Stebbins, 1957;Wright et al., 2013). The first step in this transition is a change in the breeding system with the loss of genetic self-incompatibility (SI). This loss, often considered a unidirectional transition (Barrett, 2013;Goodwillie, 1999;Herman & Schoen, 2016;Layman et al., 2017;Wright et al., 2013), has pronounced consequences for the subsequent evolution of a species. Most obviously, the loss of SI enables a change from a highly outcrossing mating system (who mates with whom and how frequently in populations) to the potential for higher rates of self-fertilization, with cascading effects on the distribution of genetic variation among individuals and populations, and on the balance between mutation and natural selection (Kelly, 1999;Lande & Schemske, 1985). Microevolutionary changes have macroevolutionary consequences (Cheptou, 2018;Igić & Busch, 2013).
Comparative studies in at least four plant families have shown that self-compatible (SC) lineages have higher extinction rates and lower net-diversification rates than closely related SI lineages (Freyman & Höhna, 2018;Gamisch et al., 2015;Goldberg et al., 2010;de Vos et al., 2014).
Comparative studies do not indicate how rapidly the loss of SI leads to mating system change. Natural populations routinely harbor genetic variation in reproductive traits that determine selfing rate and thus have the capacity to rapidly evolve selfing (e.g., Bodbyl-Roels & Kelly, 2011;Thomann et al., 2013). However, ecological and selective factors may prevent a rapid transition.
Selfing is favored both by reproductive assurance and automatic selection. The latter refers to the 3:2 transmission advantage enjoyed by selfers when they fertilize their own ovules as well as outcrossing (Fisher, 1941;Jain, 1976;. While these are strong forces, high outcrossing can be maintained if selfed progeny are much less fit than outcrossed progeny (inbreeding depression, hereafter indicated as ID, Arista et al., 2017;Husband & Schemske, 1996), if mates are abundant, or if factors like pollen or seed discounting undermine automatic selection. In this study, we examine the tempo of selfing transitions by estimating the mating system in recently evolved SC populations of Tolpis succulenta (Asteraceae: Cichorieae). This is a species with ecological and demographic characteristics that should rapidly select for high selfing (Figure 1).
The genus Tolpis (Asteraceae: Cichorieae) is a small monophyletic group occurring mainly in the Macaronesian archipelagos, especially the Canary Islands (Jarvis, 1980;Mort et al., 2015). The breeding system of Tolpis is overwhelmingly SI with little or no self-seed produced in most populations (Crawford et al., 2008. However, SI systems can be "leaky" (pseudo-self-compatible or PSC; Levin, 1996) and greenhouse experiments on Tolpis plants from some populations can produce low levels of selfseed when denied cross-pollen (Crawford et al., 2008. In this paper, we focus on T. succulenta sensu lato, a species that was described from Madeira, but later considered to also occur in the Azores (Jarvis, 1980). The Azorean plants are morphologically distinct from the Madeiran plants, and the two archipelagos are highly divergent at SSR molecular markers (Borges . Lastly, they form distinct clades with genomic markers (Crawford et al., 2019;Kerbs, 2020), indicating that the Azorean plants likely represent a distinct species. Azorean T. succulenta are known from very few populations in rocky areas and coastal cliffs on several islands (Jarvis, 1980;Schaefer, 2005). The breeding system of Azorean T. succulenta is considered largely SI or PSC based on lack of or very low self-seed set from two populations on two islands in the Azores . Recent mating system estimation using field collected progeny sets from two populations on Madeira confirm predominant outcrossing (Gibson et al., 2020). However, despite a functional SI system, two of the 75 offspring scored in that study were confirmed as products of self-fertilization (Gibson et al., 2020). Crawford et al. (2019) recently documented a breeding system shift from SI to SC within T. succulenta from the Azorean island of Graciosa ( Figure 2). In contrast to Madeiran T. succulenta and populations from other islands in the Azores, plants from Graciosa readily self in the greenhouse. This shift has a genetic basis: selfseed set segregates in a nearly Mendelian fashion in F 2 hybrids between Graciosa plants and SI T. succulenta (J. K. Kelly et al., unpubl). As a contrast to T. succulenta in Madeira, we here estimate the outcrossing rate in the field for the two known populations on Graciosa (GRSC and GRBL; Figure 2). These populations are very small, with estimated census sizes of 30-80 (GRBL) and 10-20 individuals (GRSC), respectively. Judging from the few collections of T. succulenta made from Graciosa, it appears that flowering can occur from July through September, which falls within the broad flowering period (June to September) for SI populations of Azorean T. succulenta (Jarvis, 1980;Schaefer, 2005). Both populations occur in disturbed habitats (Crawford et al., Schaefer, pers. obs.). These two populations are sister and form a strongly supported clade in a molecular phylogeny of Azorean T. succulenta (Crawford et al., 2019). Field studies conducted on Graciosa after this manuscript was accepted indicate that population GRSC may be extinct, as plants could not be located (H. Schaefer, pers. obs. 2020). The floral parts in plants from GRSC and GRBL are smaller than in SI T. succulenta (Crawford et al., 2019; L. Borges Silva & M. Moura, unpubl.) but the "selfing syndrome" (Cutter, 2019;Ornduff, 1969;Slotte et al., 2012) is not nearly as highly developed as in ostensibly more ancient transitions to predominant selfing in Macaronesian Tolpis (Crawford et al., 2008;Koseva et al., 2017;Soto-Trejo et al., 2013).
Another line of evidence supports Graciosa T. succulenta as a recent origin of SC. The transition to selfing was likely associated with the colonization of disturbed habitats (complex volcanic history and/or human colonization) on Graciosa. If so, the loss of SI occurred in the range of 400,000 and 1.05 million years as estimated by several radiometric dating studies (Larrea et al., 2014) or 700,000 years as dated by Sibrant et al. (2014) for the age of Graciosa. Secondly, the estimated divergence time between SI and SC T. succulenta in the Azores based on a dated Bayesian tree using genome-wide genotyping is 511,000 years (B. Kerbs et al., unpubl.). The range of island age estimates from radiometric dating, estimated divergence time from a dated molecular phylogeny, and the limited evolution of the floral selfing syndrome point to a recent origin of SC in Graciosa succulenta (Crawford et al., 2019). We hypothesized that the transition to SC could lead to high selfing in GRSC and GRBL because of the limited number of compatible mates and pollinators in small populations, which should produce strong selection for reproductive assurance (Pannell, 2015). Alternatively, ID could select against selfing and maintain outcrossing.

| Sampling
This study examined the two known populations of T. succulenta A total of 22 progeny (2-7 progeny per mother plant from population GRBL and 20 progeny (2-7 per mother plant) from GRSC were genotyped (Table 1). Mean and range of greenhouse self-seed set in the two populations were 60% (31%-96%) for GRBL and 41.9% (0%-99%) for GRSC (Crawford et al., 2019). Vouchers of progeny are deposited in the R. L. McGregor Herbarium of the University of Kansas (KANU).

| Cultivation and DNA extraction
Seeds from wild maternal plants were germinated and reared in greenhouses at the University of Kansas. Leaf tissue was collected, pressed, dried, and ground. Samples were then frozen using liquid nitrogen and pulverized using chromium beads. DNA was subsequently extracted from the ground, dried tissue using DNeasy Plant Mini Kits (Qiagen Inc.) and DNA quantity was validated using a Qubit fluorometer (Thermo Fisher Scientific).
DNA from samples was cut using the restriction enzyme Csp6I (syn. CviQI), 250-300 bp fragments were selected for using a BluePippin (Sage Science), and 6 bp barcodes were ligated to fragments. DNA was sequenced on an Illumina NovaSeq 6000 (Novogene) to produce 150 bp paired-end reads. Following sequencing, the demultiplexing of FastQ files was carried out using STACKS (Catchen et al., 2013) and loci were de novo assembled using the same program with parameters M = 2, m = 3, n = 1, as well as invoking the deleveraging algorithm and specifying alpha = 0.05 in ustacks. De novo assembly and SNP calling yielded a total of 111,613 variant sites.

| Mating system estimation and multiple paternity
The resultant VCF from STACKS was assessed, and SNPs were filtered using custom python scripts. We first determined the distribution of reads per SNP per individual at all SNPs called in at least 10 plants. SNPs called in fewer than 10 plants were suppressed.
We next eliminated SNPs with excessively high or low coverage: low threshold = 7.4 based on 10th percentile of distribution, high threshold = 42 based on 90th percentile. 8,530 SNPs remained.
We then suppressed SNPs that exhibited a statistically significant excess of heterozygotes (relative to Hardy-Weinberg proportions) and then thinned to the data to one SNP per RADtag. We selected the one SNP per RADtag with the most minor genotype calls. This produced the list of 516 SNPs that were formatted for input to BORICE. We simultaneously made the "CX" file that species the fraction of SNPs called for each plant, an input to the genotype uncertainty calculations in BORICE. We first determined the distribution of reads per SNP per individual at all SNPs called in at least 10 plants. SNPs called in fewer than 10 plants were suppressed.
We next eliminated SNPs with excessively high or low coverage: low threshold = 7.4 based on 10th percentile of distribution, high threshold = 42 based on 90th percentile. 8,530 SNPs remained.
We then suppressed SNPs that exhibited a statistically significant excess of heterozygotes (relative to Hardy-Weinberg proportions) and then thinned to the data to one SNP per RADtag. We selected the one SNP per RADtag with the most minor genotype calls. This produced the list of 516 SNPs that were formatted for input to BORICE. We simultaneously made the "CX" file that species the fraction of SNPs called for each plant, an input to the genotype uncertainty calculations in BORICE. The programs used to perform these operations as well as the genotype file and BORICE settings script are contained in Supplemental File 1.
BORICE was run with burn-in of 1,000 and chain length of 4,000 steps.

| Inbredness of individual parents and offspring
Individual progeny from both populations were determined to be outcrossed or selfed with strong confidence (posterior TA B L E 1 Assignment of sibship for each family across two Tolpis populations probabilities > 0.99). One offspring of the five in family 8 from GRSC was found to be selfed; hence, the posterior probability for the overall outcrossing is 0 at t = 1 ( Figure 3). All other offspring were found to be produced via outcrossing. All eight maternal plants across the two populations were determined to be outbred. The high posterior probabilities-P[IH = 0] = 100% for families 1, 3, 4, 5, and 6, 95% for family 2, 94% for family 7, and 81% for family 8. There was minimal (0.05) allele frequency divergence between populations at the markers used for BORICE.

| Sibships
Across all families, the probability that progeny are full sibs is 15.3%.
No full siblings were detected in five of the eight families (Table 1).
Families 1 and 4 contain one set of full sibs each, and family 5 contains two sets of full sibs (Table 1; Figure 4). Outputs show a very high confidence (>90%) in assignment of offspring to a sire in the vast majority (78%) of contrasts. There was moderate support (~50% to 90% probability) for 23% of contrasts.  Sporophytic self-incompatibility (SSI), which is characteristic of Asteraceae (Crowe, 1954;Gerstel, 1950;Hughes & Babcock, 1950), may appear more restrictive in terms of self-and cross-compatibility than gametophytic self-incompatibility (GSI) systems. In the latter, the haploid genotype (allele at S-locus) of pollen controls compatibility but with SSI the alleles in the parental plant occur on the pollen grain and if either of the pollen alleles are the same as the stigma, fertilization does not occur. However, dominance relationships among S-alleles of SSI in Asteraceae (Crowe, 1954;Gerstel, 1950;Hughes & Babcock, 1950) increases cross-and self-compatibility because recessive alleles are not expressed, and therefore do not prevent fertilization in the presence of more dominant alleles (Brennan et al., 2011(Brennan et al., , 2013Hiscock, 2000). Alleles may also show different dominance relationships in the stigma and pollen (Brennan et al., 2006;Hiscock & Tab ah, 2003). Unlike codominant alleles where an increase in frequency of an allele will result in it finding fewer compatible mates in a population (Byers & Meagher, 1992), more recessive alleles in a dominance hierarchy are not subjected to negative frequency dependent selection and will increase in a population. Thus, species with SSI may set seed both by outcrossing and selfing in small populations despite low S-allele diversity (Brennan et al., 2002;Silva et al., 2016).

| Mating system
The automatic selection hypothesis posits that SC mutations are selected because selfing variants, in contrast to outcrossers, can fertilize their own ovules giving them a 3:2 transmission advantage (Fisher, 1941;Pannell, 2015). The highly outcrossing mating system in the two populations of Azorean Tolpis despite having the ability to set self-seed in the greenhouse indicates that there are negative factors associated with selfing, one of which could be inbreeding depression (ID), see below. Selfing may also be disadvantageous when there is a reduction in pollen and ovules available for outcrossing because they are used for selfing, so-called pollen and seed discounting (Holsinger, 1991;Nagylaki, 1976). The mating system of plants where SC has ostensibly evolved recently may depend on a complex combination of the aforementioned factors (Voillemot et al., 2019).
The role(s) of ID and pollen/seed discounting for the outcrossing mating system in the two SC populations of Azorean Tolpis are not known.

| Inbreeding depression
ID is considered a major factor opposing the evolution of selfing subsequent to SI breakdown (Charlesworth & Willis, 2009). ID within populations can be estimated comparing the inbreeding coefficients of seeds and adults (Ritland, 1990;Scofield & Schultz, 2006)

| Paternity
Pannell and Labouche (2013) Crawford et al., unpubl.). In the present study, seed was collected in bulk from each maternal plant, thus precluding determination of intra-and intercapitular components of multiple paternity.
The percent of full sibs detected in the present study may be compared with other Asteraceae, including the recent study of T. succulenta on Madeira island (Gibson et al., 2020). Those populations had 22% full sibs, some 50% higher than detected in the two small Graciosa populations. Thus, correlated paternity is higher in the SI species in Madeira than in the SC populations on Graciosa. The reason(s) for the differences are obscure and comments would be highly speculative. Hardy et al. (2004) found that in the rare endemic SI perennial herb Centaurea corymbosa (Asteraceae), 20% of sibs were full (same sire), and Sun and Ritland (1998) found 19% full sibs in the SI annual Centaurea solstitials. Perhaps, the important point is that in the small insular populations on Madeira and Graciosa, the number of sires among fruits of maternal plant is comparable not only to the other few Asteraceae investigated but also to among-fruit values estimated for other plant families (Pannell & Labouche, 2013). There is little evidence that particular sires are contributing disproportionately to the progeny of maternal plants, making it unlikely that biparental inbreeding is occurring in the populations.

| Summary and questions for future study
There are two major findings of this study. The first is that two small insular populations are highly outcrossing in nature despite the breakdown of SI. However, families of plants grown from seeds collected in nature yield high self-seed set in the greenhouse. The origin of SC in these two small populations that have ostensibly recently colonized Graciosa would seem to favor the transition to selfing (Pannell, 2015). The reasons for the somewhat unexpected results remain to be determined, and in a real sense, this study raises more questions, than it answers. One likely reason for the retention of outcrossing is high ID (Layman et al., 2017;Pannell, 2015).
That is, selfed seeds are presumably aborted on production, do not germinate well, do not flower well, or are vegetatively noncompetitive. Whether this is the situation for T. succulenta awaits further study; two generations of selfed progeny in the greenhouse flowered and set fruit but more thorough studies are needed. Of course, it may be that fitness of selfed progeny is lower in the natural habitat than under greenhouse cultivation (Arista et al., 2017;Armbruster & Reed, 2006). Self-seed set is sometimes used to infer the breeding system in plants from oceanic archipelagos Bernardello et al., 2001;Chamorro et al., 2012;Crawford et al., 2011) and indeed likely provide useful first estimates of breeding system.  Gibson et al. (2020) discuss the advantages of the methods employed herein for studying the small populations of rare island plants.
A major advantage, with conservation implications, is that mating system can be inferred because individual progeny can be called as selfed or outcrossed with few maternal families and few progeny per family. This advantage can hardly be overstated given the low seed set available for progeny arrays and small population sizes. Although small populations provide challenges for mating system and paternity studies, they may offer certain advantages using genome-wide genotyping. In small populations such as the two examined in this study, it may be feasible to map all plants making it possible to detect genetic structure in the populations (Colicchio et al., 2020) and to infer not just the number of sires but the specific sires of progeny (Gibson et al., 2020).

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The sequence data are deposited in the NCBI Sequence Read Archive (SAMN16395570-SAMN16395611) under BioProject ID PRJNA668054.