A Combined Morphological and Molecular Evolutionary Analysis of Karst-Environment Adaptation for the Genus Urophysa (Ranunculaceae)

The karst environment is characterized by low soil water content, periodic water deficiency, and poor nutrient availability, which provides an ideal natural laboratory for studying the adaptive evolution of its inhabitants. However, how species adapt to such a special karst environment remains poorly understood. Here, transcriptome sequences of two Urophysa species (Urophysa rockii and Urophysa henryi), which are Chinese endemics with karst-specific distribution, and allied species in Semiaquilegia and Aquilegia (living in non-karst habitat) were collected. Single-copy genes (SCGs) were extracted to perform the phylogenetic analysis using concatenation and coalescent methods. Positively selected genes (PSGs) and clusters of paralogous genes (Mul_genes) were detected and subsequently used to conduct gene function annotation. We filtered 2,271 SCGs and the coalescent analysis revealed that 1,930 SCGs shared the same tree topology, which was consistent with the topology detected from the concatenated tree. Total of 335 PSGs and 243 Mul_genes were detected, and many were enriched in stress and stimulus resistance, transmembrane transport, cellular ion homeostasis, calcium ion transport, calcium signaling regulation, and water retention. Both molecular and morphological evidences indicated that Urophysa species evolved complex strategies for adapting to hostile karst environments. Our findings will contribute to a new understanding of genetic and phenotypic adaptive mechanisms of karst adaptation in plants.

Introduction

Environmental heterogeneity is one of the most important factors influencing evolutionary trajectories and ecological adaption of species (Allouche et al., 2012). Plant species have adapted to pronounced gradients of environmental conditions (temperature, drought, oxidative, or osmosis stresses etc.) along with variability in their natural habitats. Species survival is dependent on maximizing their fitness in all environmental conditions they encounter (Pratlong et al., 2015). Therefore, understanding the drivers and genetic mechanisms of species persistence and adaptation is a basis for the study of ecological adaption, evolution, and speciation (Schnitzler et al., 2011).

The area of karst in southern and southwestern China boasts over 20,000 plant species, and its flora is ranked as one of the most endemic-rich subtropical flora in the world (Kang et al., 2014). This region is not only one of the most famous models of karst landform in the world with an area of 6.2 × 10⁵ km² (Yuan, 1994; Su et al., 2002; Huang Q. et al., 2008), but also has the most varied of extreme karst environments characterized by low soil water content, periodic water deficiency and poor nutrient availability, and may have exerted strong selective forces on plant evolution (Ai et al., 2015). Although the ecological environment in karst region is extremely hostile, plant communities in this region exhibit remarkably high levels of species richness and endemism and make a large contribution to the floristic diversity of China (Gao et al., 2015; Hao et al., 2015; Feng et al., 2020). Moreover, this region possesses various ecological niches afforded by complex terrains (e.g., fissured cliffs and extensive caves) and variable climatic and edaphic conditions (Hao et al., 2015). Soils in this area are also typically shallow and are characterized by higher concentrations of calcium (Ca²⁺) and magnesium (Mg²⁺), higher pH levels, and a lower water storage capacity compared to non-karst soils in other subtropical or tropical regions (Nie et al., 2011; Luo et al., 2012). Many calcicoles (species adapted to calcareous soil) have evolved in the region, and their habitats and niches have long been regarded as ‘natural laboratories’ for ecological and evolutionary studies to understand natural selection and species evolution because of the high diversity and unique biota (Clements et al., 2006). Unfortunately, previous studies on karst-specific species are scarce (Kang et al., 2014; Ai et al., 2015; Li et al., 2015; Liu et al., 2019; Feng et al., 2020), and until now, the mechanisms of karst environmental adaptation have rarely been examined. Yet such information is important for understanding the survival and adaptive strategies of species in the special karst environmental conditions.

One particularly interesting genus endemic to the karst region is the Urophysa Ulbr., which belongs to the family Ranunculaceae. The genus consists of two species, Urophysa rockii Ulbrich and Urophysa henryi (Oliver) Ulbrich, which are distributed in allopatric regions and are adapted to remarkably uniform habitats with specific calcareous soils developed from karst limestone bedrock in southwest China (Figure 1) (Wang et al., 2004; Xie et al., 2017). U. rockii and U. henryi exhibit distinct petal morphologies: the former has sacs near the base of petals while the latter has not (Xie et al., 2017). Among them, U. rockii (also be called the Panda Grass) was first discovered in 1925 but was not found until Dr. Chunyu Li rediscovered it in Jiangyou County of western Sichuan Province in 2005 (Du et al., 2010). There are only approximately 2,000 individuals in five populations currently surviving in upstream of the Fujiang River (Wang and He, 2011; Xie et al., 2017). Hence it was listed as an endangered species and included on the China Species Red List. Only a few populations of U. henryi have been found and most are separated by high mountains and deep valleys in karst regions of southern China (Xie et al., 2016, 2017). Interestingly, the two species all flower in the winter (flourishing florescence is between December to March), and in field observations and laboratory experiments, we found that U. rockii and U. henryi can not survive outside the karst limestone, which indicated that the karst limestone plays a significant role in their growth and development. However, there have been no studies investigating the biological adaptation mechanisms in the two species, but there have been studies of their genetic diversity and differentiation (Xie et al., 2016, 2017), growing environment and conservation strategies (Liu Y.Q. et al., 2009; Du et al., 2010), and biological and ecological characteristics (Wang and He, 2011; Zhang L. et al., 2013; Zhang et al., 2013a, b). More comprehensive studies are required because populations of the two species have been lost, degraded and fragmented by human development activities (e.g., scenic spot construction, hydroelectric stations) and excessive exploitation due to their medicinal value (contusion and bruise treatment). Therefore, before more populations are lost, it is crucial to explore the influence of the karst-specific habitat on the evolution and growth of these two species.

FIGURE 1

Figure 1. The habitat and morphological characters of Urophysa rockii (A–D) and Urophysa henryi (E–H). (A,E) Habitat feature; (B,F) growth sites; (C,G) the flowers; (D,H) fruits.

With the rapid development of sequencing technologies, the whole-genome and transcriptome sequencing have become increasingly popular analytical techniques that can investigate various adaptive processes, such as in response to highland, salt, drought, and temperature stress (Li et al., 2014; Zhang et al., 2016, 2019, 2020; Mock et al., 2017; Wang K. et al., 2019; Zeng et al., 2019; Birkeland et al., 2020; Xie et al., 2020). However, the vast majority have been performed on animals and limited genome/transcriptome-based research has been devoted to plants in the karst region. Previous studies on the genus Primulina which has been studied the most in karst limestone habitats, indicated that genome size is phylogenetically conserved but its variation among species is a combination of both neutral and adaptive evolution (Kang et al., 2014; Feng et al., 2020). Tao et al. (2016) suggested that the Ca²⁺-permeable channel TPC1 may be involved in the local adaptation of Primulina to karst Ca²⁺-rich environments, and the potential positive selection in phytochrome PHYE may play an important role in the adaption of Primulina to the heterogeneous environment of karst (Tao et al., 2015). Nevertheless, we know little about the genetic bases of plants adapted to the karst environment.

In this study, we performed RNA-seq to obtain most transcriptome sequences of the karst species U. rockii and U. henryi and non-karst species in genus Semiaquilegia and Aquilegia. Single-copy genes were extracted for phylogenetic analysis, and by positive selection analysis and function annotation, we aimed to investigate the potential mechanisms of how these two-species adapted to the karst environment. This study represents the first exploration of Urophysa transcriptomes using the high-throughput RNA-seq, and will help guide further molecular systematics, population genetic, and ecological adaption studies in Urophysa and other related species.

Materials and Methods

Morphological Character Investigation, Sample Collection, and Transcriptome Sequencing

To understand the morphological characters and living conditions of the karst and non-karst species, we performed in situ and ex situ field investigation in the Sichuan, Hubei and Shanxi Provinces of China. The morphological traits and habitats of each species were recorded, including the flowering time, flower and leaf traits and color variation, seeds dispersal and living environment condition. The young and mature leaves and flowers of U. rockii and U. henryi as well as Semiaquilegia adoxoides and Aquilegia ecalcarata were collected from living plants, and stored at −80°C until processed for total RNA isolation. Total RNA was extracted using Trizol (Life Technologies Corp.) according to the manufacturer’s protocols, and cleaned up using the RNeasy mini kit (Qiagen, Valencia, CA, United States). Poly (A) mRNA was purified using Oligo (dT) magnetic beads and interrupted into short fragments. The RNA sequencing libraries were constructed using the Illumina mRNA-Seq Prep Kit after assessing RNA quality. Finally, all cDNA libraries were sequenced on the Illumina HiSeq 2000 platform in paired-end form. A total of eight sequencing libraries were constructed in this study, including each flower and leaf libraries of the four species.

