Transcriptome sequencing of Crucihimalaya himalaica (Brassicaceae) reveals how Arabidopsis close relative adapt to the Qinghai-Tibet Plateau

The extreme environment of the Qinghai-Tibet Plateau (QTP) provides an ideal natural laboratory for studies on adaptive evolution. Few genome/transcriptome based studies have been conducted on how plants adapt to the environments of QTP compared to numerous studies on vertebrates. Crucihimalaya himalaica is a close relative of Arabidopsis with typical QTP distribution, and is hoped to be a new model system to study speciation and ecological adaptation in extreme environment. In this study, we de novo generated a transcriptome sequence of C. himalaica, with a total of 49,438 unigenes. Compared to five relatives, 10,487 orthogroups were shared by all six species, and 4,286 orthogroups contain putative single copy gene. Further analysis identified 487 extremely significantly positively selected genes (PSGs) in C. himalaica transcriptome. Theses PSGs were enriched in functions related to specific adaptation traits, such as response to radiation, DNA repair, nitrogen metabolism, and stabilization of membrane. These functions are responsible for the adaptation of C. himalaica to the high radiation, soil depletion and low temperature environments on QTP. Our findings indicate that C. himalaica has evolved complex strategies for adapting to the extreme environments on QTP and provide novel insights into genetic mechanisms of highland adaptation in plants.


Results
Transcriptome assembly and annotation. In this study, we generated 132.29 million clean reads and 16.52 Gb of RNA-seq data after quality filtering (Supplementary Table S1). The clean data were submitted to the NCBI Sequence Reads Archive (SRA) database (no. SRR3138110). De novo assembly of these high-quality reads generated 66,084 transcripts and 49,438 unigenes (the longest transcript in one gene) ( Table 1). The unigenes are 896 bp on average with an N50 length of 1,641 bp.
There are 29,760 (60.19%) and 33,667 (68.09%) unigenes identified homologs in Nr and Nt databases on the basis of similarity, respectively. A total of 39,189 (79.26%) unigenes were successfully annotated using at least one database (a total of seven databases, see methods) with a significant match (e-value < 10 −5 ) (Supplementary Table S2). The top species classification hits for C. himalaica in the Nr database are Capsella rubella, A. lyrata, and A. thaliana (Fig. 1A), which all belong to Brassicaceae. The e-values are very significant, which mostly close to zero (32.1%) and 0∼ 1e-100 (16.5%) (Fig. 1B), suggesting that most unigenes of C. himalaica have very similar homologs in above relatives.
Orthogroups identified and positive selection analysis. In this study, a total of 22,572 orthogroups were detected, including 152,769 genes. Among these, 10,487 orthogroups were shared by all six species, and 4,286 orthogroups contained putative one-to-one single copy genes. To further confirm the phylogenetic position of C. himalaica, we generated a phylogenetic tree based on 1,506,379 amino acid sequences of the trimmed and concatenated 4,286 single copy gene alignments from six species (Fig. 2A). The phylogenic tree showed that C. himalaica was closely related to Capsella rubella, and then clustered with the clade of A. lyrata and A. thaliana. The placement of L. stylosa was outside of these two clades with B. rapa as outgroup.
In addition, GO annotation also showed that many PSGs are related to specific adaptation traits, including nitrogen compound metabolic process (150 PSGs), regulation of nitrogen metabolic (73 PSGs), response to radiation (48 PSGs), DNA repair (21 PSGs), and photoperiodism, flowering (13 PSGs) (Fig. 4). In the following paragraphs, we will discuss the important roles of these PSGs in the adaptation of C. himalaica to the various extreme environments in QTP, respectively.

