Extensive Introgression Failed to Erode Species Boundaries Among Multiple Sympatric Closely Related Species of Roscoea

Extensive Introgression Failed to Erode Species Boundaries Among Multiple Sympatric Closely Related Species of Roscoea | Authorea try { document.documentElement.classList.add('js'); } catch (e) { } var _gaq = _gaq || []; _gaq.push(['_setAccount', 'G-8VDV14Y67G']); _gaq.push(['_trackPageview']); (function() { var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true; ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js'; var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s); })(); Skip to main content Preprints Collections Wiley Open Research IET Open Research Ecological Society of Japan All Collections About About Authorea FAQs Contact Us Quick Search anywhere Search for preprint articles, keywords, etc. Search Search ADVANCED SEARCH SCROLL Molecular Ecology This is a preprint and has not been peer reviewed. Data may be preliminary. 26 March 2025 V1 Latest version Share on Extensive Introgression Failed to Erode Species Boundaries Among Multiple Sympatric Closely Related Species of Roscoea Authors : Hong-Fan Chen , Ryan A. Folk , Ya-Li Wang , Wen-Jing Wang , Guo-Jun Shao , Mei-Yuan Huang , Xiangqin Yu 0000-0001-5245-343X , Li Li , and Jian-Li Zhao 0000-0002-5137-7735 [email protected] Authors Info & Affiliations https://doi.org/10.22541/au.174299972.22109921/v1 Published Molecular Ecology Version of record Peer review timeline 421 views 250 downloads Contents Abstract Supplementary Material Information & Authors Metrics & Citations View Options References Figures Tables Media Share Abstract not-yet-known not-yet-known not-yet-known unknown How species boundaries are maintained among sympatric closely related species experiencing gene flow is a puzzling question in evolutionary biology. Although introgression is commonly documented, the dynamics and gene function of introgression have rarely been explored to probe why frequent introgression does not necessarily destroy species boundaries in sympatry. In this study, we employ whole-genome resequencing data to examine introgression among five closely related species of Roscoea that co-exist in a “sky island” with seventeen distinct morphological traits. Our findings reveal that introgression has led to the phylogenomic discordance between nuclear and chloroplast genomes among these morphologically distinct species. Additionally, introgression is predominantly asymmetrical in both intensity and gene function, particularly concerning recombination. Notably, the majority of gene functions associated with introgressed loci are unrelated to reproductive processes. Our results suggest that closely related species with incomplete allele assortment can coexist despite theoretical predictions, highlighting the semipermeable nature of species boundaries as reproductive isolation develops. This provides a critical conceptual framework for understanding the interplay between introgression and species persistence. Our finding offers insights into how related sympatric species boundaries can be maintained in the face of frequent asymmetrical gene introgression. Title: Extensive Introgression Failed to Erode Species Boundaries Among Multiple Sympatric Closely Related Species of Roscoea Runing title: Introgression and Species Boundaries Hong-Fan Chen 1,2 , Ryan Andrew Folk 3 , Ya-Li Wang 1 , Wen-Jing Wang 1 , Guo-Jun Shao 1 , Mei-Yuan Huang 1 , Xiang-Qin Yu 4 , Li Li 1 , Jian-Li Zhao *,1 1 State Key Laboratory for Vegetation Structure, Functions and Construction, Ministry of Education Key Laboratory for Transboundary Ecosecurity of Southwest China, and Yunnan Key Laboratory of Plant Reproductive Adaptation and Evolutionary Ecology, Institute of Biodiversity, School of Ecology and Environmental Science, Yunnan University, 650500, Kunming, China; 2 Institute of International Rivers and Eco-security, Yunnan University, Kunming, 650504, Yunnan, China; 3 Department of Biological Sciences, Mississippi State University, Mississippi State, 39762, Mississippi, USA; 4 CAS Key Laboratory for Plant Diversity and Biogeography in East Asia, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming, 650201, Yunnan, China. Correspondence: Jian-Li Zhao ( [email protected] ) Telephone/fax numbers: +8613320519502/+86087165031446 not-yet-known not-yet-known not-yet-known unknown Email and ORCID of Authors Hong-Fan Chen ( [email protected] ): https://orcid.org/0000-0001-7940-1875 Ryan Andrew Folk ( [email protected] ): https://orcid.org/0000-0002-5333-9273 Ya-Li Wang ( [email protected] ): https://orcid.org/0000-0002-1762-2805 Wen-Jing Wang ( [email protected] ): https://orcid.org/0000-0001-6886-6028 Guo-Jun Shao ( [email protected] ): https://orcid.org 0000-0002-3646-5526 Mei-Yuan Huang ( [email protected] ): https://orcid.org/0001-7820-5733 Xiang-Qin Yu ( [email protected] ): https://orcid.org/0000-0001-5245-343X Li Li ( [email protected] ): https://orcid.org/0000-0001-9657-2107 Jian-Li Zhao ( [email protected] ): https://orcid.org/0000-0002-5137-7735 Data Availability Statement The whole genome resequencing raw reads of this study are openly available in NCBI (BioProject: PRJNA904212) at https://www.ncbi.nlm.nih.gov/sra/PRJNA904212. Reference genome of Roscoea schneideriana can be available from corresponding author. Funding Statement This study was supported by the National Natural Science Foundation of China (No. 42171057 and No. 41871047); the Yunnan Fundamental Research Projects (No. 202301BF07000-001); the Yunnan Revitalization Talent Support Program “Young Talent Project” (No. YNWR-QNBJ-2019-214); the Project for Talent and Platform of Science and Technology in Yunnan Province Science and Technology Department (No. 202205AM070005); the Graduate Research Innovation Fund of Yunnan University (No. 2020237). Conflict of Interest The authors declare no conflict of interest. not-yet-known not-yet-known not-yet-known unknown Abstract How species boundaries are maintained among sympatric closely related species experiencing gene flow is a puzzling question in evolutionary biology. Although introgression is commonly documented, the dynamics and gene function of introgression have rarely been explored to probe why frequent introgression does not necessarily destroy species boundaries in sympatry. In this study, we employ whole-genome resequencing data to examine introgression among five closely related species of Roscoea that co-exist in a “sky island” with seventeen distinct morphological traits. Our findings reveal that introgression has led to the phylogenomic discordance between nuclear and chloroplast genomes among these morphologically distinct species. Additionally, introgression is predominantly asymmetrical in both intensity and gene function, particularly concerning recombination. Notably, the majority of gene functions associated with introgressed loci are unrelated to reproductive processes. Our results suggest that closely related species with incomplete allele assortment can coexist despite theoretical predictions, highlighting the semipermeable nature of species boundaries as reproductive isolation develops. This provides a critical conceptual framework for understanding the interplay between introgression and species persistence. Our finding offers insights into how related sympatric species boundaries can be maintained in the face of frequent asymmetrical gene introgression. Key words: asymmetrical introgression, closely related species, introgressed loci, species boundary, Roscoea . Introduction Introgression (or introgressive hybridization) is a prevalent and significant biological process in nature (Anderson 1953; Rieseberg and Wendel 1993; Abbott et al. 2013; Suarez-Gonzalez et al. 2018), traditionally defined as the incorporation of genes from one lineage (species or population) into the genome of another through hybridization and subsequent backcrossing (Anderson and Hubricht 1938; Harrison and Larson 2014; Jamie and Meier 2020). It is widely acknowledged that introgression substantially contributes to diversification by providing a rapid pathway to genetic variation, thereby facilitating adaptive radiation (Tigano and Friesen 2016; Suarez-Gonzalez et al. 2018; Edelman et al. 2019; Marques et al. 2019; Oziolor et al. 2019; Jamie and Meier 2020). However, simulations have predicted that species boundaries could deteriorate due to the homogenization of reproductive isolation alleles following introgression (Owens and Samuk 2020). This process could lead to the rapid fusion of one of the parent species, phenomenon known as ”speciation reversal”, for which empirical evidence exists (Seehausen 2006; Seehausen et al. 2008; Abbott et al. 2013; Yu et al. 2023). Although extinction resulting from introgression may be under-detected due to its potential involvement with cryptic diversity, it may be widespread across the Earth (Rhymer and Simberloff 1996). The risk of extinction is likely highest among closely related species with incomplete reproductive isolation that encounter secondary contact, as this scenario predicts low hybrid fitness (Irwin and Schluter 2022). Nonetheless, empirical evidence indicates that frequent introgression over time does not invariably lead to the weakening of species boundaries, such as the classic example of introgression in oak species (Khodwekar and Gailing 2017; Fu et al. 2022), leading to the unresolved question of how related species maintain boundaries under continuous introgression. Sympatric distribution is a central focus in introgression studies as it offers a common habitat for secondary contact and the survival of progeny. If introgressed loci are linked to key traits of the donor species, it is anticipated that the morphology of the recipient species will increasingly resemble that of the donor. For example, Rifkin et al. (2019) observed that the floral and inflorescence traits of Ipomoea cordatotriloba were significantly more similar to those of I. lacunosa sympatric regions, likely due to asymmetric introgression from I. lacunosa to I. cordatotriloba . Although sympatric environments evolve over time, a donor species may face extinction due to maladaptation, allowing introgressed descendants to persist with partial genomic elements of the extinct species, potentially exemplifying speciation reversal (Vonlanthen et al. 2012; Restani et al. 2018; Frei et al. 202 2). Nevertheless, species merging as a result of introgression is uncommon; instead, sympatric occurrences of closely related species with introgression are prevalent, a scenario likely sustained by divergent or disruptive selection. For instance, Larson et al. (2013) demonstrated that the maintenance of species boundaries between Gryllus firmus and G. pennsylvanicus, despite gene flow within a hybrid zone, was likely attributed to strong selection pressures and the limited extent of introgressed genes associated with pre-zygotic barriers. Similarly, following secondary contact between the invasive Helicoverpa armigera and the native H. zea in Brazil, species boundaries were robustly preserved despite extensive gene flow. This was evidence by H. zea acquiring an insecticide-resistant locus from H. armigera, while selection acted against most other introgressed variations (Valencia-Montoya et al. 2020). Additional studies have indicated that highly divergent genomic regions with low recombination rate, particularly those involved in reproduction and reproductive isolation, are instrumental in driving species persistence amidst hybridization (Malinsky et al. 2015; Gao et al. 2022). Furthermore, reinforcement can counteract gene flow, facilitating the completion of the speciation process and playing a crucial role in maintaining species boundaries between sympatric species (Marshall et al. 2002; Matute 2010; Hopkins 2013; Calabrese and Pfennig 2020). Although researches on introgression between related species in hybrid zones have advanced our comprehension of the mechanisms underlying species boundary maintenance, the comparative analysis of functions between introgressed and non-introgressed genes remains largely unexplored. This gap limits our understanding of the gene-level mechanisms that may contribute to species persistence in sympatric conditions. In the context of adaptive introgression, introgressed genes are expected be linked to specific adaptations to both biotic and abiotic environments. In contrast, no-introgressed genes may be subject to divergent selection or associated with hybrid maladaptation (Turner et al. 2005; Minder and Widmer 2008; Shaw and Mullen 2011; Harrison and Larson 2014). Therefore, a comparative analysis of the functions of introgressed and non-introgressed genes in sympatric, closely related species is essential. Such an analysis would enhance our understanding of why species boundaries are maintained despite ongoing introgression and aid in predicting whether species are face extinction threats due to introgression. Sympatric distributions of closely related species are frequently observed in mountainous regions, where encompass the majority global biodiversity hotspots (Myers et al. 2000; Mittermeier et al. 2011; Rahbek et al. 2019). This underscores the potentials of mountainous regions as prime locations for investigating species boundaries characterized by extensive introgressions. However, current studies in biodiversity hotspots have predominantly concentrated on speciation and species migration (Hughes and Atchison 2015; Chen et al. 2019; Ding et al. 2020), rather than on species boundaries maintenance in the species complexes prevalent in these areas. “Sky islands”, which are defined by their characteristically discontinuous distributions across isolated regions of suitable montane habitat, offer relatively ecologically and geologically insular environments conductive to species evolution with restricted gene flow (He and Jiang 2014; Hughes and Atchison 2015; Luo et al. 2016; Zhang et al. 2019; Yu et al. 2023). Liang-Wang Mountain, situated at approximately 2,800 m.a.s.l in the southwestern mountains of China and the Hengduan Mountain areas, represents the highest peak in central Yunnan Province. This mountain is encircled by three lakes, thereby forming a typical “sky island” (Figure 1A-C). Roscoea, an alpine endemic flowering plant of the ginger family Zingiberaceae, is widely distributed from the southern Himalayas to the Hengduan Mountains in southwestern China (Cowley et al. 2 007). Five closely related species in this genus, R. tibetica, R. debilis, R. lingbaoshanensis, R. kunmingensis, and R. praecox (Figure 1D), coexist above 2,500 m.a.s.l. on the top of Liang-Wang Mountain. These species exhibit overlapping flowering phenology and share similar habitat requirements. Notably, these five species formed a monophyletic group, with the exception of R. forrestii (Zhao et al. 2017. R. sp. in this study is R. lingbaoshanensis (Li et al. 2021)). Roscoea tibetica and R. lingbaoshanensis are widespread in the Hengduan Mountains and have several sympatric distribution sites (Li et al. 2021). Roscoea debilis is also widespread, albeit with few populations in the Hengduan Mountains (Cowley et al. 2 007). Roscoea praecox is found in the eastern part of the Hengduan Mountains with limited populations (Cowley et al. 2 007). Roscoea kunmingensis is represented by a single population on Liang-Wang Mountain, consistent with previous record (Cowley et al. 2 007). Consequently, Liang-Wang Mountain is the sole location where these five Roscoea species co-occur. Previous studies have suggested that interspecific introgression among closely related Roscoea species may frequently occur (Du et al 2012; Zhao et al. 2016; Zhao et al. 2021). Preliminary field investigations identified several individuals exhibiting suspected hybrid characteristics, displaying an intermediate phenotype between R. tibetica and either R. kunmingensis or R. lingbaoshanensis . This suggests that the area may represent a hybrid zone involving multiple species. Thus, the sympatric occurrence of these five Roscoea species within this “sky island” provides an exemplary system for in-depth exploration of the mechanisms maintaining species boundaries in a sympatric location. Given the distinct morphological differences among the five species, we hypothesize that potential introgressions have not disrupted species boundaries, positing that introgressed loci are not linked to reproductive traits or processes. To evaluate this hypothesis, we initially conducted observations of phenology and assessed morphological traits across five species. Subsequently, we reconstructed phylogenetic relationships utilizing whole genome resequencing data to verify the evolutionary independence of each species. Given that introgression is likely a causative factor in cytonuclear conflict between species, we proceeded to estimate the frequency and direction of introgression at both temporal and species-specific scales. We then compared the population recombination rates among the five species, and inferred potential gene functions for both introgressed and non-introgressed loci. This study ultimately seeks to elucidate the genetic framework of co-occurring individuals within a ”sky island” environment, aiming to determine the magnitude, direction, and functional implications of introgression among multiple sympatric, closely related species. The findings are intended to inform a conceptual framework for understanding semipermeable species boundaries, thereby addressing the interplay between sympatry, introgression, and the long-term evolutionary trajectories of species. 2 Materials and Methods 2.1 Phenological Observation and Morphological Traits Measurement The flowering duration for each species was documented, and seventeen morphological characteristics of the floral components were collected and measured. These characteristics included plant height, flower length (from the apex of the dorsal petal to the tip of the labellum), flower width (from the right to the left labellum), sepal length, dorsal petal length, lateral petal length and width, staminode length and width, labellum length, width, and depth, the angle between the style and stamen appendage, stigma diameter, length of fine hair on the stigma, and pollen diameter. Among these traits, the height and size of the flowers were measured in situ. The fresh floral parts were preserved in a mixed solution of ethanol and glycerol (ratio 9:1), subsequently dissected, and measured in the laboratory. Each trait was photographed using a fluorescent microscope (Leica DM3000, under fluorescent light), and the images were processed using ImageJ-v.1.53 (Schneider et al. 2012). To quantify the morphological differences between species, the principal component analysis (PCA) was performed with ggplot2 -v.3.3.5 (Wickham 2016) and the significantly different levels ( P < 0.05) of each trait among species were estimated by one-way ANOVA in SPSS statistics-v.22.0. not-yet-known not-yet-known not-yet-known unknown 2.2 Sampling, Resequencing and Variant Calling Thirty-seven individuals representing five species of Roscoea were collected from Liang-Wang Mountain (Figure 1B, C), with one individual of R. alpina from the Himalayas serving as the outgroup. The sample sizes for each species are detailed in Table 1. Fresh leaves from all samples were collected in the field and preserved in liquid nitrogen, subsequently stored at -80°C in the laboratory. Genomic DNA was extracted from the fresh leaves using the Plant Genomic DNA Kit (TIANGEN BIOTECH CO., LTD, Beijing, China), following the manufacturer’s protocols. To remove RNA, each DNA sample was treated with 10 μl of RNase (TIANGEN BIOTECH CO., LTD., Beijing, China). The DNA library was then prepared and sequenced on the Illumina NovaSeq 6000 platform, producing 150-bp paired-end reads with a sequencing depth of 30-fold for each individual. For quality control, the raw sequences were processed using Fastp-v0.21.0 (Chen et al., 2018) for filtering, and adapters were removed using Trimmomatic-v0.39 (Bolger et al., 2014). The quality of the paired reads was assessed with FastQC-v0.11.9 (https://github.com/s-andrews/FastQC). Subsequently, the reads were mapped to the chromosome-level whole genome of Roscoea schneideriana, which served as the reference genome (Chromosome number = 12, genome size = 1.35 Gb, contig N50 = 3.51 Mb, genome completeness = 98.20%, gene completeness = 94.80%; unpublished genome assembled and annotated by our team), using BWA-MEM (Li, 2013). All BAM files were sorted and de-duplicated using SAMtools-v1.11 (Li et al., 2009) and Picard (Broad Institute, 2019). The mapping rate for each individual was calculated using SAMtools flagstat and SAMtools stats. Although the mapping rate for the outgroup was lower than that of the ingroups (~85% vs. ~99%; Table S1), the potential bias due to missing data can be disregarded, as the SNP missing rate was low (Figure S1). Additionally, the outgroup was not included in the reconstruction of the species tree or in the estimation of PCA, genetic diversity, genetic structure, and recombination. Raw single nucleotide polymorphisms (SNPs) were extracted utilizing the GATK4 pipeline, employing the following parameters: ”QD < 2.0 || MQ 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0” (https://github.com/broadinstitute/gatk). To assess the missing rate, minor allele frequency (MAF), and the number of loci deviating from Hardy-Weinberg equilibrium (HWE), we employed PLINK-v1.90 (Chang et al., 2015). Visualization of these statistics was facilitated by the R scripts hist_miss.R, MAF_check.R, and hwe.R (https://github.com/MareesAT/GWA_tutorial). The linkage disequilibrium decay rate of the raw SNPs was calculated and visualized using PopLDdecay-v3.43 (Zhang et al., 2018). Based on the quality assessments of the raw SNPs (Figure S1), a high-quality whole-genome SNP dataset was generated using PLINK with the parameters ”–maf 0.14 –geno 0.05 –hwe 0.0001 –indep-pairwise 100 25 0.1”. For each individual, the chloroplast genome dataset was assembled using NOVOWrap-v1.20 and GetOrganelle-v1.7.7.0 (Jin et al., 2019; Wu et al., 2021), with the complete genome of R. tibetica chloroplast (accession number: MT180128.1) from NCBI serving as the chloroplast reference. Subsequently, multiple sequence alignments were performed using MAFFT version 7.490 (https://mafft.cbrc.jp/alignment/software). 2.3 Genetic Structure and Phylogenetic Reconstruction The phylogenetic tree was constructed utilizing the maximum likelihood approach implemented in IQ-TREE-v2.0.6 (Minh et al., 2020), applied to both the whole-genome SNPs dataset of nuclear DNA and the chloroplast genome dataset. The parameter ”-m MFP+MERGE” was employed to automatically determine the optimal nucleotide substitution model. The consensus tree was generated through a search involving 1,000 bootstrap replicates. Visualization of all phylogenetic trees was conducted using Figtree-v1.4.4 (https://github.com/rambaut/figtree). The consensus sequences for all individuals were extracted using ANGSD-v0.934 (Korneliussen et al., 2014) from the de-duplicated and sorted BAM files. Reads with a mapping quality (-minMapQ) and base quality score (-minQ) both below 20 were excluded from the analysis. Sequences were divided into 100-Kb windows, and the proportion of ’N’ bases within each window was calculated using Seqkit-v2.0.0 (Shen et al., 2016). Windows with an ’N’ proportion of 50% or greater were discarded. This process resulted in a dataset encompassing all chromosomes, consisting of 2,330 100-Kb consensus sequence windows. Subsequently, a gene tree was constructed for each window using IQ-TREE with consistent parameters, and a species tree was generated using multiple gene trees as input files with ASTRAL-v5.7.1 (Mirarab et al., 2014). Genetic structure was assessed using Principal Component Analysis (PCA) in PLINK and ADMIXTURE-v1.3.0 (Alexander et al., 2009) for all individuals across five species, with the range of K from 1 to 10 and 1,000 bootstrap replicates. The online tool Pophelper-v1.0.10 (http://pophelper.com) was subsequently employed to visualize the genetic component columns. 2.4 Test of Introgression To enhance the efficiency of our analysis, we selected five homozygous individuals—Rde-8 of R. debilis , Rku-23 of R. kunmingensis , Rli-10 of R. lingbaoshanensis , Rti-33 of R. tibetica , and Rpr-4 of R. praecox —based on their mapping rates (Table S1) and the homozygosity of their genetic components. These individuals, exhibiting the highest mapping rates, were chosen to represent their respective species in the introgression analysis. Additionally, R. alpina , a distantly related species from the Himalayas, was included as the outgroup. We investigated the roles of introgression and incomplete lineage sorting (ILS) in Roscoea by analyzing 2,330 100-Kb consensus windows. The five species were organized into ten distinct species triplets, each comprising different combinations (Table S2). For each triplet, gene trees were constructed using IQ-TREE. These gene trees then served as input topologies for QuIBL tests, as described by Edelman et al. (2019). QuIBL provided estimates of the proportion of the genome affected by introgression or ILS and assigned gene origins as introgressed or non-introgressed. Model selection between the Non-ILS and ILS models was based on the Bayesian information criterion (BIC), with a lower BIC value and a significant difference ( P < 0.01) indicating the preferred model. To mitigate the potential influence of window size, we also generated a dataset comprising 3,263 20-Kb consensus sequence windows, applied the same analytical pipeline, and conducted QuIBL tests to compare results across different window sizes. Both datasets yielded consistent results, as detailed in the Results section. We employed 2,330 100-Kb consensus sequences for the purpose of introgression analysis. To ascertain the direction of introgression, we utilized D foil , which is based on the f 5 -ratio test as described by Pease and Hahn (2015). D foil is a site-based statistical method designed to identify both the presence and directionality of hybridization using five-taxon trees. This method leverages the fact that ILS predicts specific patterns of discordant topologies relative to the ”true” topology, while introgression introduces imbalances in these statistical expectations. The five individuals, representing different species, were categorized into 60 groups after arrangement and combination (Table S3), encompassing all possible phylogenetic relationships. Subsequently, we employed fasta2dfoil.py to generate 139,800 count files and dfoil.py to identify introgression events based on all count files. Based on the presence of introgressive signals in the D foil outputs, the 2,330 consensus sequence windows were categorized into introgressed and non-introgressed loci. Subsequently, we calculated the frequency of each introgression direction among species. To estimate divergence times for the 2,330 windows, we employed BEAST-v2.6.7 (Bouckaert et al., 2019). The analyses utilized the substitution model selected by ModelTest-NG (Darriba et al., 2020) for each window, applying a strict clock rate of 3.18 × 10⁻⁹, which corresponds to the divergence between Roscoea and Musa (unpublished data). The Markov Chain Monte Carlo (MCMC) analysis was conducted with 50 million chains, with outputs recorded every 5,000 steps. For each window, the initial 10% of the trees were discarded as burn-in samples, and a maximum clade credibility (MCC) tree was generated using TreeAnnotator-v2.6.7 (Drummond and Rambaut, 2007). After analyzing the branch times of 2,330 MCC trees, we estimated the minimum timing of introgression events by examining the divergence times of each introgressed window between introgressive sister species. Theoretically, the earliest introgression events between sister species should occur after their divergence time. Since loci require time to exhibit variation post-introgression, we excluded introgression events involving ancestral clades from our analysis. Ultimately, we utilized the R package pheatmap (Kolde, 2019) to summarize the frequency and direction of introgression events per million years between species. 2.5 Estimation of Recombination Rates and Nucleotide Diversity Of the 37 individuals analyzed, putative hybrids were excluded, and genotype imputation for SNP sites was conducted using Beagle-v5.3 (Browning et al., 2018). We estimated the population recombination rates ( ρ ) for five species across 12 chromosomes utilizing FastEPRR-v2.0 (Hao et al., 2022). Subsequently, we extracted the corresponding SNPs within the ranges of 2,330 consensus sequence windows and estimated the recombination rates. We compared the recombination rates between introgressed and non-introgressed loci at both the species and chromosome levels. Additionally, nucleotide diversity ( π ) for the five species was calculated using VCFtools-v0.1.13 (Danecek et al., 2011), with comparisons made across different levels and windows. For this analysis, we employed SNPs without missing data and minor allele sites, applying parameter settings (–maf 0.05 –geno 0.2), and utilized a 100-kb window size for calculating π . 2.6 Gene Alignment and GO-enrichment We utilized Bedtools-v2.30.0 (Quinlan and Hall 2010) to extract complete sequences from the reference genome based on the structural annotation file. The module makeblastdb from BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) was employed to construct the subject sequence databases. Then, we executed blastn-v2.12.0 using 2,330 consensus sequence windows as query sequences, specifying the parameters ”-max_target_seqs 1 -evalue 1e-20 -outfmt 6” to obtain the optimal alignment results. In instances where a window aligned with multiple fragments, we retained only the fragment with the highest identity rate and discarded BLAST results with identity rates below 95%. As a result, we identified the genes present in the 2,330 consensus sequence windows following alignment. Gene ontology (GO) enrichment analysis was performed on multiple alignment results. Initially, EggNOG-mapper-v1.0.3 (Huerta-Cepas et al., 2017) was employed to annotate protein sequences, facilitating the construction of the non-model species OrgDB using the R packages AnnotationForge -v1.36.0 (Carlson and Pagès, 2023), dplyr -v1.0.8 (Wickham et al., 2023), and stringr -v1.4.0 (Wickham, 2023). ClusterProfiler -v4.2.2 (Yu et al., 2012) was utilized to conduct GO enrichment analysis with aligned gene IDs from introgressed and non-introgressed loci in each direction of introgression. A significance threshold of P < 0.05 was applied to determine significant GO enrichment. To ensure the GO enrichment results were not confounded by other species, the GO-IDs specific to the two species involved in introgression were extracted. Dotplot -v4.3.0 (https://www.rdocumentation.org/packages/Seurat/versions/4.3.0/topics/DotPlot) was employed to generate visualizations of the top 10 GO functions. 2.7 Reconstruction of Effective Population Sizes We estimated the effective population sizes for five species across historical periods utilizing the Pairwise Sequential Markovian Coalescent (PSMC) model (available at https://github.com/lh3/psmc). The parameters employed were “-N25 -t15 -r5 -b -p 4+25*2+4+6”. A nucleotide substitution rate of 3.18 × 10 -9 per site per year was applied, with a generation time of 2 years, as reported by Li et al. (2021). Additionally, five heterozygous individuals, identified through population structure analyses, were excluded from the study. not-yet-known not-yet-known not-yet-known unknown 3 Results not-yet-known not-yet-known not-yet-known unknown 3.1 Phenology and Morphology The flowering periods of R. kunmingensis and R. tibetica occur from late April to early June and from late May to the end of June, respectively. R. praecox and R. tibetica exhibit overlapping flowering periods with these two species. In contrast, R. debilis typically blooms in July, which is temporally distinct from the flowering periods of the other four species (Figure 1E). Morphometric PCA and multiple comparisons of seventeen traits demonstrate that the five Roscoea species are significantly distinct from one another (P < 0.05 for all traits; Figure 1F, S2, and Table S4), corroborating the morphological differentiation observed. The PCA plot positions R. kunmingensis between R. tibetica and R. praecox (Figure 1F), suggesting some degree of morphological intermediacy. The pollen size of R. tibetica is nearly identical to that of R. debilis, whereas R. kunmingensis exhibits the largest pollen size. Roscoea praecox occupies an intermediate position between R. kunmingensis and R. debilis in terms of pollen size (Figure S2). 3.2 Genetic and Phylogenetic Structure We extracted 65,922 whole-genome single nucleotide polymorphisms (SNPs) from an initial set of 53,930,487 raw SNP sites based on quality criteria (Figure S1). We also obtained a chloroplast genome consensus length of 166,959 bp. The monophyly of each species was robustly supported by the majority of individuals in the whole-genome SNP tree, species tree, and PCA (Figure 2A, C), with the exception of five individuals (Rti-47 of R. tibetica , Rku-18 of R. kunmingensis , and three putative hybrid individuals) that exhibited distinctive heterozygous components in their genetic structure (Figure 2D). Conversely, the chloroplast genome tree did not support the monophyly of any species, as numerous individuals from different species were intermixed (Figure 2B). Furthermore, the histogram of genetic components (Figure 2D) indicated that multiple individuals displayed mixed genetic components when the number of ancestral populations (K) was set between 3 and 5 in the ADMIXTURE analysis. These findings suggest the presence of hybridization, introgression, or incomplete lineage sorting (ILS) among the five species. 3.3 Interspecific Introgression The QuIBL analyses, conducted on 2,330 100-Kb windows and 3,263 20-Kb windows, revealed the similar results that the Non-ILS model exhibited a lower average Bayesian Information Criterion (BIC) value (average BIC = -7595.472 and -9653.914, respectively) compared to the ILS model (average BIC = -7043.193 and -9313.374, respectively), with statistically significant differences ( P < 0.01; see Table S5 and Table S6). Additionally, the estimated introgression ratio for all tested windows exceeded 99% and 94%, as determined by the mixed proportion (Table S5). These results provide robust evidence supporting the occurrence of interspecific introgression. D foil tests identified 5,106 instances of introgression signals among species pairs out of a total of 139,800 tests conducted. The frequency of introgression events was notably higher when R. praecox was involved, whereas it was lower for R. tibetica (Figure 3A, Table S7). Introgression from R. debilis , R. kunmingensis and R. lingbaoshanensis into R. praecox occurred more frequently than introgression in the reverse direction. Conversely, introgression between R. tibetica and the trio of R. debilis , R. kunmingensis and R. lingbaoshanensis was relatively infrequent and exhibited near parity in both directions (Figure 3B). The distribution of inferred introgressions was uneven, predominantly occurring between 1 and 4 million years ago (Mya), with a notable peak between 2 and 3 Mya, followed by a substantial decline since 1 Mya (Figure 3A). not-yet-known not-yet-known not-yet-known unknown 3.4 Recombination Rates and Nucleotide Diversity Recombination rates (ρ ) for each chromosome are presented in Figure S3. This study focuses on the differences in ρ between introgressed and non-introgressed loci. At the species level, recombination rates differ significantly among the five species examined (P < 0.01; Figure 3C, Table S8). Specifically, the ρ values for R. praecox and R. tibetica are substantially higher than those for R. debilis, R. kunmingensis and R. lingbaoshanensis . In both R. praecox and R. tibetica, the ρ values at introgressed loci are significantly greater than those at non-introgressed loci (P < 0.01; Figure 3C, Table S9). At the chromosomal level, no single chromosome consistently exhibits significant differences in ρ between introgressed and non-introgressed loci. However, in general, the ρ values at introgressed loci are significantly higher than those at non-introgressed loci (Table S9). Furthermore, there are significant differences in nucleotide diversity (π ) among the species (P < 0.01; Table S10), with R. praecox exhibiting the highest π, followed by R. kunmingensis, R. tibetica, R. lingbaoshanensis and R. debilis (Figure 3D). not-yet-known not-yet-known not-yet-known unknown 3.5 Functions of Introgressive Genes The functions of introgressed loci among most species were not associated with reproductive traits for the majority of species pairs. However, numerous genes related to reproductive processes were identified within the non-introgressed loci, demonstrating significant enrichment (P < 0.01; Table S11). These genes were involved in processes such as flower development (GO:0009908), regulation of flower development (GO:0009909), negative regulation of flower development (GO:0009910), and reproductive shoot system development (GO:0090567). These Gene Ontology (GO) terms collectively represent a comprehensive biological process, as documented in the QuickGO online database (https://www.ebi.ac.uk/QuickGO; Figure S4). An exception was identified in the study: although the gene functions of introgressed loci from R. kunmingensis into R. praecox did not show significant involvement in reproductive terms (Figure 4A), numerous reproductive terms were associated with the introgressed loci from R. praecox into R. kunmingensis . These terms include growth (GO:0040007), pollination (GO:0009856), pollen germination (GO:0009846), pollen tube development and growth (GO:0048868, GO:0009860), embryo sac development (GO:0009553), and negative regulation of post-embryonic development (GO:0048581) (Figure 4B). Collectively, these GO terms constitute a comprehensive biological process related to plant growth and reproduction, as indicated in the QuickGO online database (Figure 4C, S4). not-yet-known not-yet-known not-yet-known unknown 3.6 Changes of Effective Population Size The historical effective population sizes of five species indicate that certain individuals within R. tibetica and R. kunmingensis, as well as all individuals of R. praecox, demonstrated similar times to the most recent common ancestor (TMCRA; see Figure 5). Around 3 million years ago (Mya), the effective population sizes of R. tibetica and R. kunmingensis began to increase, reaching a peak approximately 1 Mya. In contrast, the effective population sizes of R. debilis, R. lingbaoshanensis, and some individuals of R. tibetica and R. kunmingensis exhibited a gradual increase from before 1 Mya to 0.7 Mya, followed by a pronounced rise around 3 thousand years ago (Kyr). The effective population size of R. praecox experienced a rapid increase from approximately 1 Mya to 0.4 Mya. With the exception of R. praecox, the effective population sizes of the other four species approached zero after 10 Kyr. 4 Discussion The distinctiveness among the five sympatric Roscoea species on the Liang-Wang Mountain “sky island” was robustly supported by all measured morphological traits and reciprocal monophyly in both whole-genome SNP tree and species tree based on nuclear DNA. However, the nuclear-DNA based trees showed conflict phylogenetic relationships from the chloroplast genome tree among species (Figure 2A, B), suggesting the presence of interspecific introgression or incomplete lineage soring (ILS). Despite the absence of sampling from allopatric populations, introgression tests among the five closely related sympatric species offer valuable insights into the mechanisms maintaining species boundaries in sympatric environments. Population PCA clustering (Figure 2C) and genetic structure components (Figure 2D) indicate extensive interspecific genetic ’mixing’ among the five species. The QuIBL tests suggest that this phylogenetic discordance is attributable to introgression. However, this introgression is not temporally uniform, as the majority of the detected introgression signals date back to periods before 1 million years ago (Mya). Further analysis revealed that asymmetric introgression in terms of intensity, direction, and gene function plays a crucial role in maintaining species boundaries among closely related species. In this study, we examine how closely related Roscoea species persist in sympatric locations, influenced by asymmetrical introgression. not-yet-known not-yet-known not-yet-known unknown 4.1 Asymmetrical Introgression with Intensity and Gene Functions “Species boundaries are semipermeable, with permeability (gene exchange) being a function of genome region” (Harrison and Larson 2014). This assertion implies that the introgression of functional genes and other genetic material between species is contingent upon specific genomic regions. Our findings indicate that the pattern of gene introgression between Roscoea species is asymmetrical, both in the intensity of gene flow and the involvement of functional genes, thereby confirming the semipermeable nature of species boundaries. In the context of semipermeable species boundaries, one proposed scenario posits that genes associated with reproductive isolation and those subject to divergent selection are retained within each species and exhibit resistance to introgression, whereas neutral genes are more readily exchanged between species (Barton 1979; Wu 2001; Harrison and Larson 2014). Genomic evidence further indicates that the rate of introgression is reduced for functionally significant genomic regions or genes (Calfee et al. 2020; Moran et al. 2021). This scenario is supported by the observation of bidirectional introgression of genes between Roscoea species, where the functions of introgressed genes are unrelated to reproduction. However, a notable exception was identified in the introgression from R. praecox into R. kunmingensis, which included genes with functions clearly related to reproduction. It is important to note that, as we rely on inferred functions based on gene annotations, we cannot conclusively determine whether all seemingly irrelevant introgressed genes might indirectly influence reproductive isolation. Additionally, it is plausible that some of these genes are subject to divergent selection. Nevertheless, the general absence of gene functions directly relevant to reproductive pathways strongly suggests that the maintenance of species boundaries in Roscoea was associated with the persistence of reproductive isolation and resistance of associated genes to introgression. Asymmetric introgressions can be influenced by variations in mating systems between species (Pickup et al. 2019). Mating systems affect the directionality of introgression through two primary mechanisms. Firstly, during hybridization events between self-fertilizing and outcrossing species, hybrid offspring produced by self-fertilizing species often exhibit reduced fitness, thereby impeding introgression from outcrossing species into self-fertilizing species (Pickup et al. 2019). Secondly, the fate of introgressed alleles differs depending on the direction of introgression. Alleles introgressed from self-fertilizing species into outcrossing species may be subject to higher recombination rates, which can facilitate the movement of neutral alleles and allow them to evade the purging of linked non-neutral loci by natural selection (Glémin et al. 2019). Conversely, the lower recombination rates characteristic of self-fertilizing species may hinder the transmission of introgressed alleles to subsequent generations due to maladaptive effects associated with linked non-neutral loci (Pickup et al. 2019). Our field investigation on Liang-Wang Mountain revealed that four species are capable of self-fertilization, with the exception of the outcrossing species R. praecox . Notably, R. lingbaoshanensis predominantly maintains its flowers in a closed state and exhibits a shorter flowering period in the wild, which likely facilitates self-pollination (Figure 1E). However, occasional flower openings may have allowed for gene flow, although the historical flowering patterns remain unknown. Both R. tibetica and R. kunmingensis exhibit characteristics of both outcrossing and selfing. Nonetheless, R. tibetica shows a preference for outcrossing, whereas R. kunmingensis tends towards selfing (Figure 3D), as evidenced by reduced molecular diversity and increased linkage disequilibrium in their genomes (Cutter 2019; Wang et al. 2021; Zhang et al. 2022). Additionally, the stigma of R. debilis is typically enclosed within the pollen sacs (Figure 1D), enabling delayed selfing through the secretion of stigmatic fluid in the absence of pollinators (Fan and Li 2012). Introgression tests of Roscoea indicated that the signal of introgressions from species with the ability of selfing to the outcrossing species (R. tibetica and R. praecox ) was stronger than introgression in the reverse direction (Figure 3A, B), which suggests that asymmetric introgressions were possibly the result of a higher capacity of crossing species to integrate genomic regions from selfing species through genome recombination. Theoretical models suggest that regions with high recombination rates should exhibit increased introgression, whereas regions with lower recombination rates should display reduced introgression (Barton and Bengtsson 1986; Martin and Jiggins 2017; Martin et al. 2019). Our findings corroborate this hypothesis, revealing that the recombination rates of introgressed loci were significantly higher than those of non-introgressed loci in R. tibetica and R. praecox across all chromosomes of both species (Figure 3C). Previous studies have demonstrated that genomic regions with elevated recombination rates can more rapidly disassemble heterologous fragments, eliminate deleterious alleles, and reduce linkage disequilibrium between alleles at different loci, thereby facilitating the generation of novel allele combinations that are subject to selection (Nachman and Payseur 2012; Schumer et al. 2018; Edelman et al. 2019; Feng et al. 2023). Consequently, for R. tibetica and R. praecox, maintaining higher recombination rates (Figure 3C) is advantageous for the accelerated breakdown of heterogenic genomic regions following introgression. This process promotes asymmetrical introgression, influencing both the intensity and functional implications of gene flow, thereby contributing to the maintenance of species boundaries. 4.2 Decrease in introgression over time An alternative scenario suggested by the semipermeable nature of species boundaries, which aligns with classic reinforcement theory, posits that closely related species possess a limited number of genes involved in reproductive isolation during the initial stages of speciation. This limited genetic involvement permits introgression across many other genomic regions. Over time, as the number of genes associated with reproductive isolation increases, and assuming these genes are unlinked, the proportion of the genome available for introgressive gene exchange is expected to diminish, thereby imposing constraints on the extent of introgression. Consequently, this scenario predicts a reduction in the permeability of species boundaries over time in response to introgression (Wu 2001; Harrison and Larson 2014). Conversely, genes exhibiting reduced introgression between species over time can be identified as potential reproductive isolation genes, contributing to the establishment of reproductive barriers (Harrison and Larson 2014). In the case of Roscoea species, the extent of introgression was not temporally uniform, displaying a marked decline following an initial peak. Such temporal changes suggest that reinforcement may account for the observed reduction in introgression over time. According to reinforcement theory, selection should act to prevent hybridization when divergent species come into contact, thereby avoiding the evolutionary costs associated with maladapted hybrids. Consequently, sympatric species exhibit stronger interspecific reproductive isolation mechanisms compared to allopatric species (Hopkins, 2013). Notably, selfing species can develop more robust reproductive isolation and reinforcement (Mallet, 2005; Pickup et al., 2019). As previously discussed, reinforcement may be augmented through selfing due to both maladaptation and recombination effects. This could further explain the observed substantial reduction in detected introgression after 1 Mya (Figure 3A), suggesting that semipermeability has decreased over time to reinforce species boundaries. not-yet-known not-yet-known not-yet-known unknown 4.3 Species Boundary Maintenance Introgression between morphologically distinctive species can be interpreted as a consequence of incomplete assortative mating. Theoretical models often predict that closely related species exhibiting incomplete assortative mating will be unable to coexist (Irwin and Schluter 2022). However, this prediction is contradicted by empirical evidence of sympatry and introgression in numerous plant taxa. In the genus Roscoea, introgression primarily manifests as asymmetrical among species, both in terms of intensity and the putative function of the genes involved, as well as temporal non-uniformity. Our findings align with a semipermeable model of species boundaries, which posits that certain genomic regions are more amenable to introgression between species, and that these regions vary over time. It is likely that genomic restrictions are mediated by genes associated with reproductive functions, which are expected to resist introgression. Furthermore, our results support the hypothesis that breeding systems play a significant role in shaping the asymmetry of introgression. Finally, we present clear evidence within the Roscoea system that reproductive boundaries are becoming more pronounced over time, and that this gradual accumulation of reproductive isolation is the primary mechanism facilitating the frequent occurrence of introgression and easily recognizable “good” species. In summary, this study elucidates semipermeable species boundaries as a crucial conceptual framework for addressing the conflict between sympatry, introgression, and the long-term fate of species. Acknowledgments This study was supported by the National Natural Science Foundation of China (No. 42171057 and No. 41871047); the Yunnan Fundamental Research Projects (Grant No. 202301BF07000-001); the Yunnan Revitalization Talent Support Program “Young Talent Project” (No. YNWR-QNBJ-2019-214); the Project for Talent and Platform of Science and Technology in Yunnan Province Science and Technology Department (No. 202205AM070005); the Graduate Research Innovation Fund of Yunnan University (No. 2020237). We thank Professor Qing-Jun Li even he passed away during revision of this manuscript, this paper is also in memory of him for his significant contributions to this work and plant evolutionary ecology. not-yet-known not-yet-known not-yet-known unknown not-yet-known not-yet-known not-yet-known unknown References Abbott RJ, Albach D, Ansell S, Arntzen JW, Baird SJ, Bierne N, Boughman J, Brelsford A, Buerkle CA, Buggs R et al. 2013. Hybridization and speciation. Journal of Evolutionary Biology 26:229-246. Alexander DH, Novembre J, Lange K. 2009. Fast model-based estimation of ancestry in unrelated individuals. Genome Research 19:1655-1664. Anderson E. 1953. Introgressive hybridization. Biological Reviews 28:280-307. Anderson E, Hubricht L. 1938. Hybridisation in Tradescantia . III The evidence for introgressive hybridisation. American Journal of Botany 25:396-402. Barton N, Bengtsson BO. 1986. The barrier to genetic exchange between hybridising populations. Heredity 56:357-376. Barton NH. 1979. Gene flow past a cline. Heredity 43:333-339. Bolger AM, Marc L, Bjoern U. 2014. Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics 30:2114-2120. Bouckaert R, Vaughan TG, Barido-Sottani J, Duchêne S, Fourment M, Gavryushkina A, Heled J, Jones G, Kühnert D, Maio ND et al. 2019. BEAST2.5: an advanced software platform for Bayesian evolutionary analysis. PLoS Computational Biology 15:e1006650. Broad-Institute. 2019. Picard Toolkit. Broad Institute, GitHub repository, https://broadinstitute.github.io/picard. Browning BL, Zhou Y, Browning SR. 2018. A one-penny imputed genome from Next-Generation reference panels. American Journal of Human Genetics 103:338-348. Calabrese GM, Pfennig KS. 2020. Reinforcement and the proliferation of species. Journal of Heredity 111:138-146. Calfee E, Agra MN, Palacio MA, Ramírez SR, Coop G. 2020. Selection and hybridization shaped the rapid spread of African honey bee ancestry in the Americas. PLoS Genetics 16:e1009038. Carlson, M., and H. Pagès. 2023. AnnotationForge: tools for building SQLite-based annotation data packages. Doi: 10.18129/B9.bioc.AnnotationForge. Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. 2015. Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience 4:7. Chen JH, Huang Y, Brachi B, Yun QZ, Zhang W, Lu W, Li HN, Li WQ, Sun XD, Wang GY et al. 2019. Genome-wide analysis of cushion willow provides insights into alpine plant divergence in a biodiversity hotspot. Nature Communications 10:5230. Chen SF, Zhou YQ, Chen YR, Gu J. 2018. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34:884-890. Cowley J, Wilford R, Bream R. 2007. The genus Roscoea . UK: The Royal Botanic Gardens, Kew. Cutter AD. 2019. Reproductive transitions in plants and animals: selfing syndrome, sexual selection and speciation. New Phytologist 224:1080-1094. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker R, Lunter G, Marth G, Sherry ST et al. 2011. The variant call format and VCFtools. Bioinformatics 27:2156-2158. Darriba D, Posada D, Kozlov AM, Stamatakis A, Morel B, Flouri T. 2020. ModelTest-NG: a new and scalable tool for the selection of DNA and protein evolutionary models. Molecular Biology and Evolution 37:291-294. Ding WN, Ree RH, Spicer RA, Xing YW. 2020. Ancient orogenic and monsoon-driven assembly of the world’s richest temperate alpine flora. Science 369:578-581. Drummond AJ, Rambaut A. 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Ecology and Evolution 7:214. Du GH, Zhang ZQ, Li QJ. 2012. Morphological and molecular evidence for natural hybridization in sympatric population of Roscoea humeana and R. cautleoides (Zingiberaceae). Journal of Plant Research 125: 595-603. Edelman NB, Frandsen PB, Miyagi M, Clavijo B, Davey J, Dikow RB, García-Accinelli G, Belleghem SMV, Patterson N, Neafsey DE et al. 2019. Genomic architecture and introgression shape a butterfly radiation. Science 366:594-599. Fan YL, Li QJ. 2012. Stigmatic fluid aids self-pollination in Roscoea debilis (Zingiberaceae): a new delayed selfing mechanism. Annals of Botany 110:969-975. Feng C, Wang J, Liston A, Kang M. 2023. Recombination variation shapes phylogeny and introgression in wild diploid strawberries. Molecular Biology and Evolution 40:msad049. Frei D, De-Kayne R, Selz OM, Seehausen O, Feulner PG. 2022. Genomic variation from an extinct species is retained in the extant radiation following speciation reversal. Nature Ecology and Evolution 6:461-468. Fu RR, Zhu YX, Liu Y, Feng Y, Lu RS, Li Y, Li P, Kremer A, Lascoux M, Chen J. 2022. Genome-wide analyses of introgression between two sympatric Asian oak species. Nature Ecology and Evolution 6:924-935. Gao W, Yu CX, Zhou WW, Zhang BL, Chambers EA, Dahn HA, Jin JQ, Murphy RW, Zhang YP, Che J. 2022. Species persistence with hybridization in toad-headed lizards driven by divergent selection and low recombination. Molecular Biology and Evolution 39:msac064. Glémin S, François CM, Galtier N. 2019. Genome evolution in outcrossing vs. selfing vs. asexual species. Methods in Molecular Biology 1910:331-369. Hao ZQ, Du PY, Pan YH, Li HP. 2022. Fine human genetic map based on UK10K data set. Human Genetics 141:273-281. Harrison RG, Larson EL. 2014. Hybridization, introgression, and the nature of species boundaries. Journal of Heredity 105:795-809. He K, Jiang XL. 2014. Sky islands of southwest China. I: an overview of phylogeographic patterns. Chinese Science Bulletin 59:1-13. Hopkins R. 2013. Reinforcement in plants. New Phytologist 197:1095-1103. Huerta-Cepas J, Forslund K, Coelho LP, Szklarczyk D, Jensen LJ, Mering Cv, Bork P. 2017. Fast genome-wide functional annotation through orthology assignment by eggNOG-mapper. Molecular Biology and Evolution 34:2115-2122. Hughes CE, Atchison GW. 2015. The ubiquity of alpine plant radiations: from the Andes to the Hengduan Mountains. New Phytologist 207:275-282. Irwin D, Schluter D. 2022. Hybridization and the coexistence of species. The American Naturalist 200:e93-e109. Jamie GA, Meier JI. 2020. The persistence of polymorphisms across species radiations. Trends in Ecology and Evolution 35:795-808. Jin JJ, Yu WB, Yang JB, Song Y, dePamphilis CW, Yi TS, Li DZ. 2019. GetOrganelle: a fast and versatile toolkit for accurate de novo assembly of organelle genomes. Genome Biology 21:241. Khodwekar S, Gailing O. 2017. Evidence for environment-dependent introgression of adaptive genes between two red oak species with different drought adaptations. American Journal of Botany 104:1088-1098. Kolde R. 2019. pheatmap: Pretty Heatmaps. Doi: 10.32614/CRAN.package.pheatmap. Korneliussen TS, Albrechtsen A, Nielsen R. 2014. ANGSD: analysis of next generation sequencing data. Bioinformatics 15:356. Larson EL, Becker CG, Bondra ER, Harrison RG. 2013. Structure of a mosaic hybrid zone between the field crickets Gryllus firmus and G. pennsylvanicus . Ecology and Evolution 3:985-1002. Li C, Wu YJ, Chen BC, Cai YD, Guo JZ, Leonard AS, Kalds P, Zhou SW, Zhang JC, Zhou P et al. 2022. Markhor-derived introgression of a genomic region encompassing PAPSS2 confers high-altitude adaptability in Tibetan goats. Molecular Biology and Evolution 39:msac253. Li H. 2013. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997v2 [q-bio.GN]. Li H, Handsaker B, Wysoke A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. 2009. The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics 25:2078-2079. Li L, Zhang J, Lu ZQ, Zhao JL, Li QJ. 2021. Genomic data reveal two distinct species from the widespread alpine ginger Roscoea tibetica Batalin (Zingiberaceae). Journal of Systematics and Evolution 59:1232-1243. Luo YH, Liu J, Tan SL, Cadotte MW, Xu K, Gao LM, Li DZ. 2016. Trait variation and functional diversity maintenance of understory herbaceous species coexisting along an elevational gradient in Yulong Mountain, Southwest China. Plant Diversity 38:303-311. Malinsky M, Challis RJ, Tyers AM, Schiffels S, Terai Y, Ngatunga BP, Miska EA, Durbin R, Genner MJ, Turner GF. 2015. Genomic islands of speciation separate cichlid ecomorphs in an East African crater lake. Science 350:1493-1498. Mallet J. 2005. Hybridization as an invasion of the genome. Trends in Ecology and Evolution 20:229-237. Marques DA, Meier JI, Seehausen O. 2019. A combinatorial view on speciation and adaptive radiation. Trends in Ecology and Evolution 34:531-544. Marshall JL, Arnold, ML, Howard DJ. 2002. Reinforcement: the road not taken. Trends in Ecology and Evolution 17:558-563. Martin SH, Davey JW, Salazar C, Jiggins CD. 2019. Recombination rate variation shapes barriers to introgression across butterfly genomes. PLoS Biology 17:e2006288. Martin SH, Jiggins CD. 2017. Interpreting the genomic landscape of introgression. Current Opinion in Genetics and Development 47:69-74. Matute DR. 2010. Reinforcement can overcome gene flow during speciation in Drosophila . Current Biology 20:2229-2233. Minder AM, Widmer A. 2008. A population genomic analysis of species boundaries: neutral processes, adaptive divergence and introgression between two hybridizing plant species. Molecular Ecology 17:1552-1563. Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams M, Haeseler Av, Lanfear R. 2020. IQ-TREE2: new models and efficient methods for phylogenetic inference in the genomic era. Molecular Biology and Evolution 37:1530-1534. Mirarab S, Reaz R, Bayzid MS, Zimmermann T, Swenson MS, Warnow T. 2014. ASTRAL: genome-scale coalescent-based species tree estimation. Bioinformatics 30:i541-i548. Mittermeier RA, Turner WR, Larsen FW, Brooks TM, Gascon C. 2011. Global Biodiversity Conservation: The Critical Role of Hotspots. In: Zachos FE, Habel JC, editors. Biodiversity Hotspots. Berlin, Heidelberg: Springer. Moran BM, Payne C, Langdon Q, Powell DL, Brandvain Y, Schumer M. 2021. The genomic consequences of hybridization. eLife 10:e69016. Myers N, Mittermeier RA, Mittermeier CG, Fonseca GAB, Kent J. 2000. Biodiversity hotspots for conservation priorities. Nature 403:853-858. Nachman MW, Payseur BA. 2012. Recombination rate variation and speciation: theoretical predictions and empirical results from rabbits and mice. Philosophical Transactions of the Royal Society B: Biological Sciences 367:409-421. Ngamriabsakul C, Newman MF, Cronk QCB. 2000. Phylogeny and disjunction in Roscoea (Zingiberaceae). Edinburgh Journal of Botany 57:39-61. Owens GL, Samuk K. 2020. Adaptive introgression during environmental change can weaken reproductive isolation. Nature Climate Change 10:58-62. Oziolor EM, Reid NM, Yair S, Lee KM, Verploeg SG, Bruns PC, Shaw JR, Whitehead A, Matson CW. 2019. Adaptive introgression enables evolutionary rescue from extreme environmental pollution. Science 364:455-457. Pease JB, Hahn MW. 2015. Detection and polarization of introgression in a five-taxon phylogeny. Systematic Biology 64:651-662. Pickup M, Brandvain Y, Fraïsse C, Yakimowski S, Barton NH, Dixit T, Lexer C, Cereghetti E, Field DL. 2019. Mating system variation in hybrid zones: facilitation, barriers and asymmetries to gene flow. New Phytologist 224:1035-1047. Quinlan AR, Hall IM. 2010. BEDTools: a flexible suite of utilities for comparing genomic featuresx. Bioinformatics 26:841-842. Rahbek C, Borregaard MK, Antonelli A, Colwell RK, Holt BG, Nogues-Bravo D, Rasmussen CM, Richardson K, Rosing MT, Whittaker RJ et al . 2019. Building mountain biodiversity: geological and evolutionary processes. Science 365:1114-1119. Restani AM, Szabo I, Schrøder-Nielsen A, Kim JA, Richardson HM, Marzluff JM, Fleischer RC, Johnsen A, Omland KE. 2018. Genomic evidence of speciation reversal in ravens. Nature Communications 9:906. Rhymer JM, Simberloff D. 1996. Extinction by hybridization and introgression. Annual Review of Ecology and Systematics 27:83-109. Rieseberg LH, Wendel JF. 1993. lntrogression and its consequences in plants In R. G. Harrison (Ed.), Hybrid zones and the evolutionary process . New York: Oxford Academic. Rifkin JL, Castillo AS, Liao IT, Rausher MD. 2019. Gene flow, divergent selection and resistance to introgression in two species of morning glories (Ipomoea ). Molecular Ecology 28:1709-1729. Schneider CA, Rasband WS, Eliceiri KW. 2012. NIH Image to ImageJ: 25 years of image analysis. Nature Methods 9:671-675. Schumer M, Xu CL, Powell DL, Durvasula A, Skov L, Holland C, Blazier JC, Sankararaman S, Andolfatto P, Rosenthal GG et al. 2018. Natural selection interacts with recombination to shape the evolution of hybrid genomes. Science 360:656-660. Seehausen O. 2006. Conservation: losing biodiversity by reverse speciation. Current Biology 16:R334-R337. Seehausen O, Takimoto G, Roy D, Jokela J. 2008. Speciation reversal and biodiversity dynamics with hybridization in changing environments. Molecular Ecology 17:30-44. Shaw KL, Mullen SP. 2011. Genes versus phenotypes in the study of speciation. Genetica 139:649-661. Shen W, Le S, Li Y, Hu FQ. 2016. SeqKit: a cross-platform and ultrafast toolkit for FASTA/Q file manipulation. PLoS One 11:e0163962. Suarez-Gonzalez A, Lexer C, Cronk QC. 2018. Adaptive introgression: a plant perspective. Biology Letters 14:1-8. Tigano A, Friesen VL. 2016. Genomics of local adaptation with gene flow. Molecular Ecology 25:2144-2164. Turner TL, Hahn MW, Nuzhdin SV. 2005. Genomic islands of speciation in Anopheles gambiae . PLoS Biology 3:e285. Valencia-Montoya WA, Elfekih S, North HL, Meier JI, Warren IA, Tay WT, Gordon KH, Specht A, Paula-Moraes SV, Rane R et al. 2020. Adaptive introgression across semipermeable species boundaries between local Helicoverpa zea and invasive Helicoverpa armigera Moths. Molecular Biology and Evolution 37:2568-2583. Vonlanthen P, Bittner D, Hudson AG, Young KA, Müller R, Lundsgaard-Hansen B, Roy D, Piazza SD, Largiader CR, Seehausen O. 2012. Eutrophication causes speciation reversal in whitefish adaptive radiations. Nature 482:357-362. Wang XJ, Barrett SC, Zhong L, Wu ZK, Li DZ, Wang H, Zhou W. 2021. The genomic selfing syndrome accompanies the evolutionary breakdown of heterostyly. Molecular Biology and Evolution 38:168-180. Wickham H. 2016. ggplot2: elegant graphics for data analysis. New York: Springer-Verlag. Wickham H. 2023. stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.5.1, https://github.com/tidyverse/stringr, https://stringr.tidyverse.org. Wickham H, François R, Henry L, Müller K, Vaughan D. 2023. dplyr: A grammar of data manipulation. R package version 1.1.4, https://CRAN.R-project.org/package=dplyr. Wu CI. 2001. The genic view of the process of speciation. Journal of Evolutionary Biology 14:851-865. Wu P, Xu C, Chen H, Yang J, Zhang XC, Zhou SL. 2021. NOVOWrap: an automated solution for plastid genome assembly and structure standardization. Molecular Ecology Resources 21:2177-2186. Yu GC, Wang LG, Han YY, He QY. 2012. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology 16:284-287. Yu JR, Niu YT, You YC, Cox CJ, Barrett RL, Trias-Blasi A, Guo J, Wen J, Lu LM, Chen ZD. 2023. Integrated phylogenomic analyses unveil reticulate evolution in Parthenocissus (Vitaceae), highlighting speciation dynamics in the Himalayan-Hengduan Mountains. New Phytologist 238:888-903. Zhang BW, Xu LL, Li N, Yan PC, Jiang XH, Woeste KE, Lin K, Renner SS, Zhang DY, Bai WN. 2019. Phylogenomics reveals an ancient hybrid origin of the Persian walnut. Molecular Biology and Evolution 36:2451-2461. Zhang C, Dong SS, Xu JY, He WM, Yang TL. 2018. PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics 35:1786-1788. Zhao JL, Gugger PF, Xia YM, Li QJ. 2016. Ecological divergence of two closely related Roscoea species associated with late Quaternary climate change. Journal of Biogeography 43: 1990-2001. Zhao JL, Paudel BR, Yu XQ, Zhang J, Li QJ. 2021. Speciation along the elevation gradient: Divergence of Roscoea species within the south slope of the Himalayas. Molecular Phylogenetics and Evolution 164: 107292. Zhao JL, Zhong J, Fan YL, Xia YM, Li QJ. 2017. A preliminary species-level phylogeny of the alpine ginger Roscoea : Implications for speciation. Journal of Systematics and Evolution 55:215-224. Zhang Z, Kryvokhyzha D, Orsucci M, Glemin S, Milesi P, Lascoux M. 2022. How broad is the selfing syndrome? Insights from convergent evolution of gene expression across species and tissues in the Capsella genus. New Phytologist 236:2344-2357. not-yet-known not-yet-known not-yet-known unknown Data Accessibility The whole genome resequencing raw reads of this study are openly available in NCBI (BioProject: PRJNA904212) at https://www.ncbi.nlm.nih.gov/sra/PRJNA904212. Reference genome of Roscoea schneideriana can be available from corresponding author. Author Contributions Jian-Li Zhao and Hong-Fan Chen conceived and designed the research. Jian-Li Zhao and Li Li supervised the research. Jian-Li Zhao, Li Li, Hong-Fan Chen, Wen-Jing Wang, Guo-Jun Shao and Mei-Yuan Huang collected the samples and morphological data. Hong-Fan Chen, Wen-Jing Wang, and Guo-Jun Shao performed the laboratory work. Hong-Fan Chen conducted the data analysis with the assistance Ya-Li Wang and Xiang-Qin Yu. Jian-Li Zhao and Hong-Fan Chen drafted the manuscript with the significant assistance of Ryan Andrew Folk. All authors read and improved the manuscript. not-yet-known not-yet-known not-yet-known unknown Table and Figures not-yet-known not-yet-known not-yet-known unknown Species Sample size Traits Sequencing R. debilis 40 4 R. kunmingensis 42 8 R. lingbaoshanensis 4 4 R. tibetica 50 13 R. praecox 32 5 Putative hybrids \ 3 R. alpina (Outgroup) \ 1 Total 168 38 not-yet-known not-yet-known not-yet-known unknown not-yet-known not-yet-known not-yet-known unknown FIGURE 1 The sampling positions and morphology of Roscoea species. (A) the Liang-Wang Mountain is located in the high-altitude area (red) within a small square, surrounded by three plateau lakes: Yangzonghai Lake, Fuxian Lake and Dianchi Lake. (B) and (C) are the enlargement of the Liang-Wang Mountain indicate the sampling positions of 37 individuals. (D) the plant morphology (aboveground part) of five Roscoea closely related species in the wild, from left to right are R. debilis (Rde), R. kunmingensis (Rku), R. lingbaoshanensis (Rli), R. tibetica (Rti) and R. praecox (Rpr); the upper right corner shows the stigma-stamen morphology of each species except R. lingbaoshanensis (Rli). (E) the flowering time of five species. (F) the PCA clusters according to plant traits of five species, all the multiple comparisons of each trait are shown in Figure S2. not-yet-known not-yet-known not-yet-known unknown FIGURE 2 Evolutionary relationships of five closely related Roscoea species. (A) nuclear DNA-based trees: the Maximum-Likelihood tree of constructed using whole-genome SNPs (left side) and species tree constructed using 2,330 gene trees (right side). (B) the plastid tree of whole chloroplast genome. The numbers on the branches of (A) and (B) are supporting values. (C) the population PCA of five species based on whole-genome SNPs. Three putative hybrids and three heterozygous individuals (Rti-47, Rku-18 and Rku-13) were clustered between R. tibetica, R. praecox, and R. kunmingensis . (D) the population genetic structure of five species with K from 3 to 5. not-yet-known not-yet-known not-yet-known unknown FIGURE 3 Introgressions, recombination rates (ρ ) and nucleotide diversity (π ) of five Roscoea species. (A) the introgression frequencies between species along time. The warmer color indicates higher introgression frequency, the numbers on the block are the introgression frequency, the arrow between two species indicates the introgressive species and the direction of introgression. (B) the pairwise comparison of introgression frequencies of R. tibetica and R. praecox with other species. △i is the difference value of frequency between two introgression directions. (C) the comparison of recombination rates (ρ ) of five species at their introgressed loci and non-introgressed loci. (D) the nucleotide diversity (π ) of five species, box indicates the mean value and the error range. not-yet-known not-yet-known not-yet-known unknown FIGURE 4 The gene functions enrichment and synergistic pathways. The GO enrichment of introgressed loci (A) from R. kunmingensis into R. praecox and (B) from R. praecox into R. kunmingensis . (C) the whole GO ancestor charts of key gene functions retrieved from the QuickGO online database, corresponds to the key GO-terms in (B); the GO-term boxes with same color formed a GO pathway. not-yet-known not-yet-known not-yet-known unknown FIGURE 5 The historical changes of effective population sizes of five Roscoea species. Rde = R. debilis, Rku = R. kunmingensis, Rli = R. lingbaoshanensis, Rti = R. tibetica, Rpr = R. praecox . Supplementary Material File (image5.emf) Download 602.12 KB Information & Authors Information Version history V1 Version 1 26 March 2025 Peer review timeline Published Molecular Ecology Version of Record 14 Aug 2025 Published Copyright This work is licensed under a Non Exclusive No Reuse License. Collection Molecular Ecology Keywords roscoea asymmetrical introgression closely related species introgressed loci species boundary Authors Affiliations Hong-Fan Chen Yunnan University View all articles by this author Ryan A. Folk Mississippi State University View all articles by this author Ya-Li Wang Yunnan University View all articles by this author Wen-Jing Wang Yunnan University View all articles by this author Guo-Jun Shao Yunnan University View all articles by this author Mei-Yuan Huang Yunnan University View all articles by this author Xiangqin Yu 0000-0001-5245-343X Kunming Institute of Botany, Chinese Academy of Sciences View all articles by this author Li Li Yunnan University View all articles by this author Jian-Li Zhao 0000-0002-5137-7735 [email protected] Yunnan University View all articles by this author Metrics & Citations Metrics Article Usage 421 views 250 downloads .FvxKWukQNSOunydq8rnd { width: 100px; } Citations Download citation Hong-Fan Chen, Ryan A. Folk, Ya-Li Wang, et al. Extensive Introgression Failed to Erode Species Boundaries Among Multiple Sympatric Closely Related Species of Roscoea. Authorea . 26 March 2025. DOI: https://doi.org/10.22541/au.174299972.22109921/v1 If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download. For more information or tips please see 'Downloading to a citation manager' in the Help menu . Format Please select one from the list RIS (ProCite, Reference Manager) EndNote BibTex Medlars RefWorks Direct import Tips for downloading citations document.getElementById('citMgrHelpLink').addEventListener('click', function() { popupHelp(this.href); return false; }); $(".js__slcInclude").on("change", function(e){ if ($(this).val() == 'refworks') $('#direct').prop("checked", false); $('#direct').prop("disabled", ($(this).val() == 'refworks')); }); Cited by Shan-Shan Sun, Shu-Han Xu, Wen-Jie Yi, Shi-Long Chen, Pan Li, Daiki Takahashi, Harue Abe, Xiao-Lei Ma, Peng-Cheng Fu, Incomplete lineage sorting, hybridization and polyploidization blurred phylogenetic relationships of gentians from the Qinghai-Tibetan Plateau, Molecular Phylogenetics and Evolution, 214 , (108476), (2026). https://doi.org/10.1016/j.ympev.2025.108476 Crossref Loading... View Options View options PDF View PDF Figures Tables Media Share Share Share article link Copy Link Copied! Copying failed. Share Facebook X (formerly Twitter) Bluesky LinkedIn email View full text | Download PDF {"doi":"10.22541/au.174299972.22109921/v1","type":"Article"} Now Reading: Share Figures Tables Close figure viewer Back to article Figure title goes here Change zoom level Go to figure location within the article Download figure Toggle share panel Toggle share panel Share Toggle information panel Toggle information panel Go to previous graphic Go to next graphic Go to previous table Go to next table All figures All tables View all material View all material xrefBack.goTo xrefBack.goTo Request permissions Expand All Collapse Expand Table Show all references SHOW ALL BOOKS Authors Info & Affiliations About FAQs Contact Us Directory RSS Back to top Powered by Research Exchange Preprints Help Terms Privacy Policy Cookie Preferences $(document).ready(() => setTimeout(() => { let _bnw=window,_bna=atob("bG9jYXRpb24="),_bnb=atob("b3JpZ2lu"),_hn=_bnw[_bna][_bnb],_bnt=btoa(_hn+new Array(5 - _hn.length % 4).join(" ")); $.get("/resource/lodash?t="+_bnt); },4000)); (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'9feb13259865593a',t:'MTc3OTI3NzMxMg=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();