De novo Assembly, Completeness Test, and Functional Annotation

Quality of clean data (clean reads) was assessed using the FastQC version 0.11.8¹ by removing sequencing adaptors and reads with unknown nucleotides and low quality (quality scores < 20). All subsequent analyses were based on these filtered reads. Transcriptome de novo assembly was performed using Trinity software with default parameters (Grabherr et al., 2011). To evaluate the completeness of the assembly results, we applied BUSCO version 3.0.2 (Simão et al., 2015) to assess the quality of assembled transcriptome data based on the embryophyte gene database (embryophyte_odb 10.2020-09-10)². Only contigs with lengths greater than 200 bp were kept for further analysis. Then, CD-HIT-EST version 4.6.8 (Li and Godzik, 2006) was used to remove redundant contigs with a threshold of 1, and TransDecoder version 0.36³ was employed to predict the open reading frame (ORF) with a minimum protein length of 100 bp.

Functional annotations of all assembled protein sequences were performed by searching against the following databases: Flowering plant in NCBI Non-redundant protein sequence database (NR, Taxonomy ID is 3398), Protein family (Pfam), Clusters of Orthologous Groups of proteins (COG), Swiss-Prot protein (Swiss-Prot), KEGG ortholog database (KEGG), and Gene Ontology (GO). Among them, sequences were searched against the flowering plant database, Pfam, COG, and Swiss-Prot using BLASTP version 2.2.31 with an e-value cut-off of 1 × e^–5. GO and KEGG annotations were performed by Blast2GO software (Conesa et al., 2005) with an e-value cutoff of 10^–5 and then plotted with functional classification using Web Gene Ontology Annotation Plot (WEGO) (Ye et al., 2006).

Orthologous and Paralogous Genes Identification and Phylogenetic Analysis

In order to identify the shared orthologs and then used for further analyses, genome data of Aquilegia coerulea were downloaded from Phytozome⁴ and the OrthoMCL version 2.0.9 software (Li et al., 2003) was used to extract the orthologous genes of the five species with e-value of 1e^–10 and inflation = 1.5. The one to one ortholog predictions were conducted according to the following criteria: (1) Muscle (Edgar, 2004) was used to align the sequences of each orthologous group from the flow of Agalma (Dunn et al., 2013); (2) GBlocks (Talavera and Castresana, 2007) was used to filter the conserved regions of each sequence; (3) the ML tree was built in IQ-TREE 2 (Minh et al., 2020) by using above filtered conserved sequences with parameters set as: iqtree -s ^∗.phy -m MFP -bb 1,000 -bnni -cmax 20 -redo; (4) the ML tree was trimmed by implementing the Dendropy (Sukumaran and Holder, 2010) and only one copy of each species was finally obtained in each orthogroup. Additionally, PAL2NAL version 14 (Suyama et al., 2006) was employed to map the nucleotide sequences of each ortholog to the corresponding protein sequences and remove the highly variable nucleotide sites to obtain the final nucleotide sequences. A total of 2,271 single-copy genes (SCGs) were generated and concatenation and coalescent methods were implemented for phylogenetic analysis. Meanwhile, 243 clusters of paralogous genes were detected. For detailed analysis procedures, please see Figure 2.

FIGURE 2

Figure 2. The detailed analysis procedures performed in this study. (A) The phylogenetic analysis workflow including concatenate and coalescent methods; (B) procedures of positive selection analysis.

For the concatenation method, all SCGs were concatenated as a super data matrix after alignment and manual checking. Two conventional approaches were used to perform the phylogenetic analysis: (1) A Maximum Likelihood (ML) tree was reconstructed using IQ-TREE 2 (Minh et al., 2020), with the parameters applied as: iqtree -s ^∗.phy -m MFP -bb 1,000 -bnni -cmax 20 -redo. (2) A Bayesian inference (BI) analysis was performed in MrBayes v3.2 (Ronquist and Huelsenbeck, 2003) with the best-fit model GTR + G being selected under the Akaike information criterion (AIC). Four independent Markov chain Monte Carlo (MCMC) runs were conducted, each chain run for 5 × 10⁸ generations using initial trees with every 1,000 generations being sampled, and the initial 20% of the samples were discarded as burn-in to confirm the stationarity. Tracer v1.5 (Rambaut and Drummond, 2009) was used to assess the quality of the MCMC simulations and the stability of runs. Effective sample size (ESS) values were greater than 200 for all parameters, indicating sufficient sampling. For the coalescent method, we performed the ML analysis for each of 2,271 SCGs using IQ-TREE 2 (Minh et al., 2020) with the same parameters and processes as in the concatenation method. Urophysa species were rooted as outgroup based on previous studies (Xie et al., 2017; Zhai et al., 2019; Duan et al., 2020). Furthermore, in order to improve the average positive branch rates and achieve believable phylogenetic results, the programe STAR (Liu L. et al., 2009) and MP-EST (Liu et al., 2010) were used to estimate the species tree based on 2,271 gene trees. Finally, the concatenate and coalescent phylogenetic results were used for species tree derivation.

Positive Selection Analysis

To detect positive selection on genes along a specific lineage, the optimized branch-site model (BSM) (Yang and Dos, 2011) in the Codeml program of the PAML 4 package (Yang, 2007) was performed. This model was conducted based on (1) a confirmed phylogenetic tree; (2) gene sequences with Paml format; (3) parameter-setting file; (4) Bayesian empirical Bayes (BEB) test and (5) Chi-square test. In this model, branches in the tree were divided a priori into foreground and background categories, and a likelihood ratio test was constructed to compare an alternative model that allows for some sites under positive selection on the foreground branches with a null model that does not. Tree topology was analyzed based on the concatenate and coalescent phylogenetic results, and the final species tree was used for the next analysis. The branch with U. rockii and U. henryi was selected as the foreground branch, and other branches were as a background branch. The rate (ω) of the non-synonymous substitution rate (dN) to synonymous substitution rate (dS) was used to measure the selective pressure. ω > 1 implies positive selection, ω = 1 means neutral selection and ω < 1 indicates negative selection (Yang and Nielsen, 2002). A likelihood ratio test (LRT) was constructed to compare the alternative model that allows sites to vary on the foreground branch with a null model that confined codon sites under neutral selection. The BEB method (Yang et al., 2005) was implemented to estimate the posterior probabilities for codon sites, and sites with a posterior probability larger than 0.9 were regarded as positively selected sites. The p-values were calculated with a Chi-square test and adjusted by the FDR method, and genes with p-values less than 0.05 were considered as positively selected genes (PSG). The Jalview software (Clamp et al., 2004) was used to view the amino acid sequences of some PSGs by highlighting the positively selected sites. The SWISS-MODEL (Schwede et al., 2003) was used to predict the three-dimensional protein structures and were further visualized with PyMOL software⁵. Furthermore, in order to detect more positive selected genes that involved in karst environment adaptation, a adaptive evolution was estimated by pairwise calculation of K_a/K_s between all members of each clusters paralogous genes. Two paris of models: M1 (neutral) vs. M2 (selection) and M7 (beta) vs. M8 (beta + ω) were used in PAML package. For each pair of nested model the log likelihood values are compared using the likelihood ratio test (LRT).

Gene ontology (GO) enrichment analyses were conducted with detected Urophysa PSGs and adapted paralogous genes to infer their functional information using the DAVID function toll (Huang D. W. et al., 2008). The KOBAS software (Xie et al., 2011) was used to test the statistical enrichment of PSGs in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Ogata et al., 2000). Moreover, the Cluster of Orthologous Groups of proteins (COG) database (Tatusov et al., 2000) was also used to predict the function of these PSGs. We then used a BLASTX search of these PSGs and Mul_genes with Arabidopsis thaliana to explore their homologies and functions.

Results

Morphological Characteristics and Growth Habit of Urophysa Species