Discussion
Our phylogenetic analysis based on genomic level highly supported that C. himalaica was most closed to Capsella rather than Arabidopsis. The relationship among these species agrees with previous study 20 . Although previous reviews suggested that Crucihimalaya is a specialized group in Himalayas derived from the Arabidopsis genus 18,19 . Our result indicates that it should be derived from the common ancestor of the Arabidopsis lineage and Crucihimalaya lineage. In fact, C. himalaica grouped with C. rubella, but not A. thaliana, is conflict with morphological data, especially the principal character for taxonomy (dehiscent siliques, linear fruits in Crucihimalaya & Arabidopsis vs. dehiscent silicles, obdeltoid fruits in Capsella). Maybe that is why Crucihimalaya was classified into Arabidopsis in previous taxonomic system. This conflict between the morphological and molecular result highlights the need for further studies on the intrinsic mechanism of speciation.
Orthologs are genes that have evolved from a common ancestral gene via speciation 25 . To investigate the selective pressures at the branch level in C. himalaica and its relatives, we estimated the substitution rates for each orthogroups. The median of dN/dS ratio in C. himalaica was significantly larger than that of other relatives, which strongly supported the accelerated evolution in C. himalaica after splitting from its ancestor lineage. Furthermore, the frequency distribution of dN/dS ratios evidently showed that C. himalaica has more genes with elevated dN/ dS ratios than other five species. These accelerated evolution of genes is often driven by positive selection or relaxed selection pressure. Combined with results of branch-site model, which revealed that many functional genes in C. himalaica undergo positive selection in extreme environment of QTP, we speculated that the evolutionary rate increased in C. himalaica are due to positive selection rather than relaxation of select.
There are very high light radiation on QTP, especially the UV-B radiation during the summer (approximately 65 kJ/m 2 ) is among the highest worldwide 26 . The light radiation can influence plant growth and development, such as photoperiodism, flowering, and DNA damage. Our results showed that as many as 43 PSGs in C. himalaica were annotated as response to light stimulus according to GO category, among which, 13 PSGs in the C. himalaica genome were involved in photoperiodism, flowering process (Fig. 4, Supplementary Table S5). Among these, the gene encoding early flowering 6 protein (lysine-specific demethylase, ELF6) is a repressor in the photoperiodic flowering pathway and its loss-of-function mutation causes early flowering 27 . Similarly, another gene encoding sensitive to freezing 6 protein (mediator of RNA polymerase II transcription subunit 16), plays an important role in regulation of the circadian clock and in the control of flowering time 28 . In the short growing season condition of QTP, flowering time is particularly critical and affect both the life cycles and reproductive success of alpine flora 29,30 . C. himalaica starts flowering relatively early (April), and has a very long duration of the flowering period (from April to September). This pattern suggests that C. himalaica has evolved specific reproductive strategies via positively selected more than a dozen genes associated with photoperiodism and flowering process as adaptation to the short growing season in the harsh environments of QTP.
In addition, highly energetic UV radiation causes direct damage to DNA, RNA, and proteins 31 . Our results showed that 25 PSGs were involved in response to DNA damage stimulus (Fig. 4, Supplementary Table S5). This result was consistent with 21 PSGs enriched in DNA repair pathway (Fig. 4, Table 2). As an essential system for correcting UV-induced DNA damage, nucleotide excision repair (NER) play more important roles in UV-resistance for plants 32 . It is worth to mention that several essential genes participating in NER process showed positive selection. Such as gene encoding DNA excision repair protein ERCC-1 (or ultraviolet hypersensitive 7), which has been reported to function as a part of a structure-specific endonuclease that cleaves 5′ to UV photoproducts in DNA 33 . And gene encoding DNA damage-binding protein 2, which forms a complex with DDB1 to recognize damaged DNA and initiation of NER process after exposure to UV light 34 . Another gene, encoding TFIIH protein, is an essential transcription initiation factor that is also pivotal for NER 35 . Besides NER, other DNA repair mechanism, such as base excision repair (BER), mismatch excision repair (MMR) and homologous recombination (HR) also participated in response to light radiation and other external environment stimulation in QTP (Table 2). Accordingly, we found PSGs encoding DNA repair protein XRCC4, DNA mismatch repair protein MUTS2, MUTS protein homolog 2/5 (MSH2/5), MUTL protein homolog 3 (MLH3), DNA repair and recombination protein RAD54, and DNA damage repair/toleration protein DRT111 (Supplementary Table S3). KEGG enrichment showed that these genes covered almost all aspects of the DNA repair mechanism (Supplementary Table S4), suggesting that C. himalaica has evolved an integrated DNA repair mechanism to adapt to the harsh habitat following the uplift of QTP.
Nitrogen compound metabolism is the basic physiological processes, which can generate components of cells (e.g. nucleic acids and proteins) as well as convert sugars and proteins into energy. Moreover, nitrogen compound metabolism is even more important for C. himalaica, because nutrient deficiency, especially nitrogen deficiency is a typical feature of the sandy soil habitats for its living 36,37 . It means nitrogen assimilation under low-nitrogen condition is also a challenge for C. himalaica living in QTP. From our results, it was evident that as many as 150 PSGs involved in nitrogen compound metabolism in C. himalaica, which included GO categories associated with DNA metabolic process (41genes) and regulation of nitrogen compound metabolic process (73 genes) (Fig. 4, Supplementary Table S5). The latter included many transcription factors involved in the nitrogen assimilation process, including the primary assimilation of ammonia to carbon skeletons to biosynthesize amino acids and other organic compounds. Remarkably, three members of DOF (DNA BINDING WITH ONE FINGER) family, a class of zinc finger domain TFs, showed significant positive selection. DOF factors could play an important role in nitrogen regulation, which was supported by the fact that Dof1 over-expressed in Arabidopsis promote the nitrogen assimilation, thus improved plant growth under low-nitrogen conditions 38 . Moreover, a PSG encoding GATA zinc finger protein homolog NTL10, also has been reported activates the expression of nitrogen-catabolic enzymes during conditions of nitrogen limitation 39 . The functions of these PSGs suggested that C. himalaica may have adaptively sped up evolutionary rates of genes involved in nitrogen metabolic regulation to better adapt to the nitrogen limitation in its habitat, which seems to be the driving force behind rapid evolution of these genes.
Low temperatures and rapid changes of temperature are also prevalent features of QTP. Organisms that live in this alpine ecosystem face various growth-related challenges from low temperatures, such as reduced fluidity of lipid membranes 40 . The lipid bilayers transmit external signals to its interior as well as protect the integrality of various organelle membranes in unfavorable environments. Our GO category results indicated that 22 genes related to lipid biosynthesis were under significant positive selection, including 14 PSGs genes involved in phospholipid biosynthetic process and 5 PSGs in glycerolipid biosynthetic process (Supplementary Table  S3). Among these, gene encoding lipid-A-disaccharide synthase with function of condense the lipid A, could play a structural role by stabilizing the plasma membrane lipid bilayer in plants, and involved in signal transduction or plant defense responses 41 . Other PSGs, such as gene encoding Phosphatidyl-N-methylethanolamine N-methyltransferase, was reported to produce the abundant membrane lipid phosphatidylcholine 42 . Previous findings also showed that the degree of unsaturated lipids remains unchanged in C. himalaica under rapid temperature changes 43 , suggesting that the remodeling of membrane lipids might protect membranes against frequent temperature changes. Based on above results, we speculate that these membranes lipids and proteins might play major roles in maintaining membrane integrity under low temperatures and rapid changes of temperature conditions.

