Comparative Transcriptome Analysis Reveals Different Mechanisms of Adaptation to Environment among Three Species of Saussurea DC.

Saussurea medusa, Saussurea hypsipeta and Saussurea obvallata are typical alpine snowline plants growing in the Qinghai-Tibet plateau. They are characterized by a specialized morphology. S. medusa and S. hypsipeta have very dense trichomes on whole plant, whereas S. obvallata has transparent bracts covered inflorescences. The different forms reflect their adaptation to cold environments. To investigate the different mechanisms of adaptation of these species to cold temperatures, transcriptome sequencing was performed in three species of Saussurea DC. A total of 116394 137237 and 113879 Unigenes were identified from S. medusa, S. hypsipeta and S. obvallata, respectively. Of these, 55909 (48.03%), 65519 (47.74%) and 51889 (45.56%) Unigenes were matched in public databases. GO analysis identified that most of annotated Unigenes in the three species of plants were related to cellular, metabolic, and single−organism processes, and binding and catalytic activities. The differential expression of 37 genes related to environmental adaptation were discovered by pairwise comparisons. Of these, two candidate genes (Interaptin-like and CSLB3) related to trichome development were identified only in S. medusa and S. hypsipeta, which was consistent with their distinct morphology. Our data can provide a valuable resource for the further studies on the adaptive mechanisms of molecular and functional ecology in Saussurea DC.


Introduction
Saussurea medusa Maxim., Saussurea hypsipeta Diels. and Saussurea obvallata (DC.) Edgew. are perennial herbaceous plants in the family Asteraceae, distributed mainly in the alpine zone of the Qinghai-Tibet Plateau at altitudes from 3800 to 5000 m [1,2]. Many of the high-altitude Saussurea are known for their spectacular growth forms. S. medusa and S. hypsipeta, which are called "woolly" plants, have dense layers of woolly trichomes on their stems, leaves, bracts, and inflorescences [2]. S. obvallata is called called "glasshouse" plant, which encloses their inflorescences in large translucent bracts [3]. These spectacular growth forms are thought to represent an evolutionary response to cold and windy environments [4,5]. However, how are the morphological differences among three the species of Saussurea DC. produced? What are the molecular mechanisms behind these adaptations?
Transcriptomes are the whole RNA transcripts in cells or tissues. They reflect expressed genes at different life stages, tissue types, physiological state and environmental conditions. Transcriptome-based techniques provide useful means of studying gene expression and gene structure, and revealing the molecular mechanism involved in a specific biological process [6]. Recently, this technique has been widely used in the study of molecular mechanisms of plant responses to drought [7], waterlogging [8], low temperature [9], salt [10], high light [11], and irradiation stresses [12]. By comparative transcriptome analysis, new genes are found and molecular mechanisms are revealed in these plants.
Most studies of Saussurea plants growing in the Qinghai-Tibet plateau have focused on their adaptive shapes and anatomies as well as physiological functions [1][2][3]. Little is known about molecular mechanisms behind these plant adaptations [13]. The present work was performed to examine the comparative transcriptome analysis by transcriptome sequencing and gene function annotation on three species of Saussurea plants growing in the Daban Mountain in the northeastern Qinghai-Tibet Plateau. The studies of molecular mechanisms in alpine plant species will contribute to further understand the adaptation of plants to alpine environments since their responses are associated with cold climates; it is probably not useful to know them in the context of a warming climate.

Materials Collection
S. medusa, S. hypsipeta and S. obvallata were collected at flowering stag from Qilian Mountains (37°5′-59′N, 100°55′-102°41′E, 4,200 m a.s.l.) in the northeastern Qinghai-Tibet plateau, China. Growth form of these species is depicted in Fig. 1. This field site features typical plateau continental climate with an annual average (a) temperature of -2°C, (b) rainfall of 482 mm, (c) solar radiation of 6.46 × 109 Jm −2 , and (d) barometric pressure of 684.2 hPa. Three plants were collected from each species. The roots, stems, leaves, and flowers obtained from one individual were pooled into a single sample. These materials were flash-frozen in liquid nitrogen.

RNA Extraction and Transcriptome Sequencing
Trizol ® (Invitrogen-Thermo Fisher Scientific, Carlsbad, CA, USA) was used to extract total RNA on the three species. All RNA was treated with DNase I. The NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA) was used to test RNA purity; the Assay Kit in Qubit ® 2.0 Flurometer (Life Technologies, CA, USA) was used to determine RNA concentrations; the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) was used to assess RNA integrity, and RNA Integrated Number (RIN) values ≥8.
We used NEBNext ® Ultra™ RNA Library Prep Kit from Illumina ® (NEB, USA) to construct the sequences library and the Agilent Bioanalyzer 2100 system to evaluate library quality. We sequenced library preparation to obtain raw reads by Illumina HiSeq TM 2500; clean reads were obtained by removing reads containing adapter, reads containing poly-N and low-quality reads from raw reads. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality. Transcriptome was obtained by splicing clean reads with Trinity [14]. The longest transcript in each gene was taken as the Unigene for subsequent analysis.