U. rockii and U. henryi inhabit the karst cliffs and bloom in the winter, with the flourishing florescence between December to March. These species’ leaves possess obvious cuticular wax and the leaves color changes with the stages of development, which is green in the early stage of growth and purple or red in late development. Although these two species are herbaceous, they have long and stout rhizomes and are deeply rooted in the rock cracks. During field observations, we found that the pedicels of these two species bend toward the rock cracks in the fruiting period, and disperse seeds into the cracks of the karst rock. In contrast, the non-karst species S. adoxoides, A. ecalcarata, and A. coerulea grow in forest or roadside, although the three species also have rhizome (Tamura, 1995; Li, 2006). Leaves of the three species have not cuticular wax, and flowering times are from March to June and their seeds are mainly dispersed by gravity.

Assembly Statistics and Function Annotation

In this study, a total of 46.19 million clean reads were generated, and a total of 702,823 contigs were collected in the non-redundant assembled transcriptomes (Table 1). There were 75,477 and 85,640 unigenes detected in U. rockii and U. henryi and the N50 value was 1,715 and 1,652 bp, respectively. The number of unigenes found in A. ecalcarata and S. adoxoides were 69,978 and 100,099 with the length of N50 ranging from 1,655 to 1,582 bp. We detected more unigenes and longer N50 length compared to other species in Ranunculaceae, such as Helleborus thibetanus, Ranunculus spp. and Coptis chinensis (Chen et al., 2015, 2017; Shi et al., 2019). The assembly completeness evaluation results showed that all assembled transcriptomes had high BUSCO completeness, which varied from 86.73 to 91.59% in the four species (Supplementary Figure 1). The clean data were submitted to the NCBI Sequence Reads Archive (SRA) database (accession number: SRR6816169, SRR6816171-SRR6816172, SRR6816174, SRR12805323, and SRR12805326). By annotation, there were 37,200 (49.29%, U. rockii) and 40,493 (47.28%, U. henryi) unigenes matches on the Flowering plant database, while there were 37,121 (53.05%) and 58,408 (58.35%) in A. ecalcarata and S. adoxoides, respectively. Additionally, there were 41.46–47.41% unigene matches on Swiss-Prot database and ranged from 26.69 to 28.03% on the Pfam database, please refer to Table 2 for detailed annotation information for the four species.

TABLE 1

Table 1. Overview of the de novo assembly of the transcriptome of all species.

TABLE 2

Table 2. Statistics of annotations for assembled.

Identification of Orthologs, Paralogs, and Phylogenetic Analysis

From the orthologous analysis, a total of 40,487 putative orthologs were identified from the four species and A. coerulea based on OrthoMCL. There were 12,293 orthologs shared by all five species, and 4,588 orthologs were only found in the two Urophysa species (Figure 3). According to GO enrichment analysis, many orthologs shared by U. rockii and U. henryi were involved in response to stimuli, transport, ion binding, transporter activity, translation, membranes and others (Supplementary Figure 2). Among them, several metabolic processes were significantly enriched, for example, translation, RNA modification, proteolysis, DNA recombination, cytoplasmic translation, regulation of calcium-mediated signaling, nutrient reservoir (Supplementary Figure 3), and especially for the regulation of transmembrane transporter activity and calcium-mediated signaling (Figure 4 and Supplementary Figure 4). Meanwhile, 811 orthologs were exclusively found in the three non-karst species. All orthologs were non-significantly enriched and most of them were involved in metabolic process, molecular process, and macromolecule metabolic process (Supplementary Figure 5). In addition, 243 clusters of paralogous genes were detected among the five species, and by function annotation, many paralogs were related to transmembrane transport, salicylic id mediated signaling pathway, gene silencing by RNA, ATPase tivity (Figure 5 and Supplementary Table 1).

FIGURE 3