Conclusion
Organisms that live in QTP face a variety of external stress from their harsh environments. As such, it is often believed that these organisms have undergone a series of adaptive evolutionary changes. Low oxygen was reported as the most challenge for animals living in high-altitude, thus several genes response to hypoxia showed signature of positive selection in different species [2][3][4][5][6][7][8] . Unlike animals, plants as sessile organisms are constantly exposed to these stresses in QTP. It is conceivable that they have developed more sophisticate mechanisms to protect themselves. Our results of transcriptome sequence of C. himalaica and evolutionary genomic analysis provided evidence for this belief. In this study, we identified 487 significantly PSGs in transcriptome sequences of C. himalaica. Functions of these PSGs elucidate that they potentially possess specific traits of adaptive significance, such as response to radiation, DNA repair, membrane stabilization and organic metabolism. These functions are likely responsible for the adaptation of C. himalaica to the high radiation, low temperature and soil depletion environments on QTP. Due to the RNA-seq technology is difficult to obtain the entire and full length transcripts, further studies are required. Nevertheless, our findings indicate that sophisticated genetic mechanisms have evolved in C. himalaica to survive the harsh conditions in QTP.

Materials and Methods
Sample collection and transcriptome sequencing. The original seedlings of C. himalaica were sampled in August 24, 2013 in Xiangcheng County (alt. 4,104 m, 29°02′ 47″ N, 99°42′ 55″ E) of QTP. These seedlings from the same individual were cultivated in greenhouse in Kunming Institute of Botany. In order to obtain more expression transcripts, two developmental stages (15 days, 30 days) of different organs (leaves, stems, and roots) were sampled and stored at − 80 °C.
High quality total RNA was extracted using the TRIZOL reagent (Sigma Aldrich) following the manufacturer's instructions. A total amount of 3 μ g RNA per developmental stage was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext Ultra ™ RNA Library Prep Kit for Illumina Inc. (NEB, USA) following manufacturer's recommendations. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina Inc.). After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and paired-end reads were generated. The whole step of library construction and Illumina sequencing was performed at Novogene Bioinformatics Technology Co., Ltd (Beijing, China).
De novo assembly and functional annotation. Clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. The sequenced left files (read1 files) from two samples were pooled into one big left.fq file, and right files (read2 files) into one big right.fq file. Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity program (trinityrnaseq_ r20140413) with minimum k-mer coverage set to 2 and all other parameters set by default 44 .
Functional annotations of all assembled unigenes were conducted by searching against the following databases: NCBI non-redundant protein (Nr), NCBI non-redundant nucleotide (Nt), Protein family (Pfam), Clusters of Orthologous Groups of proteins (KOG), Swiss-Prot protein (Swiss-Prot), KEGG Ortholog database (KO), and Gene Ontology database (GO).

