Rapid detection of RNase-based self-incompatibility in Lysimachia monelli (Primulaceae)

Received: 7 June 2024 | Accepted: 4 December 2024 DOI: 10.1002/ajb2.16449 RESEARCH ARTICLE Rapid detection of RNase‐based self‐incompatibility in Lysimachia monelli (Primulaceae) Karolis Ramanauskas 1 |Francisco J. Jiménez‐López 2 |Mercedes Sánchez‐Cabrera 2 | Marcial Escudero 2 |Pedro L. Ortiz 2 |Montserrat Arista 2 |Boris Igić 1 1 Department of Biological Sciences, University of Illinois at Chicago, Chicago, IL 60607, USA 2 Departamento de Biología Vegetal y Ecología, Universidad de Sevilla, Apdo. 1095, Sevilla 41080, Spain Correspondence Boris Igić, 840 West Taylor St., M/C 067, Chicago, IL 60607, USA. Email: [email protected] Abstract Premise: Primroses famously employ a system that simultaneously expresses distyly and filters out self‐pollen. Other species in the Primulaceae family, including Lysimachia monelli (blue pimpernel), also express self‐incompatibility (SI), but involving a system with distinct features and an unknown molecular genetic basis. Methods: We utilize a candidate‐based transcriptome sequencing (RNA‐seq) approach, relying on candidate T2/S‐RNase Class III and S‐linked F‐box‐motif‐ containing genes and harnessing the unusual evolutionary and genetic features of SI, to examine whether an RNase‐based mechanism underlies SI in L. monelli. We term this approach "SI detection with RNA‐seq" (SIDR). Results: The results of sequencing, crossing, population genetics, and molecular evolutionary features each support a causal association linking the recovered genotypes with SI phenotypes. The finding of RNase‐based SI in Primulaceae (Ericales) all but cements the long‐held view that this mechanism was present in the ancestral pentapetal eudicot, whose descendants now comprise two‐thirds of angiosperms. It also significantly narrows the plausible maximum age for the heterostyly evolution within the family. Conclusions: SIDR is powerful, flexible, inexpensive, and most critically enables work in often‐neglected species. It may be used with or without candidate genes to close enormous gaps in understanding the genetic basis of SI and the history of breeding system evolution. KEYWORDS breeding systems, molecular evolution, RNA‐seq, self‐incompatibility, SIDR INTRODUCTION An overwhelming majority of flowering plants are co‐ sexual, simultaneously expressing both male and female gametes. Despite the resulting potential for self‐ fertilization, they predominantly or exclusively mate with other individuals. The most widespread and commonly employed contrivance usedbyplantstoensurecross‐ fertilization is self‐incompatibility (SI; de Nettancourt, 1977). Generally defined, SI is any mechanism of recognition and rejection of a plant's own pollen, employed between pollination and fertilization, that prevents zygote formation (Stout, 1917). SI mechanisms affect the rules of mating and modulate the amount and distribution of genetic variation, with reverberating effects on selection across levels of biological hierarchy (Holsinger, 2000). Indeed, obligate outcrossing may be a key feature leading to diversification rate differences among angiosperms, allowing lineages increased persistence and diversification (Goldberg et al., 2010; de Vos et al., 2014). Characterization of SI mechanisms is additionally significant because their expression may strongly Am J Bot. 2025;112:e16449. wileyonlinelibrary.com/journal/AJB | 1of16 https://doi.org/10.1002/ajb2.16449 This is an open access article under the terms of the Creative Commons Attribution‐NonCommercial‐NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non‐commercial and no modifications or adaptations are made. © 2025 The Author(s). American Journal of Botany published by Wiley Periodicals LLC on behalf of Botanical Society of America. affect yield in a variety of crops and generally alter the propensity for domestication (Muñoz‐Sanz et al., 2020). In the context of plant conservation, information about the genetic basis of SI may be employed to alleviate mate limitation following severe population decline (DeMauro, 1993; Weekley et al., 2002). SI is present in over 100 plant families of angiosperms (Raduski et al., 2012). Conservative estimates suggest that it has independently evolved dozens of times (Bateman, 1952; Igićet al., 2008), not including analogous mechanisms in fungi and hermaphrodite animals (Hicks et al., 1979; Harada et al., 2008). For example, Brassicaceae (Nasrallah et al., 1985) and Papaveraceae (Foote et al., 1994) express distinct SI mechanisms, underlain by unrelated genetic components (Figure 1). They have repeatedly convergently accomplished the task of self‐recognition and rejection utilizing distinct molecular pathways. Yet, a comparable weight of evidence supports the ancestry of a single homologous SI mechanism, broadly distributed across eudicots (Igićand Kohn, 2001; Steinbachs and Holsinger, 2002; Ramanauskas and Igić,2017). Many species, including unrelated and disparate ones, such as those in Cactaceae, Rutaceae, and Solanaceae, employ ribonucleases to enforce SI (Anderson et al., 1986;Liangetal.,2020; Ramanauskas and Igić,2021). This mechanism is commonly referred to as “RNase‐based SI”(RSI). Individuals with RSI express an orthologous gene, a subtype of class III T2/S‐RNases, in their stylar transmitting tracts (McClure et al., 1989; Sassa et al., 1993). Their pollen tubes each express many paralogous tandem‐duplicated F‐box‐motif‐containing genes (S‐locus F‐ box genes or SLFs), which interact with S‐RNases to cause SI (Lai et al., 2002;Sijačićet al., 2004;Kuboetal.,2010;Liang et al., 2020). Superficially, the mechanism appears fairly simple. Pollen, which expresses the contents of one SI locus allele (S‐allele or S‐haplotype), lands and grows through diploid style, which expresses two S‐alleles. If the haploid pollen grain contains an S‐allele (e.g., S1 ) matching either of two (e.g., S S 1 2 ) foundinthediploidstyle,fertilization is prevented and no seed is made. Class III T2/S‐RNases are absent in monocots and magnoliid dicots (Ramanauskas and Igić,2017). Therefore, RSI is widely thought to have originated in the common ancestor of eudicots and to have been continually present in many descendant lineages (Igic and Kohn, 2001;Steinbachsand Holsinger, 2002; Charlesworth, 2003; Allen and Hiscock, 2008; Niu et al., 2017;Liangetal.,2020; Fujii and Takayama, 2020). The apparent simplicity of this mechanism collides with an astonishingly complex reality. The system instead accomplishes self‐rejection by non‐self recognition (Kubo et al., 2010; McClure et al., 2011; Fujii et al., 2016). By default, S‐RNases enter growing pollen tubes and degrade their rRNA (McClure et al., 1989), disabling the pollen grain's journey to fertilization. However, if S‐RNases are recognized by any of a number of pollen‐expressed SLFs located on non‐self S‐haplotypes, the S‐RNases are inactivated and degraded via a ubiquitin‐mediated pathway (Kubo et al., 2010). SLFs physically linked to a specificS‐ RNase on the same S‐haplotype (self) are unable to defuse them. Haplotype‐dependent degradation relies on the interaction between the S‐RNase, expressed from a particular S‐haplotype, and one or more non‐self SLFs, expressed on different S‐haplotypes (Kubo et al., 2010). SLFs form a part of multi‐subunit E3 ubiquitin ligase complexes, alongside Skp1, Cullin1, and Rbx1 (Li et al., 2016). Non‐self S‐RNases are therefore ubiquitinated and degraded by the 26S proteasome, and only an unrecognized self S‐RNase prevents self‐pollination. The details of the RSI mechanism are daunting, remain partly unclear, and include significant evolutionary questions about the initial function and assembly of new S‐ alleles, particularly regarding the plausible origin of the system (Uyenoyama, 2000; Fujii et al., 2016). In species with a functional RSI system, each S‐haplotype consists of a single S‐RNase tightly linked to dozens of SLF genes (Qiao et al., 2004; Kubo et al., 2010). Rare‐allele advantage at the S‐locus is often so strong that even a modestly sized population can contain dozens of S‐haplotypes (Lawrence, 2000), but the system cannot function with fewer than three S‐haplotypes (Wright, 1939). Finally, at least in the genus Prunus (Rosaceae), the operation of the system departs considerably from others where it is found (Tao and Iezzoni, 2010; McClure et al., 2011; Fujii FIGURE 1 A summary of the known distribution and history of RNase‐based SI (RSI) systems in selected angiosperms. Tip colors indicate the type of molecular mechanisms found in each family, with key molecular components as identified in the legend. Colored branches mark a plausible evolution ordering of SI systems. The families Rubiaceae, Primulaceae, and Passifloraceae are shown as pied, with some species in each family employing homomorphic and independently derived heteromorphic SI. (The genetic bases of heterostyly in Rubiaceae are unknown. Perhaps two uncharacterized mechanisms of SI in Passifloraceae, as well as their likely evolutionary history, also remain unknown.) The detection of new instances of RSI across the core eudicots, and its presence in the common ancestor of two‐thirds of all land plants, increasingly places the weight of evidence on the homology of this mechanism. 2of16 | RAPID DETECTION OF RNASE‐BASED SELF‐INCOMPATIBILITY IN LYSIMACHIA MONELLI (PRIMULACEAE) 15372197, 2025, 1, Downloaded from https://bsapubs.onlinelibrary.wiley.com/doi/10.1002/ajb2.16449 by Readcube (Labtiva Inc.), Wiley Online Library on [31/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License et al., 2016), potentially as a byproduct of paralogous locus copies (Aguiar et al., 2015; Zhao et al., 2022a), as is the case elsewhere, in the Brassicaceae SI system (Chantha et al., 2017), and perhaps expected given opportunity for divergence across 100 million years time scales and thousands of lineages. Nevertheless, a vast suite of characteristics, including the evolutionary relationships between the components of the RSI mechanism, are consistent with homology (reviewed in Ramanauskas and Igić,2017). This hypothesis is rapidly solidifying, as evidenced by the successful candidate gene‐based detection of RSI across distant flowering plant orders (Xue et al., 1996; Asquini et al., 2011; Liang et al., 2020; Ramanauskas and Igić,2021). One clear prediction is that RSI ought to be found in many other species across eudicots, particularly those that display similar superficial features. An interesting and generally underappreciated group of candidates for the expression of RSI is species in the primrose family (Talavera et al., 2001). At least two unrelated mechanisms are thought to co‐ occur within this family (Primulaceae): one consistent with RSI and another entirely different (Gibbs and Talavera, 2001). Species in the genus Primula, famously express the latter kind of mechanism, displaying heteromorphic SI, which couples heterostyly and self‐pollen rejection (Darwin, 1877; Richards, 2003). Heterostylous species produce two kinds of individuals, or morphs, distinguished by reciprocal positioning of stigmas and anthers in their flowers. One morph has long styles and the other short styles. Matings between different morphs are compatible, but self‐and within‐morph matings result in reduced fruit and seed sets (Darwin, 1877; Richards, 2003). Recent work has shown that a cause of SI in Primula involves a brassinosteroid‐ inactivating cytochrome P450 gene (CYP734A50 or S‐CYP; Huu et al., 2022), which is also involved in determining the relative pistil and anther position (Huu et al., 2016; Li et al., 2016). The S‐CYP mechanism appears to be relatively young, confined to primroses (de Vos et al., 2014), and unrelated to perhaps two dozen other occurrences of heterostyly and SI across angiosperms (Weller, 2009; Matzke et al., 2021; Potente et al., 2022). Other species in the family, including those in the genus Lysimachia, appear to express SI that bears close resemblance to the phenotypic and hereditary characteristics of RSI. However, the molecular genetic basis underlying their SI response is presently unknown. We have recently shown that a candidate gene screen, which combines an experimental design harnessing many unusual features of SI along with phylotranscriptomic and biosystematic approaches, can be successfully employed to uncover the genetic basis of SI (Ramanauskas and Igić,2021;Table1). This powerful, simple, and flexible approach for SI detection with RNA‐seq (SIDR) may be used with or without candidate genes to close enormous gaps in understanding the genetic basis of SI and history of breeding system evolution. Although simultaneously less laborious and inexpensive than filial mapping, functional molecular characterization, or in vitro confirmation, SIDR has the distinct advantage of permitting the discovery of mechanism(s) in non‐model systems. It allows the genetic basisofSItobeinterrogatedinspecieswithdelayedmaturation time and limited genomic resources. Here, we employ SIDR to uncover the molecular genetic basis of SI in Lysimachia monelli (Primulaceae). Concentrating on the single‐copy class III T2/S‐RNases, we examine whether they possess the key features of genes involved in RSI, including heterozygosity, high polymorphism, tissue‐specific, and high expression in pistils, as well as membership in a subgroup of the T2/S‐RNases family with conserved features of RSI shared across eudicots. We use a series of controlled hand pollinations to confirm that a subset of putatively identified S‐allele sequences can predict the pollen‐style interaction, demonstrating a key link between genotype and phenotype. We additionally attempt to recover RSI genes in Primula to piece together and account for their fates, as well as more precisely infer the history of breeding system evolution. We find strong support for RSI as the cause of single‐locus gametophytic SI in L. monelli. Along the way, we demonstrate the value of our approach (Ramanauskas and Igić,2021), whose application can rapidly establish the genetic basis of SI across angiosperms, whether or not species employ RSI. Finally, we consider the implications of these findings in the context of recent work on the genetics of heteromorphic SI (Huu et al., 2020; Matzke et al., 2021; Zhao et al., 2022b), which jointly supports dozens of independent cases of the re‐evolution of heteromorphic SI across eudicots. MATERIALS AND METHODS Study organism and materials L. monelli (L.) U. Manns and Anderberg (= Anagallis monelli L.) is a short‐lived, perennial, diploid species distributed around the western Mediterranean area (Gibbs and Talavera, 2001). In selecting this species for a closer examination and characterization of the genes that might cause SI, we relied on the previous work by Gibbs and Talavera (2001) and Talavera et al. (2001). They established the genetic mode of action as gametophytic—meaning the pollen haplotype is recognized and rejected by the pistil genotype—and uncovered styles as the site of pollen inhibition. This closely resembles the phenomenology observed in other instances of RSI. In the present study, we use a species‐wide sample of 19 individuals from natural populations, along with eight greenhouse‐grown seeds from an additional maternal accession, for a total of 27 individual plants (Table S1). The greenhouse‐grown individuals (all SI) are at least half‐sibs, sharing the maternal parent, and possibly full sibs—also sharing the paternal parent. All samples from natural populations consisted of pistil tissues only. We also collected pistils for each of the greenhouse individuals and added samples of anther and pollen tissue (five individuals) and leaf and root tissue(oneindividual),astissue‐specific expression controls. RAPID DETECTION OF RNASE‐BASED SELF‐INCOMPATIBILITY IN LYSIMACHIA MONELLI (PRIMULACEAE) | 3of16 15372197, 2025, 1, Downloaded from https://bsapubs.onlinelibrary.wiley.com/doi/10.1002/ajb2.16449 by Readcube (Labtiva Inc.), Wiley Online Library on [31/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License RNA extraction and sequencing Two different protocols were used for short‐term storage of the tissues, depending on whether they were field‐collected or growninthegreenhousesattheUniversityofIllinoisat Chicago. Pistils (without ovaries) from 19 field‐collected samples—labeled 148‐19a through 165‐19a (Table S1)—were submerged in 1.5 mL of RNAlater solution (Invitrogen). Pistils (without ovaries), pollen, leaf, and root tissues collected from the greenhouse (labeled 17K68 through 19K02) were placed in 1.5 mL tubes and stored on dry ice. Both field and greenhouse‐ collected samples were moved to a ∘ − 80 C freezer for long‐ term storage. Approximately 100 mg of tissue was ground to a fine powder in 1.5 mL tubes submerged in liquid nitrogen. Total RNA was isolated using the Total RNA Mini Kit (Plant Kit, Cat. No. IB47341, IBI Scientific, Dubuque, Iowa) following the manufacturer's instructions. We assessed RNA concentration and purity with a NanoDrop Lite Spectrophotometer (Thermo Fisher Scientific, Waltham, Massachusetts). The 34 samples used in this study were sequenced as part of a larger sequencing effort, which included four separate sequencing runs and additional samples from other plant species. Sequencing libraries were prepared using the KAPA Stranded mRNA‐Seq (Roche, Basel, Switzerland), and these libraries were sequenced on a single lane of the Illumina (San Diego, California) HiSeq 4000 or Illumina NovaSeq 6000 platform (paired‐end 150 bp reads) at the Duke University Center for Genomic and Computational Biology. Sequence assembly and analyses Raw paired‐end Illumina reads were processed with Rcorrector v1.0.4 (Song and Florea, 2015) to correct sequencing errors, trimmed with Trimmomatic v0.39 (Bolger et al., 2014) to remove low‐quality reads, and filtered with Kraken 2 (Wood et al., 2019) to remove Small and Large Subunit ribosomal RNA (SILVA database) (Quast et al., 2013) and contaminating reads (minikraken2_v2 database). We used custom‐built databases, derived from RefSeq libraries to minimize the number of non‐plant and non‐nuclear reads. Filtered reads were combined across samples into a single RNA‐seq dataset. We conducted a de novo transcriptome assembly using Trinity v2.8.5 (Grabherr et al., 2011)togenerate a single reference assembly. We used default options except for SS_lib_type RF and min_contig_length 200. Transcript quantification was performed using RSEM (Li and Dewey, 2011), and the Trimmed Mean of TABLE 1 Rationale for self‐incompatibility detection with RNA‐seq (SIDR). The overall strategy for a systematic and rapid candidate‐based search for the molecular causes of SI, as outlined in Ramanauskas and Igić(2021), relies on collecting a weight of supporting evidence matching the expectations of functional RNase‐based SI. Key features include the establishment of functional RSI, obligate heterozygosity of S‐genotypes (Figure 2), appropriate tissue expression of candidate gene(s) (Figure 2), genotype–phenotype matching (Figure 3), ancestral polymorphism maintained by negative frequency‐dependent selection (Figure 4), and a close phylogenetic relationship of candidate gene(s) with others known to cause SI, strongly supporting homology (Figure 4). Genetic transformations are the gold standard for functional study. They remain impractical, however, as tools for the initial discovery of new SI mechanisms and assessment of the phylogenetic distribution of previously characterized mechanisms in the majority of angiosperm species. Feature Expectation Outcome Self‐incompatibility (SI) Pre‐zygotic rejection by a dedicated molecular mechanism results in predictable pollen tube arrest following self‐pollinations and some cross‐pollinations ✓ Gibbs and Talavera (2001) and Talavera et al. (2001) found strong SI response (low fruit set), stylar pollen arrest, and genetic pattern matching gametophytic SI Obligate heterozygosity Functioning SI system prevents homozygosity at S‐loci, yielding two distinct alleles per individual ✓ All plants contain heterozygous candidate S‐RNases Stylar expression of class III T2/S‐ RNases S‐RNases are generally most highly expressed in the top third of styles ✓ Two copies of style‐expressed (among top genes) class III T2/S‐RNases Pollen expression of S‐linked F‐box family Dozen(s) of SLFs are expressed in pollen, linked in a haplotype with an S‐RNase ✓ Many copies of pollen‐expressed SLFs (phasing within haplotypes awaiting genome sequence) Putative S‐locus genotypes predict crossing phenotypes Functional link with S‐genotypes predicting crossing behavior (phenotype); expect the failure of matching genotypes, the success of non‐matching genotypes ✓ Crossing among sibs with matching genotypes fails, whereas crossing with unmatched genotypes succeeds Ancestral polymorphism Negative frequency‐dependent selection causes extremely old and very high allelic polymorphism at the S‐locus ✓ Scandalous polymorphism, with species‐wide S‐ allele number likely N> 150 Close phylogenetic relationship and shared features with other S‐RNases Despite the size of the T2/S‐RNase family, functional S‐ RNases are expected to be closely related and share many features ✓ Putative alleles recovered group with class III T2/S‐ RNases show features of functional S‐RNases Mutant‐generation (transformation) Molecular knock‐out/knock‐in experiments to confirm functional link; practical in a limited subset of plant species (Not performed) 4of16 | RAPID DETECTION OF RNASE‐BASED SELF‐INCOMPATIBILITY IN LYSIMACHIA MONELLI (PRIMULACEAE) 15372197, 2025, 1, Downloaded from https://bsapubs.onlinelibrary.wiley.com/doi/10.1002/ajb2.16449 by Readcube (Labtiva Inc.), Wiley Online Library on [31/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License M‐values (TMM, mean of log‐expression ratios (Robinson and Oshlack, 2010)) was calculated using SARTools v1.7.3 (Varet et al., 2016) with edgeR v3.30.2 in R v4.0.2 (Robinson et al., 2010; R Core Team, 2020). S‐RNases are the best‐understood portion of the RSI recognition and rejection response. We consequently concentrated on studying the expression of this key component of RSI. We identified candidate S‐RNases and control/reference genes (AC: actin, EF: Elongation factor 1‐ α 1oreEF1a1,andTU: α ‐tubulin) in transcriptome assemblies with the KAKAPO v0.7.9 pipeline (Ramanauskas and Igić,2023). The initial identification of putative S‐ RNases relied on high similarity to known S‐RNases, as well as their sequence features and phylogenetic relationships. We calculated isoelectric point values using the ProteinAnalysis tool in Biopython v1.79 (Cock et al., 2009). Sequences were initially classified as putative S‐RNases if they were heterozygous, contained an active‐ site histidine, showed basic pI values, and clustered within the T2/S‐RNase class III clade. It is not unusual to recover remnant sequences, particularly of S‐RNases, as molecular relicts in genomes that once contained RSI (Ramanauskas and Igić,2017). We attempted to find such relictual S‐RNase sequences with BLAST searches of existing assemblies from the published draft genomes of the Primulaceae species: Primula veris (accession GCA_000788445.1), P. vulgaris (GCA_00 1077355.1), and Embelia ribes (GCA_001753 735.1), which do not express RSI. Phylogenetic analyses Additional T2/S‐RNase and F‐box domain‐containing sequences were obtained using BLAST by searching GenBank nucleotide collection (nr/nt) database and the genomes of Arabidopsis thaliana (Brassicaceae), Coffea eugenioides (Rubiaceae), Prunus avium (Rosaceae), Spinacia oleracea (Amaranthaceae), Solanum chilense (Solanaceae), and Stenocereus thurberi (Cactaceae), downloaded from NCBI. The genome of Citrus maxima (Rutaceae) was downloaded from the Citrus Genome Database website (https://citrusgenomedb. org;Wangetal.,2017). Sequence alignment T2/S‐RNase sequences were translation‐aligned using MAFFT v7.474 (Katoh and Standley, 2013), with globalpair algorithm, and BLOSUM 45 matrix. The alignment was trimmed to include conserved regions 1 through 5 and mapped back to nucleotide sequences. F‐box domain‐containing sequences were translation‐ aligned using MAFFT v7.474, with genafpair algorithm, and BLOSUM 45 matrix. The alignment was trimmedtoexclude5'and3'regionswithcoveragelower than 100 sequences. Phylogenetic tree inference T2/S‐RNase gene tree was inferred with MrBayes version 3.2.5 (Ronquist et al., 2012). The analysis consisted of 20 independent MrBayes runs. Each run had one cold and three heated chains. We implemented a general time‐ reversible (GTR) nucleotide substitution model with Gamma‐distributed among‐site rate variation and four Gamma rate categories. The runs were allowed to complete 300 million generations, with trees and parameters sampled every 10,000 generations. We used the last 50 million generations as the posterior distribution draw, which consisted of 100,000 samples: 50,000,000 generations × 1sample 10,000 generations ×20 runs. The sampled trees were summarized using the MrBayes sumt command. F‐box domain‐containing gene tree was inferred with IQ‐TREE version 2.2.0 (Nguyen et al., 2015) using the GTR model of nucleotide substitution; base frequencies were optimized by maximum‐likelihood -m GTR+FO. Branch support was estimated using the ultrafast bootstrap approximation -bnni -B 10,000. For clarity, we display a subsampled T2/S‐RNase gene tree (Figure 4) and place a full version in the Supplementary Information (Figure S1), including a separate figure for SLFs (Figure S2). Controlled crosses: tests of genotype–phenotype association We performed controlled greenhouse crosses to confirm the association between an individual plant's S‐allele function (recognition phenotype) and S‐RNase sequences (genotype). For gametophytic SI (Gibbs and Talavera, 2001;Talavera et al., 2001), we expect each plant to express two S‐RNase alleles (S‐alleles) in pistils and an unknown complement of F‐box genes in pollen. Individuals with matching S‐alleles ought to be cross‐incompatible and all other individuals cross‐compatible. For example, we expect crosses with a plant bearing the recovered RNA‐seq S‐genotype S S 1 2 to be reciprocally sterile with each matching genotype and reciprocally fertile with other genotypes. Crossing patterns consistent with this expectation constitute evidence supporting the involvement of S‐RNase genotypes in causing pistil‐side SI. The only other explanation for a matching pattern would involve linked loci in the same gene family, but linked loci are not expected to show obligate heterozygosity. A full diallel crossing design was not possible because of the limited number of co‐flowering individuals. Nevertheless, we self‐ fertilized each plant and conducted as many cross‐pollination as practical—often using all flowers for either crosses or tissue collection. RAPID DETECTION OF RNASE‐BASED SELF‐INCOMPATIBILITY IN LYSIMACHIA MONELLI (PRIMULACEAE) | 5of16 15372197, 2025, 1, Downloaded from https://bsapubs.onlinelibrary.wiley.com/doi/10.1002/ajb2.16449 by Readcube (Labtiva Inc.), Wiley Online Library on [31/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Species‐wide S‐allele number Spectacular polymorphism is a hallmark of long‐term negative frequency‐dependent selection at S‐loci (Richman et al., 1996; Paape et al., 2008). We expected a large number of S‐alleles, as found in other species with RSI. We used a maximum likelihood (ML) estimator N ˆ(Paxman, 1963) obtained by equating E n n ()= and numerically solving                                E nN N ()= 1− 1−2, m where nis the number of alleles in a sample of mgenotypes. We computed an approximate 95% confidence interval by finding the values of Nfor which NN l ik( ) ≥ lik( ˆ)−1.9 2 . The 95% confidence interval consists of all values of Nfor which the log‐likelihood function drops offby no more than 1.92 units below the value of N ˆ. We obtained two sets of estimates for total S‐allele numbers. The first estimate used 21 plants, including all 19 wild‐collected samples and two (maternal and paternal) genotypes that could be inferred from the greenhouse‐raised population, for a total of m=2 1 genotypes and n=3 7 alleles. We could not infer the full paternal genotype of the greenhouse individual 17K75, which expressed a unique allele S 3 4 . The second estimate did not use the inferred S‐ alleles of parents from our greenhouse‐raised population, therefore totaling m=19 genotypes and m=3 4 alleles. Detection and estimation of diversifying selection Another characteristic of long‐term negative frequency‐ dependent selection (imposed by functional RSI) is a detectable excess of non‐synonymous substitutions at S‐ RNase codon sites that encode the recognition domain (Clark and Kao, 1991). Nucleotide changes that result in amino acid substitutions are disproportionally retained at residues that participate in interaction with pollen‐ determinant genes. We implemented a series of models in PAML v4.8 (Yang, 2007) to detect the signature of diversifying selection on the putative S‐RNases in L. monelli. Sites with values of ω< 1 indicate an excess of synonymous substitutions and therefore can be inferred to evolve under purifying selection that constrains certain nucleotide changes at those sites. Alternatively, those sites for which a significant excess of non‐synonymous substitutions is recorded, or ω> 1 , suggest molecular evolution under positive selection pressure. Values of ω≈ 1 are inferred as consistent with neutrality. We estimated non‐synonymous to synonymous rate ratio ( ω ) and calculated posterior probabilities that codon sites are under “positive”selection using CODEML in PAML v4.8. We compared values of model log‐likelihoods in two nested tests, whose nomenclature follows (Yang, 2007). We report the results of likelihood ratio tests (LRTs) and consider a strong preference for models with positive selection to indicate support for the operation of negative frequency‐ dependent selection at the sequenced S‐alleles. We considered codons to be persuasively likely to have experienced diversifying selection if they had a significant posterior probability ( p p>0. 9 ) of having been drawn from the positively selected class ( ω> 1 ; Yang, 2007). RESULTS RNA‐seq and transcriptome assembly The number of resulting read pairs ranged from 3,323,529 to 10,803,599 with a median of 6,171,496 and an average of 6,692,478 (Table S1). Sequencing of L. monelli pistil mRNA resulted in a median transcriptome assembly of 14,620 genes and 39,407 isoforms (range 6391–48,865 and 12,641–119,138, respectively; detailed in Table S1). Pistil expression of S‐RNases Using minimally invasive sampling, relying on as few as two pistils per individual, we were able to recover 38 unique candidate S‐RNase sequences. They have coding region sequences from 705 to 858 bp, similar to known S‐RNases from other plant families. We recovered exactly two putative S‐alleles per individual in all sequenced individuals. In total, we found 38 putative unique S‐RNase sequences in 27 individuals (genets). The mean and median coverage for each assembled S‐allele sequence was high, ca. 3000x. We carefully examined each of the 38 sequences to rule out the possibility of mischaracterization of paralogs, chimeras, or otherwise artifactual gene assemblies in several ways. We checked all individual reads in assemblies and retained only those with one or no differences. Because S‐allele divergences are large, the occurrence of intra‐locus chimeras is relatively easily spotted, characterized by high‐identity and low‐identity regions. Additionally, we inspected all recovered class III T2/S‐RNase sequences for other artifacts, and unusual variation in coverage, identity, and features. A single isoform was recovered for all coding region sequences (occasionally, we recovered immature mRNA, with non‐spliced introns, designated as an isoform). Eight putative S‐alleles were recovered in more than one individual (different genets). In each case, coding region sequences were identical or extremely highly conserved ( > 99. 5% nucleotide identity), indicating the absence of assembly artifacts. All of the putative S‐RNase sequences display the features of functional S‐RNases (Ramanauskas and Igić,2017). First, each individual is heterozygous at the inferred S‐locus. Two unique putative S‐allele sequences are present in each individual's pistils, and the remaining 36 are absent (Figure 2). This strongly suggests that these are products of a single 6of16 | RAPID DETECTION OF RNASE‐BASED SELF‐INCOMPATIBILITY IN LYSIMACHIA MONELLI (PRIMULACEAE) 15372197, 2025, 1, Downloaded from https://bsapubs.onlinelibrary.wiley.com/doi/10.1002/ajb2.16449 by Readcube (Labtiva Inc.), Wiley Online Library on [31/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License diverse S‐locus, and not unrelated paralogs, which would be present in each individual, or minimally cause the inference of additional “alleles.” Second, each sequence is highly expressed in pistils and absent in other tissues (Figure 2). Combining both putative S‐locus products, or S‐alleles, the expression is very high indeed (TMM = 7276.4): 36‐fold, 10‐fold, and 37‐fold higher than actin (TMM = 204.6); αEF1 1 (TMM = 741.8); and tubulin (TMM = 196.8), respectively. Third, each sequence contains highly conserved residues (Vieira et al., 2008), closely resembling the structure of S‐ RNases. This includes active‐site histidines present across the T2/S‐RNase protein family, multiple predicted disulfide bond‐ forming cysteine residues at conserved sites, and a highly basic isoelectric point ( p I ¯=8.7 1 , range [7.72, 9.18 ] ), similar to the known S‐RNases expressed in other families with RSI (Table S2 and Figure S6; Ramanauskas and Igić,2017). Finally, sequencing of occasional immature mRNAs allows us to detect the location of a single intron, which is present at the same position and in the same reading frame (intron phase) as is canonically found in functional S‐ RNases of other euasterid RSI representatives (Igićand Kohn, 2001; Ramanauskas and Igić,2017). Other members of this protein family in our samples did not satisfy any of the criteria mentioned earlier, and their expression is much lower than that of candidate S‐RNases. Pollen expression of putative SLF genes We found 18 candidate SLFs, which we infer to be S‐locus‐ related principally due to their expression in appropriate tissue (anthers) and their phylogenetic placement (Figure S2). Compared to the single‐or low‐copy RNase family genes, the F‐box‐containing gene family is enormous, with hundreds of members in plant genomes (Yang et al., 2008) and possibly rampant gene conversion obscuring their phylogenetic relationships (Newbigin et al., 2008; Kubo et al., 2015). It is difficult to clearly discern individual S‐linked SLFs in expression analyses because their sequences are not as distinctive as those of S‐RNases (Kubo et al., 2010; Wu et al., 2020). Therefore, it is not clear whether the recovered S‐RNases are linked to (m)any candidate SLFs. Generally, a single S‐RNaseislinkedtodozensofSLFs within a single S‐haplotype (Kubo et al., 2015;Wuetal.,2020), FIGURE 2 Geographic origin and expression of putative S‐RNases in Lysimachia monelli. (A) Map of sampled individuals' population locations (green circles) around the Western Mediterranean, in relation to the known occurrences (dark circles). (B) Illustration of the study species and tissues used in transcriptome sequencing. (C) Expression scale, colored as shown, giving log 2 (TMM + 1). (D) Putative S‐allele (S‐RNase) expression for individual plant accessions with sequenced transcriptomes shown in rows and alleles shown in columns, which are numbered from 1 to 38. Each sequenced plant is heterozygous, expressing two alleles. The expression of each allele pair, when present, is very high, among the top 1% of pistil‐expressed genes. Conversely, in each case, none of the remaining 36 alleles are found, and their expression cannot be detected in root, leaf, and pollen, where examined (top). Subsequent crossing experiments, conducted on eight individuals (accessions with “K”)confirm that at least a subset of 38 putative S‐alleles in fact comprise S‐genotypes, which predict crossing phenotypes (success/failure; Figure 3). Expression patterns of three control genes are also shown: AC, actin; EF, elongation factor EF‐1alpha; TU, alpha‐tubulin. RAPID DETECTION OF RNASE‐BASED SELF‐INCOMPATIBILITY IN LYSIMACHIA MONELLI (PRIMULACEAE) | 7of16 15372197, 2025, 1, Downloaded from https://bsapubs.onlinelibrary.wiley.com/doi/10.1002/ajb2.16449 by Readcube (Labtiva Inc.), Wiley Online Library on [31/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License and distinct haplotypes may share a similar complement of SLFs. Each obligately heterozygous plant expresses two S‐ haplotypes, so many of these are identical in a given heterozygous genotype. Taken together, this makes ascertainment of a particular SLF candidate as belonging to individual S‐ haplotypes very difficult, even with deep chromosomal‐scale genome sequencing, because misassembly is common (Kubo et al., 2015; Ramanauskas and Igić,2017;Wuetal.,2020). The expression of each 18 candidates in pollen is, however, both consistent with—and indicative of—a role in a functional RSI system. S‐RNase genotypes predict crossing phenotypes Confirmation of an association between putative S‐RNase genotypes and the predicted crossing phenotypes is critical for inferring their role in causing RSI. L. monelli employs a gametophytic SI mechanism (Talavera et al., 2001). Therefore, individual plants ought not to be able to produce fruit when self‐fertilized or cross‐fertilized with other individuals, whose S‐genotypes match (e.g., dam S S 1 2 xsire S S 1 2 fails to set fruit or seed). Crosses among S‐genotypes with single matching alleles (e.g., dam S S 1 2 xsire S S 1 3 ) generally result in a complete fruit set, indistinguishable from fully compatible crosses, which involve unmatched alleles (e.g., dam S S 1 2 xsire S S 34 ). As expected, in each case, crossing results are consistent with the expectation that the S‐RNase sequences we obtained represent functional pistil‐part S‐alleles in an operational single‐locus gametophytic SI (Figure 3). All self‐fertilizations failed to produce fruit, as did all crosses involving SI individuals with fully matching S‐allele genotypes (Figure 3). This includes four S S 14 1 6 genotypes and two S S 18 2 1 genotypes. In contrast, crosses involving individuals not sharing full genotypes were successful, and all such combinations yielded fruit sets (not rejected in half‐ compatible crosses). Therefore, we provide a strong confirmation of linked function for these four alleles ( S SS S,,,and 14 16 18 2 1 ). We were unable to directly test the remaining putative S‐alleles, from wild collections, whose identity is only inferred from sequence features, expression data, and obligate heterozygosity. One allele tested in crossing experiments, S 1 6 , was independently recovered in a plant from Portugal, more than 300 km away from our Spain source population (Figure 2D and Table S1). For brevity, relying on heretofore accumulated evidence, we omit the qualification “putative”in the subsequent references to the identified S‐RNases and S‐alleles. Phylogenetic relationships among S‐RNases and putative SLFs Information regarding relationships among S‐RNases appears reasonably well‐preserved, despite the strong selection pressure and enormous time scales involved (Figure S1). For example, we can glean a functional and structural distinction between three “classes”of T2/S‐RNase protein family members. This assessment is further supported by comparisons with data independent of DNA sequence alignments used in phylogenetic inference, such as intron presence, phase, and mature protein isoelectric points (Ramanauskas and Igić,2017). Our sample of S‐RNases belongs to class III T2/S‐RNases, which share conserved structural and biochemical characteristics (Table S2). All genes in this clade sequenced to date are closely related, although not strictly monophyletic (Figures 4and S1). The inference is complicated by the presence of frequently generated relicts, paralogs and pseudogenes, loss‐of‐function mutants, as well as the deep divergences between S‐RNase alleles (Figures 4and S1). We also recovered a group of genes, members of T2/S‐RNase family class III, closely related to S‐RNases, which we labeled as C3A‐and C3B‐RNases. These genes are closely related to S‐RNases, and some are even expressed in styles. Unlike S‐ RNases, they are expressed either in all or most individuals, and they are not associated with S‐phenotypes in our controlled crosses. Moreover, they are recovered at far lower expression levels (on average 22‐fold to 1142‐fold lower). Although plant genomes only contain a dozen, or so, T2/S‐RNase family genes (Ramanauskas and Igić,2017), they may contain several hundred genes with an N‐terminal F‐box motif (Xu et al., 2009). This necessarily complicates FIGURE 3 Crossing relationships among individuals with RNA‐seq identified genotypes. Individual genotypes are given for each parent, using the same S‐allele naming as shown in Figure 2. Rows show dams (pistil parents) and columns show sires (pollen parents), in the format of S#/S#. Immediately below, we list the unique trailing digits of plant accession. Successes are check‐marked, whereas failures are cross‐marked. Empty squares were not attempted and are lightly shaded to indicate expectations. Self‐fertilization attempts are indicated with a deep red color across the diagonal. Crossing between dam × sire combinations involving identical genotypes, shown in red color, ought to result in fruit set failure (pollen alleles are recognized and rejected if they match alleles in the maternal genotype). Otherwise, when the parental plants possess non‐identical genotypes (green), crosses ought to be successful. In each case, self‐ fertilizations fail and the observed crossing data match the expected genotype–phenotype mapping predictions, indicating that the identified sequences involved in controlled crosses most likely correspond to functional S‐alleles, causing SI in this species. 8of16 | RAPID DETECTION OF RNASE‐BASED SELF‐INCOMPATIBILITY IN LYSIMACHIA MONELLI (PRIMULACEAE) 15372197, 2025, 1, Downloaded from https://bsapubs.onlinelibrary.wiley.com/doi/10.1002/ajb2.16449 by Readcube (Labtiva Inc.), Wiley Online Library on [31/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License the detection of putative SLFs and the related functional inference (Xu et al., 2013; Williams et al., 2014,2015; Zhao et al., 2022a). SLFs found across eudicots apparently do not form a single clade. However, within smaller taxonomic groups, SLFs form closely related, although not necessarily monophyletic, clades (Figure S2). We recover a group of 18 pollen‐expressed F‐box domain‐containing genes in L. monelli, possibly involved in RSI response (see Figures S2 and S9 for gene genealogy and expression patterns). The group belongs to a larger clade that includes SLFs from other plant families with RSI, as well as anonymous SLFs, with uncharacterized functions. Among these 18 (Figure S2, red tips), 12 belong to a L. monelli‐ specific clade. Eleven of these are most likely to be SLFs. It is difficult to speculate on the identity and function of these recovered transcripts, in the absence of additional data, but they minimally comprise SLF candidate genes. One sequence in a neighboring clade is related to sequences that are non‐ specifically expressed (in tissues other than pollen). A distantly related clade harbors six additional sequences. Few complete S‐loci have been assembled to date, partly because they are regions of suppressed recombination, composed of closely related SLFs, are very large (often exceeding 10 Mb; Figure S5), and are replete with repetitive elements and transposons (Kubo et al., 2015;Wuetal.,2020). S‐RNase diversity and signature of diversifying selection Functional S‐loci are extraordinarily polymorphic because of strong negative frequency‐dependent selection, favoring rare S‐alleles, which can fertilize mates more reliably than common S‐alleles segregating in populations (Wright, 1939; Clark and Kao, 1991). We recover likely S‐RNases in L. monelli, with all associated hallmarks of this manner of balancing selection. They show high polymorphism and sequence divergence—often exceeding one substitution per site—as well as a high ratio of non‐synonymous to synonymous mutations. In total, we found 38 S‐alleles in 27 plants (unique genets), which is the minimum estimate for allele number in L. monelli. We also used the known sampling effort and allele recovery to generate ML estimates of allele numbers (Paxman, 1963). In 19 wild‐sampled individuals, each a distinct genet, we sequenced 34 unique S‐alleles (Figure 2). Excluding the greenhouse population entirely, 19 wild‐ collected individuals with 34 alleles yield N ˆ= 159 alleles (95% CI = [77,486]; Figure S3). In total, we recovered 21 independent genotypes composed of 37 alleles, which give N ˆ=15 5 (95% CI = [81,408]; Figure S3). This is likely a slightly downward biased estimate of species‐wide S‐allele diversity for this species, given that we discarded one allele ( S 3 4 ) from a paternal plant of uncertain complete genotype. These estimates are imprecise but clear in one respect: L. monelli harbors extremely high diversity at its S‐locus, perhaps higher than any locus recovered in plants so far, and very likely reflects large and stable long‐term effective population sizes (Richman et al., 1996;Lawrence,2000). Pairwise amino acid identities for the 38 S‐RNases were as low as 45% and upward of 93%, comparable to S‐RNases FIGURE 4 Gene tree with a sample of characterized S‐RNases sequenced to date. Lysimachia monelli sequences uncovered in the present study are shown alongside previously studied S‐RNases from six families with RSI. Pairwise divergences are exceedingly high among S‐RNases because rare‐allele advantage preserves polymorphism and favors non‐synonymous changes at many sites (see the scale at bottom‐left). This is in addition to the high expected divergence due to the time elapsed since the common ancestor of eudicots, which jointly complicates accurate inference of branching order. A larger gene tree displaying amuchgreatersampleofplantT2/S‐RNasegenefamilymembers,alongwith posterior branch support, is available in Figure S1. RAPID DETECTION OF RNASE‐BASED SELF‐INCOMPATIBILITY IN LYSIMACHIA MONELLI (PRIMULACEAE) | 9of16 15372197, 2025, 1, Downloaded from https://bsapubs.onlinelibrary.wiley.com/doi/10.1002/ajb2.16449 by Readcube (Labtiva Inc.), Wiley Online Library on [31/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Figure S7. Frequency distribution of isoelectric points (pI) for T2/S‐RNases in land plants. Figure S8. Crossing relationships among individuals with RNA‐seq identified genotypes. Figure S9. Expression of the F‐box domain‐containing gene family sequences recovered from Lysimachia monelli. Table S1. RNA‐seq and assembly summary statistics. Table S2. Assembled Lysimachia monelli T2/S‐RNase transcript information. Table S3. Detailed information about RNA‐seq sequence read filtering (removal). How to cite this article: Ramanauskas K, Jiménez‐López FJ, Sánchez‐Cabrera M, et al. 2025. Rapid detection of RNase‐based self‐incompatibility in Lysimachia monelli (Primulaceae). American Journal of Botany 112(1): e16449. https://doi.org/10.1002/ajb2.16449 16 of 16 | RAPID DETECTION OF RNASE‐BASED SELF‐INCOMPATIBILITY IN LYSIMACHIA MONELLI (PRIMULACEAE) 15372197, 2025, 1, Downloaded from https://bsapubs.onlinelibrary.wiley.com/doi/10.1002/ajb2.16449 by Readcube (Labtiva Inc.), Wiley Online Library on [31/01/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License