Figure 3. Venn diagram showing the numbers of orthologs identified in genus Urophysa and close species using OrthoVenn2 (https://orthovenn2.bioinfotoolkits.net/home).

FIGURE 4

Figure 4. Transmembrane transporter activity related pathway detected from the positively selected genes of Urophysa.

FIGURE 5

Figure 5. Function annotation for the 30 clusters of paralogous genes under significant adaptive evolution. (A–E) indicate the detail information of U. rockii, U. henryi, S. adoxoides, A. ecalcarata and A. coerulea.

A total of 2,271 SCGs were filtered from the five species and used for phylogenetic analysis. Alignments of concatenated SCGs were 2,229,258 bp in length, including 126,095bp (5.65%) variable sites and 31,439bp (1.41%) parsimony-informative sites. ML and BI phylogenetic analyses showed that S. adoxoides was firstly differentiated when rooted Urophysa species as outgroup, and A. ecalcarata tightly cluster with A. coerulea (Figure 6A). The coalescent analysis detected 2,271 gene trees, and the final species tree was consistent with the result from concatenate method (Figure 6B). Phylogenetic signal exploration of each gene tree topology found that 1,930 gene trees (84.95%) shared the same topology (Figure 6B1), which was consistent with the topology of species tree (Figures 6A,B). Gene trees conflicts were also detected (Figure 6B2,B3), and the final species tree was used in next positive selection analysis.

FIGURE 6

Figure 6. The phylogenetic analysis using concatenate (A) and coalescent method (B) based on single-copy genes. Asterisk (^∗) above the branches in left (A) concatenate tree indicate the maximum support values and posterior probabilities in ML and BI analyses. ^∗Above branches in right (B) coalescent tree are support values of MP-EST and STAR with bootstrapping analyses and local posterior probabilities. Panels (B1–B3) are the coalescent results with the 3 tree topologies and correspondent SCG numbers.

Positive Selection Analysis and Functional Annotation

A total of 335 SCGs exhibited significant positive selection (PSGs) (P-value < 0.05). We found that most of the PSGs were enriched in (GO terms) biological processes including metabolic process, stimulus-response, nitrogen compound metabolic process, membrane, nucleic acid/ion binding, transporter/transferase/hydrolase activity, ion/acid, and types of biological compound transmembrane transport (Figure 7A and Supplementary Table 2). KEGG function analysis indicated that most PSGs participated in pathways that related to signal transcription, transport and catabolism, translation, cell growth and death, carbohydrate metabolism, environmental adaptation, and aging (Figure 7B and Supplementary Table 3). The results of COG enrichment analysis indicated that many PSGs were involved in replication, recombination and repair, transcription, defense mechanisms, translation, ribosomal structure and biogenesis, signal transduction mechanisms, carbohydrate transport and metabolism, secondary metabolites biosynthesis, transport and catabolism, energy production and conversion (Figure 7C). Besides, through adaptive analysis for the 243 clusters of paralogous genes, we found 30 paralogous clusters were obvious enriched, and 3 clusters with p < 0.001, 5 clusters with p < 0.01 and 22 clusters with p < 0.05 in both M1–M2 and M7–M8 test, which were involved in transmembrane transport, salicylic id mediated signaling pathway, gene silencing by RNA, ATPase activity, wax biosynthetic process, ATP binding (Figure 5 and Supplementary Table 1).

FIGURE 7

Figure 7. The function enrichment analyses of PSGs in Urophysa species. (A) GO annotation results; (B) KEGG annotation results; (C) COG annotation results.

By further functional annotation of positively selected genes and BLASTX with A. thaliana, many genes were related to the transmembrane transporter activity of ion/anion/organic anion/carboxylic acid and some were enriched with transmembrane transport (Figure 7 and Table 3), such as At3g20240 (PSG_140) and SKOR (PSG_282). Other genes, included At2g23790 (PSG_16), CAMRLK (PSG_168), are involved in the homeostasis of calcium ions were detected. We also found many genes related to water retention or water deprivation response and contributed to the osmotic stress response [such as HAT22 (PSG_130), CPK13 (PSG_95), DRIP2 (PSG_162) and LEA6 (PSG_305)]. Many of PSGs were found related to genetic information process [e.g., PCMP-H14 (PSG_13), FAS2 (PSG_175), EXO1 (PSG_275)], stress response [ATPK2 (PSG_198), OXI1 (PSG_274), At5g58480 (PSG_323)], photosynthesis and energy metabolism [ndhS (PSG_143), GLYK (PSG_154), PSB28 (PSG_183)]. Detailed annotation information is listed in Supplementary Table 2. Several PSGs with their positively selected sites and three-dimensional protein structures were visualized and are shown in Figures 8, 9. The illustrated PSGs were involved in transmembrane transport, ion homeostasis and water retention (Figures 8A,B, 9A,B), genetic information process (Figures 8D, 9D), stress response (Figures 8E,F, 9E,F), photosynthesis and energy metabolism (Figures 8G,H, 9G,H). Most of the positively selected sites that were detected had BEB posterior probabilities of more than 0.90. The three-dimensional protein structures of PGSs mainly constituted an α-helix and β-sheet. Additionally, some paralogs were also involved in transmembrane transport, calcium ion transport (e.g., ABC family genes), and energy metabolism (e.g., HSP70) (Table 3).

TABLE 3

Table 3. Positively selected genes (PSGs) and clusters of paralogous genes (Mul_gene) related to ion homeostasis and water retention, genetic information process, stress response, photosynthesis, and energy metabolism in Urophysa species.

FIGURE 8

Figure 8. Partial alignment of some positively selected genes. (A–C) Positively selected genes (PSGs) related to transmembrane transport, ion homeostasis and water retention. (D) PSGs related to genetic information process. (E,F) PSGs related to the stress response. (G,H) PSGs related to photosynthesis and energy metabolism. Double red asterisks in the top stand for the amino acids is a positively selected site with a BEB posterior probability of more than 0.95, and one asterisk stands for the sites with a posterior probability larger than 0.90, but lower than 0.95. Numbers at the base of the asterisks stand for the site position in the PSGs. (AC, A. coerulea; AE, A. ecalcarata; SA, S. adoxoides; UH, U. henryi and UR, U. rockii).

FIGURE 9

Figure 9. The secondary protein structures of some positively selected genes. (A–C) Positively selected genes (PSGs) related to transmembrane transport, ion homeostasis and water retention. (D) PSGs related to genetic information process. (E,F) PSGs related to the stress response. (G,H) PSGs related to photosynthesis and energy metabolism. Red color represented the α-helix, the yellow color stand for the β-sheet, and the green color stands for the β-turn.

Discussion

Adaptation is an evolutionary process that allows organisms to survive in their habitat better through natural selection (Lan et al., 2017). Recently, transcriptome data of many non-model species have been obtained and employed for studies of species polyploid, species divergence and differential gene expression (Wang H.F. et al., 2019; Hina et al., 2020; Li et al., 2020), phylogenetic and phylogenomic analysis (Pouchon et al., 2018; Opatova et al., 2019) as well as for detecting natural selection and exploring adaptive evolution in closely related species (Ma et al., 2019; Zeng et al., 2019; Peng et al., 2020; Xie et al., 2020). To date, however, most transcriptome studies have been undertaken on model species and studies of adaptation to the karst environment are not well understood. This is the first study to explore the adaptation mechanisms of species in the karst floristic region through transcriptome data, which will offer excellent opportunities to explore the processes of adaptive evolution.

Stresses or Stimuli Resistance

Stresses and stimuli, such as drought stress, salinity stress, extreme temperatures, chemical toxicity and oxidative stress are serious threats to species and result in the deterioration of the environment (Wang and Altman, 2003). Karst limestones are sedimentary rock outcrops that consist primarily of calcium carbonate and are generally characterized by poor soil development, low soil water content, periodic water deficiency and heat stress (Kang et al., 2014). Under such hostile conditions, plant uptake of nutrients (such as nitrogen and phosphorus) may be restricted, which greatly threaten plant growth, development and survival (Wang et al., 2014). We found that through the functional annotation of shared orthologs (4,588) (Figure 3) and PSGs in two karst species (U. rockii and U. henryi) (Figure 7), genes associated with stress resistance, DNA repair, ion transmembrane transport and homeostasis maintenance were overrepresented (Table 3 and Supplementary Tables 1–3). Among them, a large proportion of PSGs were involved in various stresses or stimulus resistance, such as salt stress, oxidative stress, cold and wounding.

Abiotic stresses are often interconnected and may induce similar cellular damage (Wang and Altman, 2003). For example, oxidative stress, which frequently accompanies high temperature, salinity, or drought stress, may cause denaturation of functional and structural proteins (Smirnoff, 1998). Through a BLASTX search of these PSGs in A. thaliana, PSG_146 was homologous to SODCP (Table 3), which plays an important role in superoxide dismutase activity (Abarca et al., 2001) and it can destroy radicals that are normally produced within the cells and toxic to biological systems (Sunkar et al., 2006). We found that some PSGs were involved in DNA repair, recombination, and RNA modification, which are essential for maintaining genomic stability in organisms (Zhang L. et al., 2013). For example, PSG_175 and PSG_275 were, respectively, homologous to FAS2 and EXO1, which are involved in the repair of DNA double-strand breaks (DBSs) via homologous recombination (Mozgová et al., 2010) and DNA mismatch repair (Schaetzlein et al., 2007), respectively. These PSGs may allow the Urophysa species to overcome the stresses and stimuli in the karst environment.

Homeostasis Related Genes in Karst Adaptation

Because of the relatively higher concentrations of calcium (Ca²⁺) and magnesium (Mg²⁺), higher pH levels and a lower water storage capacity in karst limestone bedrock compared with non-karst soils (Nie et al., 2011; Hao et al., 2015), U. rockii and U. henryi must confront highly alkaline conditions, thin soil layers, and desiccation on porous limestone bedrock. To address this issue, the two species must regulate ion and water transport to maintain osmotic balance. We found more orthologs in the two karst species (4,588) than other three non-karst species (811) (Figure 3), Additionally, our GO enrichment analysis found more orthologs involved in calcium-mediated signaling, nitrogen/phosphorus compound metabolic process, membrane, transport, ion binding, and nutrient reservoir activity (Figure 4 and Supplementary Figure 3). Many PSGs and Mul_genes that were related to the membrane, transmemberane transport, calcium ion transport, transport and ion homeostasis were detected (Figures 5, 7, Table 3, and Supplementary Figure 4). The three-dimensional protein structures of these PGSs were mainly composed of an α-helix and β-sheet, and numerous positively selected sites existed in their amino acid sequences, which might suggest that these genes undergo strong natural selection and contribute to the karst environment adaptation of Urophysa species.

The homeostasis of intracellular ion concentrations is fundamental to the physiology of living cells. Proper regulation of ion flux is necessary for cells to keep the concentrations of toxic ions low and to accumulate essential ions (Zhu, 2001). Positively selected gene SKOR (Table 3), a highly selective outward-rectifying potassium channel, plays an important role in the regulation of ion transmembrane transport (Frédéric et al., 1998). Previous studies have revealed that the outwardly rectifying K⁺ channel, which is coded by the SKOR gene, involves in K⁺ release into the xylem sap toward the shoots, a critical step in the long-distance distribution of K⁺ from roots to the upper parts of the plant (Johansson et al., 2006; Liu et al., 2006). Thus, the potassium channel coded by the SKOR gene may assist in regulating homeostasis in U. rockii and U. henryi in karst habitat. Notably, we found genes PSG_16 and PSG168 that coded the calcium uniporter protein 2 and calmodulin-binding receptor kinase CaMRLK, respectively. Furthermore, we also found some paralogs (e.g., Mul_92 is homologous to gene LCA1) that were associated with calcium ion transport from the cytosol to an endomembrane compartment (Wimmers et al., 1992) (Table 3). The calcium uniporter mediates calcium uptake into mitochondria and regulates the mitochondrial calcium ion homeostasis, which thereby plays a key roles in cellular physiology and regulates cell bioenergetics, cytoplasmic calcium signals and activation of cell death pathways (De Stefani et al., 2011). The Calmodulin-binding receptor kinase CaMRLK is involved in auxin and osmotic stress response (Charpenteau et al., 2004). Moreover, PSGs involved in water deprivation response and hyperosmotic tolerance were detected in this study, such as HAT22, CPK13, DRIP2, LEA6, and fcpA. Among them, CPK13 plays an important role in cellular water deprivation (Saijo et al., 2010), which might function in signal transduction pathways that positively regulate responses to cold, salt and drought stresses. LEA6 is involved in the adaptive response of water deficit for vascular plants (Olvera-Carrillo et al., 2010), which may partially explain why the two Urophysa species were able to inhabit habitats with a high level of aridity. All the calcium-regulator and water-transporter-related genes may play key roles for U. rockii and U. henryi in adapting to the extremely hostile karst environment with severe water deprivation and violent osmosis. Moreover, some paralogs (e.g., Mul_46, Mul_132, and Mul_178) belong to ABC transporter family, and relate to transmembrane transport (Bessire et al., 2011). We also found many membrane related genes from functional annotation of orthologs shared by Urophysa species and PSGs (Supplementary Figure 2), which also suggests that membrane systems play an important role in osmotic response to the karst environment.

Morphological Adaptation of Urophysa Species

Previous studies suggested that the genus Urophysa has a close relationship with Semiaquilegia and Aquilegia and shared a common ancestor (Wang, 1992). Among them, Urophysa species possess obvious cuticular wax on their leaves’ surfaces, which was absent in other three non-karst species. The cuticular wax plays a vitally important role in protecting plant tissue from environmental stresses, limiting non-stomatal water loss, and helping to prevent the germination of pathogenic microbes (Raven and Edwards, 2004; Samuels et al., 2008; Zhang et al., 2010). Furthermore, through scanning electron micrographs of Urophysa species’ leaf epidermis in our previous study (Xie et al., 2017), numerous epidermal hairs were detected. Previous studies suggested that the epidermal hairs can increase the resistance of plants to herbivorous insects (Reymond et al., 2004; Clauss et al., 2006), and others found that epidermal hairs play important roles in the defense and protection of plants by adjusting the content of defense-related compounds (Traw and Bergelson, 2003; Jakoby et al., 2008). Additionally, species of Aquilegia, Semiaquilegia, and Urophysa have rhizomes, however, compared to the non-karst species, Urophysa species possess longer and stout rhizomes and are deeply rooted into the rock cracks, which can help them absorb nutrients and water efficiently (Xie et al., 2017). Moreover, through field investigation, we found that these two species disperse seeds into the karst cracks by bending their pedicels toward the rock cracks in the fruiting period, which allows more seeds to enter the cracks. This reproductive approach might ensure Urophysa species better adapt to the karst habitat. Phenotypic plasticity of Urophysa species associated with karst environment adaptation was also found by Xie et al. (2017) who identified extensive footprints of local adaptation from U. rockii and U. henryi specialized morphologies, including unusual floral organs, such as petaloid sepals that can display various colors in different flower phases, petals with a nectar spur that could attract more pollinators, and a mass of small seeds (thousand seed weight is 0.6684 ± 0.0038 g) (Zhang et al., 2013b). All of these specialized morphologies are likely to contribute to enhance reproductive efficiency, but also be favorable to adaptating to the harsh karst environment.

Conclusion

Our transcriptomic analysis detected many PSGs and clusters of paralogous genes in Urophysa species. Among them, genes related to stress and stimulus resistance, transmembrane transport, calcium ion transport, cellular ion homeostasis, calcium signaling regulation and water retention were detected, which allow better adaptation to the karst environment. Several phenotypic characteristics, such as obvious cuticular wax, numerous epidermal hairs, long and stout rhizomes and unusual flowers, may also contribute to habitat adaptation. Our results showed that Urophysa species may have evolved successful strategies for karst environment adaptation.

Data Availability Statement

The datasets generated for this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

D-FX, R-YC, and X-JH designed the study. D-FX and X-YZ collected the materials and planned the experimental setup. D-FX, R-YC, and XF conducted bioinformatic analyses and drafted the manuscript. D-FX, MP, Y-LL, and C-BW revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the China Postdoctoral Science Foundation (2020M683303), the Fundamental Research Funds for the Central Universities (20826041E4158), the National Natural Science Foundation of China (Grant Nos. 31872647 and 32070221), and the Chinese Ministry of Science and Technology through the National Science and Technology Infrastructure Platform Project (Grant No. 2005DKA21 403-JK).

Conflict of Interest

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.

Acknowledgments

We acknowledge Ting Ren, Jun Wen, and Ling-Jian Gui for their help in materials collection, and we thank Yi-Qi Deng and Jun-Pei Chen for their help in data analysis.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.667988/full#supplementary-material

Supplementary Figure 1 | The assembly completeness test for the transcriptome of each species.

Supplementary Figure 2 | Gene Ontology (GO) annotation for shared orthologs from Urophysa rockii and Urophysa henryi.

Supplementary Figure 3 | Gene Ontology function analysis of the significantly enriched orthologs shared from U. rockii and U. henryi.

Supplementary Figure 4 | The directed acycling graph revealed from the shared orthologs from U. rockii and U. henryi.

Supplementary Figure 5 | Gene Ontology annotation for shared orthologs from non-karst species Semiaquilegia adoxoides, Aquilegia ecalcarata, and Aquilegia coerulea.

Supplementary Table 1 | The function annotation and adaptive significance test for the 243 clusters of paralogous genes.

Supplementary Table 2 | The GO annotation for all 335 PSGs detected from Urophysa species.

Supplementary Table 3 | The KEGG annotation results for 335 PSGs detected from Urophysa species.

Footnotes

1. ^ http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
2. ^ https://busco-data.ezlab.org/v4/data/lineages
3. ^ https://github.com/TransDecoder/TransDecoder/releases
4. ^ http://www.phytozome.net/
5. ^ https://pymol.org/2/

References

Abarca, D., Roldán, M., Martín, M., and Sabater, B. (2001). Arabidopsis thaliana ecotype Cvi shows an increased tolerance to photo-oxidative stress and contains a new chloroplastic copper/zinc superoxide dismutase isoenzyme. J. Exp. Bot. 52, 1417–1425. doi: 10.1093/jexbot/52.360.1417

PubMed Abstract | CrossRef Full Text | Google Scholar

Ai, B., Gao, Y., Zhang, X., Tao, J., Kang, M., and Huang, H. (2015). Comparative transcriptome resources of eleven Primulina species, a group of ‘stone plants’ from a biodiversity hot spot. Mol. Ecol. Resour. 15, 619–632. doi: 10.1111/1755-0998.12333

PubMed Abstract | CrossRef Full Text | Google Scholar

Allouche, O., Kalyuzhny, M., Moreno-Rueda, G., Pizarro, M., and Kadmon, R. (2012). Area–heterogeneity tradeoff and the diversity of ecological communities. Proc. Natl. Acad. Sci. U. S. A. 109, 17495–17500.

Google Scholar

Bessire, M., Borel, S., Fabre, G., Garraca, L., Efremova, N., Yephremov, A., et al. (2011). A member of the PLEIOTROPIC DRUG RESISTANCE family of ATP binding cassette transporters is required for the formation of a functional cuticle in Arabidopsis. Plant Cell. 23, 1958–1970. doi: 10.1105/tpc.111.083121

PubMed Abstract | CrossRef Full Text | Google Scholar

Birkeland, S., Gustafsson, A. L. S., Brysting, A. K., Brochmann, C., and Nowak, M. D. (2020). Multiple genetic trajectories to extreme abiotic stress adaptation in Arctic Brassicaceae. Mol. Biol. Evol. 37, 2052–2068. doi: 10.1093/molbev/msaa068

PubMed Abstract | CrossRef Full Text | Google Scholar

Charpenteau, M., Jaworski, K., Ramirez, B. C., Tretyn, A., Ranjeva, R., and Ranty, B. (2004). A receptor-like kinase from Arabidopsis thaliana is a calmodulin-binding protein. Biochem. J. 379, 841–848. doi: 10.1042/bj20031045

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H., Deng, C., Nie, H., Fan, G., and He, Y. (2017). Transcriptome analyses provide insights into the difference of alkaloids biosynthesis in the Chinese goldthread (Coptis chinensis Franch.) from different biotopes. PeerJ. 2017, 1–16.

Google Scholar

Chen, L. Y., Zhao, S. Y., Wang, Q. F., and Moody, M. L. (2015). Transcriptome sequencing of three Ranunculus species (Ranunculaceae) reveals candidate genes in adaptation from terrestrial to aquatic habitats. Sci. Rep. 2015:10098.

Google Scholar

Clamp, M., Cuff, J., Searle, S. M., and Barton, G. J. (2004). The Jalview Java alignment editor. Bioinformatics 20, 426–427. doi: 10.1093/bioinformatics/btg430

PubMed Abstract | CrossRef Full Text | Google Scholar

Clauss, M. J., Dietel, S., Schubert, G., and Mitchell-Olds, T. (2006). Glucosinolate and trichome defenses in a natural Arabidopsis lyrate population. J. Chem. Ecol. 32, 2351–2373. doi: 10.1007/s10886-006-9150-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Clements, R., Sodhi, N. S., Schilthuizen, M., and Ng, P. K. L. (2006). Limestone Karsts of Southeast Asia: Imperiled Arks of Biodiversity. Bioscience 56, 733–742. doi: 10.1641/0006-3568(2006)56[733:lkosai]2.0.co;2

CrossRef Full Text | Google Scholar

Conesa, A., Götz, S., Garcíagómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610

PubMed Abstract | CrossRef Full Text | Google Scholar

De Stefani, D., Raffaello, A., Teardo, E., Szabò, I., and Rizzuto, R. (2011). A forty-kilodalton protein of the inner membrane is the mitochondrial calcium uniporter. Nature 476, 336–340. doi: 10.1038/nature10230

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, B. G., Zhu, D. Y., Yang, Y. J., Shen, J., Yang, F. L., and Su, Z. Y. (2010). Living situation and protection strategies of endangered Urophysa rockii. Jiangsu J. Agri. Sci. 1, 324–325.

Google Scholar

Duan, X., Zhao, C., Jiang, Y., Zhang, R., Shan, H., and Kong, H. (2020). Parallel evolution of apetalous lineages within the buttercup family (Ranunculaceae): outward expansion of AGAMOUS1, rather than disruption of APETALA3-3. Plant J. 104, 1169–1181. doi: 10.1111/tpj.14985

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunn, C. W., Howison, M., and Zapata, F. (2013). Agalma: an automated phylogenomics workflow. BMC Bioinform. 14:330. doi: 10.1186/1471-2105-14-330

PubMed Abstract | CrossRef Full Text | Google Scholar

Edgar, R. C. (2004). MUSCLE: multiple sequence alignment with improved accuracy and speed. Comput. Syst. Bioinform. Conf. 32, 1792–1797. doi: 10.1093/nar/gkh340

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, C., Wang, J., Wu, L., Kong, H., Yang, L., Feng, C., et al. (2020). The genome of a cave plant, Primulina huaijiensis, provides insights into adaptation to limestone karst habitats. N. Phytol. 227, 1249–1263. doi: 10.1111/nph.16588

PubMed Abstract | CrossRef Full Text | Google Scholar

Frédéric, G., Pilot, G., Lacombe, B., Bouchez, D., and Hervé, S. (1998). Identification and disruption of a plant shaker-like outward channel involved in K+ release into the xylem sap. Cell 94, 647–655. doi: 10.1016/s0092-8674(00)81606-2

CrossRef Full Text | Google Scholar

Gao, Y., Ai, B., Kong, H. H., Kang, M., and Huang, H. W. (2015). Geographical pattern of isolation and diversification in karst habitat islands: a case study in the Primulina eburnea complex. J. Biogeogr. 11, 2131–2144. doi: 10.1111/jbi.12576

CrossRef Full Text | Google Scholar

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883

PubMed Abstract | CrossRef Full Text | Google Scholar

Hao, Z., Kuang, Y., and Kang, M. (2015). Untangling the influence of phylogeny, soil and climate on leaf element concentrations in a biodiversity hotspot. Funct. Ecol. 29, 165–176. doi: 10.1111/1365-2435.12344

CrossRef Full Text | Google Scholar

Hina, F., Yisilam, G., Wang, S., Li, P., and Fu, C. (2020). De novo transcriptome assembly, gene annotation and SSR marker development in the moon seed genus Menispermum (Menispermaceae). Front. Genet. 11:380. doi: 10.3389/fgene.2020.00380

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, D. W., Sherman, B. T., and Lempicki, R. A. (2008). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Q., Cai, Y., and Xing, X. (2008). Rocky desertification, antidesertification, and sustainable development in the karst mountain region of southwest China. AMBIO 37, 390–392. doi: 10.1579/08-s-493.1

CrossRef Full Text | Google Scholar

Jakoby, M. J., Doris, F., Mader, M. T., Brininstool, G., Wischnitzki, E., Platz, N., et al. (2008). Transcriptional profiling of mature Arabidopsis trichomes reveals that NOECK encodes the MIXTA-like transcriptional regulator MYB106. Plant Physiol. 148, 1583–1602. doi: 10.1104/pp.108.126979

PubMed Abstract | CrossRef Full Text | Google Scholar

Johansson, I., Wulfetange, K., Porée, F., Michard, E., Gajdanowicz, P., Lacombe, B., et al. (2006). External K+ modulates the activity of the Arabidopsis potassium channel SKOR via an unusual mechanism. Plant J. 46, 269–281. doi: 10.1111/j.1365-313x.2006.02690.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, M., Tao, J., Wang, J., Ren, C., Qi, Q., Xiang, Q. Y., et al. (2014). Adaptive and nonadaptive genome size evolution in Karst endemic flora of China. New Phytol. 202, 1371–1381. doi: 10.1111/nph.12726

PubMed Abstract | CrossRef Full Text | Google Scholar

Lan, Y., Sun, J., Tian, R., Bartlett, D. H., Li, R., Yue, H. W., et al. (2017). Molecular adaptation in the world’s deepest-living animal: Insights from transcriptome sequencing of the hadal amphipod Hirondellea gigas. Mol. Ecol. 26, 3732–3743. doi: 10.1111/mec.14149

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C. Y. (2006). Classification and Systematics of the Aquilegiinae Tamura, Vol. 7. Beijing: Chinese Academy of Sciences, 3.

Google Scholar

Li, D. M., Zhao, C. Y., Liu, X. R., Liu, X. F., Lin, Y. J., Liu, J. W., et al. (2015). De novo assembly and characterization of the root transcriptome and development of simple sequence repeat markers in Paphiopedilum concolor. Genet. Mol. Res. GMR. 14, 6189–6201. doi: 10.4238/2015.june.9.5

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Milne, R. I., Ru, D., Miao, J., and Mao, K. (2020). Allopatric divergence and hybridization within Cupressus chengiana (Cupressaceae), a threatened conifer in the northern Hengduan Mountains of western china. Mole. Ecol. 29, 1250–1266. doi: 10.1111/mec.15407

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Stoeckert, C. J., and Roos, D. S. (2003). OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 13, 2178–2189. doi: 10.1101/gr.1224503

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., and Godzik, A. (2006). Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22, 1658–1659. doi: 10.1093/bioinformatics/btl158

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Wu, D. D., Boyko, A. R., Wang, G. D., Wu, S. F., Irwin, D. M., et al. (2014). Population variation revealed high-altitude adaptation of Tibetan mastiffs. Mol. Biol. Evol. 31, 1200–1205. doi: 10.1093/molbev/msu070

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, K., Li, L. G., and Luan, S. (2006). Intracellular K+ sensing of SKOR, a shaker-type K+ channel from Arabidopsis. Plant J. 46, 260–268. doi: 10.1111/j.1365-313x.2006.02689.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, L., Yu, L. L., and Edwards, S. V. (2010). A maximum pseudo-likelihood approach for estimating species trees under the coalescent model. BMC Evol. Biol. 10:302. doi: 10.1186/1471-2148-10-302

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y. Q., Xu, Z. Y., Zhao, X., and Zhao, G. (2009). Living environment and cultivation experiments of Urophysa rockii. China Seed Ind. 2, 69–70.

Google Scholar

Liu, L., Yu, L. L., Pearl, D. K., and Edwards, S. V. (2009). Estimating species phylogenies using coalescence times among sequences. Sys. Biol. 58, 468–477. doi: 10.1093/sysbio/syp031

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z. J., Zhang, L. Y., Yan, Z. Z., Ren, Z. J., Han, F. M., and Tan, X. X. (2019). Genomic mechanisms of physiological and morphological adaptations of limestone langurs to karst habitats. Mol. Biol. Evol. 37, 952–968. doi: 10.1093/molbev/msz301

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, X. Q., Wang, C. Y., Yang, H. Y., and Liao, X. R. (2012). Studies on adaptive mechanisms of karst dominant plant species to drought and high calcium stress. Chin. Agric. Bull. 28, 1–5.

Google Scholar

Ma, Y. Z., Wang, J., Hu, Q., Li, J., and Mao, K. S. (2019). Ancient introgression drives adaptation to cooler and drier mountain habitats in a cypress species complex. Commun. Biol. 2:213.

Google Scholar

Minh, B. Q., Schmidt, H. A., Chernomor, O., Schrempf, D., and Lanfear, R. (2020). IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol. 37, 1530–1534. doi: 10.1093/molbev/msaa015

PubMed Abstract | CrossRef Full Text | Google Scholar

Mock, T., Otillar, R. P., Strauss, J., Mcmullan, M., Paajanen, P., Schmutz, J., et al. (2017). Evolutionary genomics of the cold-adapted diatom Fragilariopsis cylindrus. Nature 541, 536–540.

Google Scholar

Mozgová, I., Mokros, P., and Fajkus, J. (2010). Dysfunction of chromatin assembly factor 1 induces shortening of telomeres and loss of 45S rDNA in Arabidopsis thaliana. Plant Cell. 22, 2768–2780. doi: 10.1105/tpc.110.076182

PubMed Abstract | CrossRef Full Text | Google Scholar

Nie, Y. P., Chen, H. S., Wang, K. L., Tan, W., Deng, P. Y., and Yang, J. (2011). Seasonal water use patterns of woody species growing on the continuous dolostone outcrops and nearby thin soils in subtropical China. Plant Soil. 341, 399–412. doi: 10.1007/s11104-010-0653-2

CrossRef Full Text | Google Scholar

Ogata, H., Goto, S., Sato, K., Fujibuchi, W., Bono, H., and Kanehisa, M. (2000). KEGG: Kyoto Encyclopaedia of Genes and Genomes. Nucl. Acids Res. 28, 27–30.

Google Scholar

Olvera-Carrillo, Y., Campos, F., Reyes, J. L., Garciarrubio, A., and Covarrubias, A. A. (2010). Functional analysis of the group 4 late embryogenesis abundant proteins reveals their relevance in the adaptive response during water deficit in Arabidopsis. Plant Physiol. 154, 373–390. doi: 10.1104/pp.110.158964

PubMed Abstract | CrossRef Full Text | Google Scholar

Opatova, V., Hamilton, C. A., Hedin, M., De Oca, L. M., Král, J., and Bond, J. E. (2019). Phylogenetic systematics and evolution of the spider infraorder Mygalomorphae using genomic scale data. Syst. Biol. 69, 671–707. doi: 10.1093/sysbio/syz064

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, C., Ren, J. L., Deng, C., Jiang, D. C., Wang, J. C., Qu, J. Y., et al. (2020). The Genome of shaw’s sea snake (Hydrophis curtus) reveals secondary adaptation to its marine environment. Mol. Biol. Evol. 37, 1744–1760.

Google Scholar

Pouchon, C., Fernández, A., Nassar, J. M., Boyer, F., Aubert, S., Lavergne, S., et al. (2018). Phylogenomic analysis of the explosive adaptive radiation of the Espeletia Complex (Asteraceae) in the Tropical Andes. Syst. Biol. 67, 1041–1060. doi: 10.1093/sysbio/syy022

PubMed Abstract | CrossRef Full Text | Google Scholar

Pratlong, M., Haguenauer, A., Chabrol, O., Klopp, C., Pontarotti, P., and Aurelle, D. (2015). The red coral (Corallium rubrum) transcriptome: a new resource for population genetics and local adaptation studies. Mol. Ecol. Res. 15, 1205–1215. doi: 10.1111/1755-0998.12383

PubMed Abstract | CrossRef Full Text | Google Scholar

Rambaut, A., and Drummond, A. (2009). Tracer v1.5. Available online at: http://beast.bio.ed.ac.uk/Tracer

Google Scholar

Raven, J. A., and Edwards, D. (2004). “Physiological evolution of lower embryophytes: adaptations to the terrestrial environment,” in The Evolution of Plant Physiology: From Whole Plants to Ecosystems, eds A. Hemsley and I. Poole (Amsterdam: Elsevier), 17–41. doi: 10.1016/b978-012339552-8/50003-2

CrossRef Full Text | Google Scholar

Reymond, P., Bodenhausen, N., Van Poecke, R. M. P., Krishnamurthy, V., Dicke, M., and Farmer, E. E. (2004). A conserved transcript pattern in response to a specialist and a generalist herbivore. Plant Cell 16, 3132–3147. doi: 10.1105/tpc.104.026120

PubMed Abstract | CrossRef Full Text | Google Scholar

Ronquist, F., and Huelsenbeck, J. P. (2003). MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574. doi: 10.1093/bioinformatics/btg180

PubMed Abstract | CrossRef Full Text | Google Scholar

Saijo, Y., Hata, S., Kyozuka, J., Shimamoto, K., and Izui, K. (2010). Over-expression of a single Ca2+-dependent protein kinase confers both cold and salt/drought tolerance on rice plants. Plant J. 23, 319–327. doi: 10.1046/j.1365-313x.2000.00787.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Samuels, L., Kunst, L., and Jetter, R. (2008). Sealing plant surfaces: cuticular wax formation by epidermal cells. Annu. Rev. Plant Biol. 59, 683–707. doi: 10.1146/annurev.arplant.59.103006.093219

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaetzlein, S., Kodandaramireddy, N. R., Ju, Z. Y., Lechel, A., Stepczynska, A., and Lilli, D. R. (2007). Exonuclease-1 deletion impairs DNA damage signaling and prolongs lifespan of telomere-dysfunctional mice. Cell 130, 863–877. doi: 10.1016/j.cell.2007.08.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Schnitzler, J., Barraclough, T. G., Boatwright, J. S., Goldblatt, P., Manning, J. C., Powell, M. P., et al. (2011). Causes of plant diversification in the Cape biodiversity hotspot of South Africa. Syst. Biol. 60, 343–357. doi: 10.1093/sysbio/syr006

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwede, T., Kopp, J., Guex, N., and Peitsch, M. C. (2003). SWISS-MODEL: an automated protein homology-modeling server. Nucl. Acids Res. 31, 3381–3385. doi: 10.1093/nar/gkg520

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, X. H., Ma, G. Y., Zou, Q. C., and Tian, Q. D. (2019). Preliminary Analysis of Transcriptome Information and SSR Markers in Young Leaves of Helleborus thibetanus. Mol. Plant Breed. 10, 3297–3304.

Google Scholar

Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V., and Zdobnov, E. M. (2015). BUSCO: Assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210–3212. doi: 10.1093/bioinformatics/btv351

PubMed Abstract | CrossRef Full Text | Google Scholar

Smirnoff, N. (1998). Plant resistance to environmental stress. Curre. Opin. Biotech. 9, 214–219. doi: 10.1016/s0958-1669(98)80118-3

CrossRef Full Text | Google Scholar

Su, W. C., Zhu, W. X., and Xiong, K. N. (2002). Stone desertification and eco-economics improving model in Guizhou Karst Mountain. Carsol. Sin. 21, 19–25.

Google Scholar

Sukumaran, J., and Holder, M. T. (2010). DendroPy: a Python library for phylogenetic computing. Bioinformatics 26, 1569–1571. doi: 10.1093/bioinformatics/btq228

PubMed Abstract | CrossRef Full Text | Google Scholar

Sunkar, R., Kapoor, A., and Zhu, J. K. (2006). Posttranscriptional induction of two Cu/Zn superoxide dismutase genes in Arabidopsis is mediated by downregulation of miR398 and important for oxidative stress tolerance. Plant Cell. 18, 2051–2065. doi: 10.1105/tpc.106.041673

PubMed Abstract | CrossRef Full Text | Google Scholar

Suyama, M., Torrents, D., and Bork, P. (2006). PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments. Nucl. Acids Res. 34, 609–612.

Google Scholar

Talavera, G., and Castresana, J. (2007). Improvement of phylogenies after removing divergent and ambiguously aligned blocks from protein sequence alignments. Syst. Biol. 56, 564–577. doi: 10.1080/10635150701472164

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamura, M. (1995). Phylogeny and classification of the Ranunculaceae. Vienna: Springer.

Google Scholar

Tao, J., Feng, C., Ai, B., and Kang, M. (2016). Adaptive molecular evolution of the two-pore channel 1 gene TPC1 in the karst-adapted genus Primulina (Gesneriaceae). Ann. Bot. 118, 1257–1268. doi: 10.1093/aob/mcw168

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, J., Qi, Q., Kang, M., and Huang, H. (2015). Adaptive Molecular Evolution of PHYE in Primulina, a Karst Cave Plant. PLoS One. 10:e0127821. doi: 10.1371/journal.pone.0127821

PubMed Abstract | CrossRef Full Text | Google Scholar

Tatusov, R. L., Galperin, M. Y., Natale, D. A., and Koonin, E. V. (2000). The COG database: a tool for genome-scale analysis of protein functions and evolution. Nucl. Acids Res. 28, 33–36. doi: 10.1093/nar/28.1.33

PubMed Abstract | CrossRef Full Text | Google Scholar

Traw, M. B., and Bergelson, J. (2003). Interactive effects of jasmonic acid, salicylic acid, and gibberellin on induction of trichomes in Arabidopsis. Plant Physiol. 133, 1367–1375. doi: 10.1104/pp.103.027086

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H. F., Guo, C. C., Ma, H., and Qi, J. (2019). Reply to Zwaenepoel et al.: Meeting the Challenges of Detecting Polyploidy Events from Transcriptomic Data. Mol. Plant. 12, 137–140. doi: 10.1016/j.molp.2018.12.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, K., Shen, Y. J., Yang, Y. Z., Gan, X. N., Liu, G. H., Hu, K., et al. (2019). Morphology and genome of a snailfish from the Mariana Trench provide insights into deep-sea adaptation. Nat. Ecol. Evol. 3, 823–833. doi: 10.1038/s41559-019-0864-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Kang, M., and Huang, H. W. (2014). Long-distance pollen dispersal ensures genetic connectivity of the low-density tree species, Eurycorymbus cavaleriei, in a fragmented karst forest landscape. Conserv. Genet. 15, 1163–1172. doi: 10.1007/s10592-014-0608-x

CrossRef Full Text | Google Scholar

Wang, J. X., and He, X. J. (2011). Preliminary study on Urophysa rockii. I. Literature research and biological characteristics of Urophysa rockii. J. Sichuan Sci. Technol. 32, 69–73.

Google Scholar

Wang, S. J., Liu, Q. M., and Zhang, D. F. (2004). Karst Rocky Desertification in Southwestern China: Geomorphology, Landuse, Impact and Rehabilitation. Land Degrad. Dev. 15, 115–121. doi: 10.1002/ldr.592

CrossRef Full Text | Google Scholar

Wang, W., and Altman, A. (2003). Plant responses to drought, salinity and extreme temperatures: Towards genetic engineering for stress tolerance. Planta 218, 1–14. doi: 10.1007/s00425-003-1105-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W. C. (1992). On some distribution patterns and some migration routes found in the Eastern Asiatic region. J. Syst. Evol. 30, 1–24. doi: 10.4324/9780429500190-1

CrossRef Full Text | Google Scholar

Wimmers, L. E., Ewing, N. N., and Bennett, A. B. (1992). Higher plant Ca(2+)-ATPase: primary structure and regulation of mRNA abundance by salt. Proc. Natl. Acad. Sci. U. S. A. 89, 9205–9209. doi: 10.1073/pnas.89.19.9205

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, C., Mao, X. Z., Huang, J. J., Ding, Y., Wu, J. M., Dong, S., et al. (2011). KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucl. Acids Res. 39, 316–322.

Google Scholar

Xie, D. F., Li, M. L., Tan, J. B., Price, M., Xiao, Q. Y., Yu, Y., et al. (2017). Phylogeography and genetic effects of habitat fragmentation on endemic Urophysa (Ranunculaceae) in Yungui Plateau and adjacent regions. PLoS One. 12:e0186378. doi: 10.1371/journal.pone.0186378

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, D. F., Yu, Y., Wen, J., Huang, J., and He, X. J. (2020). Phylogeny and highland adaptation of Chinese species in Allium section Daghestanica (Amaryllidaceae) revealed by transcriptome sequencing. Mol. Phylogenet. Evol. 146:106737. doi: 10.1016/j.ympev.2020.106737

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, D. F., Zhang, L., Hu, H. Y., Guo, X. L., and He, X. J. (2016). Fragmented habitat drives significant genetic divergence in the Chinese endemic plant, Urophysa henryi (Ranuculaceae). Biochem. Syst. Ecol. 69, 76–82. doi: 10.1016/j.bse.2016.07.021

CrossRef Full Text | Google Scholar

Yang, Z. H. (2007). PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586–1591. doi: 10.1093/molbev/msm088

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Z. H., and Dos, R. M. (2011). Statistical properties of the branch-site test of positive selection. Mol. Biol. Evol. 28, 1217–1228. doi: 10.1093/molbev/msq303

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Z. H., and Nielsen, R. (2002). Codon-substitution models for detecting molecular adaptation at individual sites along specific lineages. Mol. Biol. Evol. 19, 908–917. doi: 10.1093/oxfordjournals.molbev.a004148

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Z. H., Wong, W. S. W., and Nielsen, R. (2005). Bayes Empirical Bayes Inference of amino acid sites under positive selection. Mol. Biol. Evol. 22, 1107–1118. doi: 10.1093/molbev/msi097

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, J., Fang, L., Zheng, H., Zhang, Y., Chen, J., Zhang, Z., et al. (2006). WEGO: a web tool for plotting GO annotations. Nucl. Acids Res. 34, 293–297.

Google Scholar

Yuan, D. X. (1994). Karstology of China. Beijing: Geological Publishing House.

Google Scholar

Zeng, L., Tu, X. L., and Dai, H. (2019). Whole genomes and transcriptomes reveal adaptation and domestication of pistachio. Genome Biol. 20:79.

Google Scholar

Zhai, W., Duan, X. S., Zhang, R., Guo, C. C., Li, L., Xu, G. X., et al. (2019). Chloroplast Genomic Data Provide New and Robust Insights into the Phylogeny and Evolution of the Ranunculaceae. Mol. Phylogenet. Evol. 135, 12–21. doi: 10.1016/j.ympev.2019.02.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J. Y., Broeckling, C. D., Blancaflor, E. B., Sledge, M. K., Sumner, L. W., and Wang, Z. Y. (2010). Overexpression of WXP1, a putative Medicago truncatula AP2 domain-containing transcription factor gene, increases cuticular wax accumulation and enhances drought tolerance in transgenic alfalfa (Medicago sativa). Plant J. 42, 689–707. doi: 10.1111/j.1365-313x.2005.02405.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, T. C., Qiao, Q., Novikova, P. Y., Wang, Q., Yue, J. P., Guan, Y. L., et al. (2019). Genome of Crucihimalaya himalaica a close relative of Arabidopsis, shows ecological adaptation to high altitude. Proc. Natl. Acad. Sci. U. S. A. 116, 7137–7146. doi: 10.1073/pnas.1817580116

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Sun, Y. X., Landis, J. B., Zhang, J. W., Yang, L., Lin, N., et al. (2020). Genomic insights into adaptation to heterogeneous environments for the ancient relictual Circaeaster agrestis (Circaeasteraceae, Ranunculales). New phytol. 228, 285–301. doi: 10.1111/nph.16669

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y. X., Hu, H. Y., and He, X. J. (2013a). Genetic Diversity of Urophysa rockii Ulbrich, an endangered and rare species, detected by ISSR. Acta Bot. Boreali Occident. Sin. 33, 1098–1105.

Google Scholar

Zhang, Y. X., Hu, H. Y., Yang, L. J., Wang, C. B., and He, X. J. (2013b). Seed Dispersal and Germination of an endangered and rare species Urophysa rockii (Ranunculaceae). Acta Bot. Boreali Occident. Sin. 35, 303–309.

Google Scholar

Zhang, L., Yan, H. F., Wu, W., Yu, H., and Ge, X. J. (2013). Comparative transcriptome analysis and marker development of two closely related Primrose species (Primula poissonii and Primula wilsonii). BMC Genomics. 14:329. doi: 10.1186/1471-2164-14-329

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z. F., Li, Y. Y., and Xiao, B. Z. (2016). Comparative transcriptome analysis highlights the crucial roles of photosynthetic system in drought stress adaptation in upland rice. Sci. Rep. 6:19349.

Google Scholar

Zhu, J. K. (2001). Plant salt tolerance. Trends Plant Sci. 6:66.

Google Scholar