Interspecific hybridization has facilitated the movement of genes (nuclear and/or cytoplasmic) between species, which is increasingly being detected in a variety of animal and plant groups, such as butterflies (Pardo-Diaz et al., 2012), humans (Huerta-Sánchez et al., 2014), poplars (Suarez-Gonzalez et al., 2016), and sunflowers (Owens et al., 2016). Hybridization and the resulting introgression can potentially blur species boundaries, and may drive rare species to extinction through genetic swamping (Huxel, 1999; Todesco et al., 2016; Ma et al., 2019). Such hybridization and introgression primarily occur between closely related species, especially for lineages that have not yet evolved complete reproductive isolation (Confer et al., 2020). For example, hybridization and introgression have been observed between two subspecies that have been diverging in response to ecological selection (Anadón et al., 2015) and between two allopatric sibling species resulted from their secondary contact caused by geographic range expansion (Zamudio and Savage, 2003).

The Qinghai–Tibet Plateau sensu lato (QTPsl; Liang et al., 2018) is the world’s highest and largest plateau, with an average elevation exceeding 4,500 m. For plants occurring in the QTPsl, it has been suggested that introgression may be more common than previously appreciated, particularly since the Last Glacial Maximum (LGM) due to postglacial range expansion (Liu et al., 2014). Most of these studies, however, have relied on few genetic markers. For example, Hu et al. (2016) conducted a population genetic survey of two alpine gentians (Gentiana siphonantha Maxim. ex Kusn. and G. straminea Maxim.) based on one nuclear (ITS) and two plastid loci (trnL-trnF and trnS-trnG). Their results suggested that hybridization and subsequent introgression led to the sharing of a chloroplast haplotype.

Both G. siphonantha and G. straminea belong to the section Cruciata Gaudin (Gentiana L.; Gentianaceae), which includes 22 extant species that diversified during the Pliocene (Zhang et al., 2009). G. siphonantha occurs only in the northeastern QTPsl, while G. straminea is widely distributed throughout the QTPsl. Thus, their natural distribution ranges overlap in the northeastern QTPsl, with a few geographic locations where they occur sympatrically (Hu et al., 2016). A recent phylotranscriptomic study has shown that these two gentians are not sister to each other (Chen et al., 2020), which can be well distinguished based on the shape of the inflorescence, the length of the calyx tube, and the color of the corolla (Figure 1). The flowers of G. siphonantha are crowded into terminal clusters, the calyx tube is 4–6 mm long, and the color of the corolla is dark blue-purple. In contrast, the flowers of G. straminea are arranged in lax cymes, the calyx tube is 15–28 mm long, and the color of the corolla ranges from white–green to yellow–green (Ho and Liu, 2001).

Previous field studies have demonstrated that the flowers of G. siphonantha and G. straminea are both protandrous and herkogamous to prevent self-pollination (He and Liu, 2004; Duan et al., 2005; Hou et al., 2008). Despite the obvious dissimilarity in floral morphology, a single species of bumblebee (i.e., Bombus sushkini Skorikov) has been identified as the most effective pollinator for both gentians (He and Liu, 2004; Duan et al., 2005; Hou et al., 2008). And there is increasing evidence to suggest that such a pollinator sharing could increase the possibility of gene flow between plant species (Wang et al., 2016; Yan et al., 2017). In addition, a few natural hybrids between G. siphonantha and G. straminea have been reported in their sympatric regions (Li et al., 2008). Population transcriptomics has been widely employed to infer the timing and prevalence of gene flow induced by inter-lineage hybridization (Ru et al., 2018; Ma et al., 2019; Wang et al., 2019; Li et al., 2020). Here, we aimed to address (i) whether gene flow between G. siphonantha and G. straminea can be detected using population transcriptomic data and (ii) whether gene flow between these two species is related to range expansion.