Orthology Genes Screening and Analysis
Orthology genes were searched for full-length CDS sequences by OrthoMCL, and one-to-one orthology genes were screened out [15]. Nonsynonymous substitution rate (Ka), synonymous substitution rate (Ks) and Ka/Ks of orthology genes were calculated by PALM [16]. GO gene enrichment analysis was carried out for the orthology genes (Ka/Ks>1) that were positively selected by GOseq [17]. Pathway significant enrichment analysis was carried out on the orthology genes positively selected by the hypergeometric test, with KEGG Pathway as the unit [18].

Analysis of Environmental Adaptation Related Unigenes
By using the classification results of Unigene in the GO database and the functional annotations in the seven public databases, the environmental adaptation-related Unigenes of the three species were analysed by interspecies reduction.

Transcriptome Data Assembly of the Three Species
S. medusa, S. obvallata and S. hypsipeta obtained 150217, 162036 and 192813 transcripts, respectively ( Table 1). The average length of the transcript of S. medusa was 748 bp and N50 was 1292 bp. The average length of the transcript and N50 of S. obvallata and S. hypsipeta were 760 bp, 1270 bp and 698 bp, 1156 bp, respectively. The longest transcript in each gene was taken as Unigene for subsequent analysis. S. medusa obtained 116394 Unigenes with an average length of 623 bp and N50 of 966 bp. S. obvallata showed 113879 Unigenes with an average length of 635 bp and N50 of 1000 bp. S. hypsipeta obtained 13723 Unigenes with an average length of 581 bp and N50 of 857 bp. The results showed that the sequencing quality was sufficient for future analysis needs.

Unigene Annotation in the Public Database
S. medusa, S. hypaipeta and S. obvalata had 7620, 7245 and 9366 Unigenes, respectively, that can be matched in all databases ( Table 2).

GO Functional Classification of Unigene
The Unigenes of S. medusa, S. obvallata and S. hypsipeta were divided into 55, 56 and 56 types of function respectively, and the 55 types of function were the same (Fig. 2). The Unigene of S. obvallata and S. hypsipeta found a gene related to cell aggregation. Cell agglutination is associated to the resistance of organisms to environmental UV stress [19].

KOG Functional Classification of Unigene
The KOG annotation results of Unigene of S. medusa, S. hypsipeta and S. obvallata were similar. There were 21053 Unigenes of S. medusa, 17117 Unigenes of S. hypsipeta and 24003 Unigenes of S. obvallata annotated on 25 KOG classifications. Firstly, the general function prediction only accounted for most of the Unigenes, which were 16.41%, 18.05% and 16.46% in S. medusa, S. hypsipeta and S. obvallata, respectively. Next, there were posttranslational modification, proteins turnover and chaperones. The proportion of S. medusa, S. hypsipeta and S. obvallata in this category was 14.32%, 14.41% and 11.64%, respectively. Signal transduction mechanisms were 8.74%, 8.54% and 7.65%, in S. medusa, S. hypsipeta and S. obvallata, respectively (Fig. 3). They accounted for the least Unigene in the classification of Cell motility, Extracellular structures and Nuclear structures.

KEGG Functional Classification of Unigene
The Unigene annotated by the three species belonged to Cellular Processes, Environmental Information Processing, Genetic Information Processing, Metabolism and Organismal Systems, including 19 secondary pathways (Fig. 4). Thereinto, 7 types of pathway Unigene accounted for a large proportion (>50%), and they included Cellular Processes, Folding and sorting and degradation, Translation, Amino acid metabolism, Carbohydrate metabolism, Energy metabolism, and Overview, respectively. The KEGG classification results indicated that the Unigene functions were mainly related to protein translation and material metabolism.

Orthology Genes Screening of the Three Species
10811 sets of orthologous genes were obtained through OrthoMCL in the three species. There were 463 sets of orthologous genes which Ka/Ks>1, indicating that these genes were rapidly divergent (Fig. 5).
The GO functional enrichment classification was carried out for the orthologous genes of positive selection in the three species. We did not find functional enrichment genes, but 78 genes were related to plant resistance. The resistance genes included peroxidase, ubiquitin, zinc-finger protein and  acyl-CoA-binding protein, which play an important role in plant environmental stress. The KEGG gene enrichment analysis revealed that Ribosome had the largest number of genes enriched in this pathway, with a total of 14 genes, revealing that the genes related to Ribosome synthesis in these three species are subject to strong positive selection (Fig. 6).