Orthologous genes identified and phylogenetic analysis. Based on previous studies in
Brassicaceae 10,20 , we selected genomes of five relatives (A. thaliana, A. lyrata, Capsella rubella, Leavenworthia alabamica and Brassica rapa) and C. himalaica to identify orthologs. Furthermore, to define a set of conserved genes for cross-taxa comparison, we used OrthoMCL software 45 to identify homologous gene clusters (orthogroups) among the six genomes. Genes with lengths less than 50 amino acids were excluded. OrthoMCL was run with  an e-value cut-off of 1e-15 and an inflation parameter of 2.0 due to the close genetic relationship between six relatives. Orthologroups with only single copy genes (one-to-one orthologs) that were shared by all six genomes were retained for further analysis. Each orthogroups was aligned using MUSCLE v3.8.31 46 with default parameters. The poorly aligned regions were further strictly trimmed by using the trimAl v1.4 software 47 with the parameter "-gt 0.8 -st 0.001". Alignments of all orthogroups were concatenated by our python script. Then maximum likelihood (ML) trees were generated using RAxML v7.0.4 48 with PROTCATJTT model, the maximum likelihood criteria.
Positive selection analysis. In the positive selection analysis, only above single copy genes were considered.
To calculate the nonsynonymous (Ka) and synonymous (Ks) substitution rates between pairs of orthogroups, above amino acid alignments were reverse-translated to the corresponding codon-based nucleotide alignments by PAL2NAL 49 . For each alignment, a gene tree was constructed by RAxML 48 using GTR + GAMMA model.
To estimate lineage-specific evolutionary rates for each branch of the six species, the codeml program in the PAML 4 package 24 with the free-ratio model (model = 1) was run on each orthogroups. We conducted the boxplot analysis using the dN/dS ratio derived from free-ratio model results and filtered dS > 3 or dN/dS > 3. Significances of the deviations from the median dN/dS ratio between six species branches were detected using Wilcoxon rank sum test. We also established frequency distribution plots of all dN/dS ratios of six species.
To increase the power of our tests for positive selection, we applied the improved branch-site model 50 implemented in codeml program 24 to estimate the dN/dS substitution rates (ω value). We also deleted all gaps (clean data = 1) from the alignments to lower the effect of ambiguous bases on the inference of positive selection. A foreground branch was specified as the clade of C. himalaica. A significant likelihood ratio test (LRT) was conducted to determine whether positive selection is operating in the foreground branch. In this study, the extremely significant positively selected genes (PSGs) were inferred if the P-value was less than 0.01.
For each PSG in C. himalaica, functional information was inferred based on its ortholog in A. thaliana. Gene Ontology (GO) enrichment analyses of PSGs were conducted using web-based agriGO (bioinfo.cau.edu. cn/agriGO) 51 with singular enrichment analysis (SEA) method and TAIR10 database. The KOBAS software 52 was also used to test the statistical enrichment of PSGs in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways 53 .