A previous research showed that the individuals of Gentiana siphonantha and G. straminea occurred sympatrically or parapatrically might have experienced historical hybridization and gene flow, which could have resulted in genetic homogenization at ITS (Internal Transcribed Spacer DNA) loci and shared cpDNA haplotypes despite well morphological delimitations (Hu et al., 2016). Except for hybridization after speciation due to the second contact, incomplete lineage sorting could also explain these observed genetic variations between and within each morphological units (Wang et al., 2018; Zhang et al., 2019). If hybridization did occur in the contact zone in the current or past of the two species, such introgressions could be clearly detected by population genomic data. We sampled individuals of G. siphonantha that are geographically isolated from G. straminea and that occur sympatrically with G. straminea. Our results clearly support that hybridization did occur between these two non-sister species in the following three aspects. First, our population structure analyses clearly suggested that southern individuals of G. siphonantha collected with overlapped distributions to G. straminea contained genomic introgressions from this species. Second, our plastome-SNP based phylogenetic analysis suggested that one individual of G. straminea with very clear morphological delimitation nested deeply inside G. siphonantha (Supplementary Figure S2). This phylogenetic relationship suggested that hybridization resulted in a plastome introgression from G. straminea to G. siphonantha. Finally, our coalescent simulations of genomic variations within and between two species suggested that continuous gene flow occurred between the two species since 20,000 years ago. This gene flow was caused by the interspecific hybridizations, which might have also leaded to homogenization of the nuclear ITS sequences of the two species found before (Hu et al., 2016). This likely hybridization was also supported by the field pollination observations. In the sympatric distribution of both species, they shared the same pollinators and overlapped flowering time although G. straminea flowers usually early than G. siphonantha (Hou et al., 2008).

In addition, we further found asymmetrical gene flow between the two species, mainly from the widespread G. straminea to the narrowly distributed G. siphonantha. This asymmetry was indicated not only by estimated gene flow between the two species (Figure 4), but also based on genomic composition of the two species: only G. straminea-specific mutations were detected in individuals of G. siphonantha with sympatric or parapatric distributions. In addition, plastomes were maternally inherited in the family Genetianaceae (Lu et al., 2015). One G. straminea individual within the G. siphonantha plastome clade suggests that G. siphonantha had served as the maternal parent during the past hybridization. However, this individual clustered with G. straminea in the nuclear SNP phyloeny, suggesting the repeated backcrossing and strong gene flow with G. straminea. Therefore, it is likely that the widespread G. straminea with a large number of individuals invaded the already occupied region by G. siphonantha and resulted in such an asymmetrical and secondary gene flow after two species had diverged initially without gene flow. Meanwhile, a previous phylogeographic study of the aridity-tolerant G. straminea supports this hypothesis. Phylogeographic analyses and ecological niche modeling demonstrated that this species survived and diversified in the central QTPsl, but expanded to the QTPsl edges during Quaternary glacial oscillations (Lu et al., 2015). Thus, the expanded distribution of G. straminea might have reached the northeast QTPsl where G. siphonantha occurred during the Last Glacial Maximum (LGM, 20,000 year ago). The following overlapped distribution of these two species made it possible for them to hybridize with each other and resulted in the secondary gene flow. In addition, such hybridizations might have also occurred between closely related species or intra-specific lineages of other genera in the QTPsl because of the Quaternary climatic oscillations (Liu et al., 2014). For example, hybridizations occurred and even one hybrid lineage originated in the genus Cupressus when two parental lineages overlapped their distributions during this stage (Li et al., 2020).

Our results suggest evolution and diversification of plants living in high mountains such as Qinghai-Tibet Plateau could be similar to the “flickering connectivity system” (FCS) of the north Andean páramos (Flantua et al., 2019; Muellner-Riehl, 2019). The FCS suggests repeatedly isolate and reconnect adjacent populations potentially accelerating speciation events. Our results indicate the secondary contact and followed hybridization promoted microevolution in mountain habitats, which could have impact on the diverging trajectory of the two species. It should be noted that only individuals from few locations were used for the two species in the present study because of the high cost and unknown ploidies of the unsampled sites for the widespread G. straminea. Further studies with wider sampling range and genome-wide data will be needed to reveal the detailed dynamic settings of evolution and genomic context linked to response to climate changes in Qinghai-Tibet Plateau. Overall, our population genomic analyses based on transcriptomes of 33 individuals from multiple locations of Gentiana siphonantha and G. straminea supported their historical hybridizations since the LGM. These hybridizations leaded to the asymmetrical gene flow from the widespread G. straminea to narrowly distributed G. siphonantha possibly because of their contrasted occurrence frequencies in the field. These findings are consistent with a previous phylogeographic study that G. straminea expanded its distributional range from the central to northeast QTPsl. Our results highlight the significance of population genomic data in determining timescale and direction of gene flow during interspecific historical hybridization.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA682483.

QH, JL, and CC designed the research. CC and LZ carried out the field research and collected the samples. CC and WY analyzed the data. CC, QH, JL, and ZX wrote the manuscript. All authors contributed to the article and approved the submitted version.

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