Analysis of Environmental Adaptation Related Unigenes
Based on the classification of the three species Unigenes in GO database, and analysis of Unigene about environmental adaptation, the results showed that the genes related to environmental adaptation in the three species were mainly molecular chaperone, ubiquitin, calmodulin, enzyme and ribosomes (Table 3). HSC82, CLPB and HSP17.9A can enhance the effect of plants on high temperature and salt stresses [20][21][22]. Dna J2 and osigba0134h18.3 were annotated by GO, which are involved in the adaptation of plants to high  temperature stress. PER50 is one of the enzymes with which plants respond to oxidative stress [23], while S2P corresponds to salt stress [24]. ATPK2 is involved in plant adaptation to cold and high salt stresses [25], while CRK29 and CIPK7 are in connection with plant immune responses and adaptation to cold stress [26,27]. FEN1 is one of the key enzymes in DNA replication and repair in eukaryotes [28]. RIN1 annotated by GO is involved in plant defensive responses to fungi. Malate oxidoreductase can regulate plant growth and respiration [29]. Gene expression is regulated by SAHH1 through gene methylation modification [30]. GABA-TP1 participates in plant adaptation to temperature stress [31]. SAMT is related to plants' defense responses to biologic stimulation [32,33]. Calm3 mediates the regulation of enzymes and ion channels through calcium ions [34]. IAA13 regulates growth and development of plants through ARFs [35]. H0901F07.20 annotated by GO involves in the binding of lipoyl coenzyme A. ATG8 is related to the formation of autophagosomes in plants [36]. NIP3 is correlated with the transport of the heavy metal arsenic in plants [37]. OJ1754_E06.16 is associated with the processes of protein secretion and signal transmission. FAF takes part in plant ABA activation signaling pathway and phosphorylation regulation. MutS participates in DNA mismatch repair. GF14A is connected with plant adaptation to drought stress [38]. NIMIN-2 is related to the acquisition of systemic resistance in plants [39], while ATL4M is related to abiotic stress responses [40]. Interaptin-like links the intracellular system to the cytoskeleton, which regulates the morphogenesis of the multicellular trichomes [41,42]. CSLB3 is related to the formation of non-cellulosic polysaccharide skeletons in plant cell walls. GBF is one of the transcription factors in which cells respond to signals [43]. AFG2 is associated with the maturation of 60S ribosomal subunit [44]. PGPS/NH15 annotated by GO participates in the regulation of oxidationreduction enzyme activities.

Discussion
Our sequencing results showed that Unigene N50 of the three species were 857~1000 bp, the base quality value Q30 was over 92%, and the GC content was over 43%. The amount of sequencing data and the number of Unigene spliced were relatively large, which showed that the sequencing quality was good, and the sequence information generated was sufficient and effective [45]. The results of Go, KOG and KEGG of three species Unigenes were very similar, indicating that their genetic relationship was very similar.
The temperature in the collecting area was between 0-30°C, UV radiation up to 7000 μW·cm −2 , and the oxygen partial pressure was 13.25 kPa for the three species. Stresses of low oxygen, strong radiation and large temperature differences made that makes the plants growing in this area formed special stress resistance mechanisms in the long-term to achieve collaborative resilience, such as trichomes, ABA signaling pathway, molecular chaperones and ubiquitin pathway. In our study, there was a gene expression related to cell aggregation in both S. hypsipeta and S. obvallata. Some research shows that Gloeocapsa sp. cells with low metabolic activity aggregate to form a lamellar biofilm to protect organisms against ultraviolet radiation [19]. Then, we compared leaf structure among the three species; S. medusa had the longest, thickest and white trichomes on the leaf epidermis. This structure can effectively resist and reflect ultraviolet radiation. S. hypsipeta and S. obvallata may resist the ultraviolet radiation stress through cell aggregation, but the specific mechanism needs to be verified in future experiments. Interaptin protein was the binding protein between the components of the inner membrane and the actin of the cytoskeleton, Microtubules and actin promote the development of trichomes [41,42]. CSLB3 participates in the formation of the noncellulose polysaccharide skeleton in plant cell wall. In this study, we found that the Interaptin-like and CSLB3 genes were both expressed in S. hypsipeta and S. medusa. It is speculated that the expression of these two genes can promote the occurrence of trichomes and improve the resistance of these two plant species to radiation and temperature stresses. This finding was also consistent with the structural characteristics of the three species. S. medusa and S. hypsipeta are densely covered with trichomes. The ground tissues of S. obvallata showed no trichomes, and the inflorescences were enclosed in translucent bracts. At the same time, we found that the three plant species expressed different molecular chaperones, ubiquitin, protease and other genes to allow the adaptation of plants to the environmental stresses, and their mechanisms need further study.
Funding Statement: This work was financially supported by the National Natural Science Foundation of China (31960222, 31360095).

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.