Transcriptomic Analysis Reveals the Mechanism of Picea crassifolia Survival for Alpine Treeline Condition

: The physiological mechanisms driving treeline formation succession captured the attention of ecologists many years ago, yet they are still not fully understood. In this study, physiological parameters (soluble sugars, starch, and nitrogen) were investigated in combination with transcriptomic analysis in the treeline tree species Picea crassifolia . The study was conducted in the middle of Qilian Mountain Reserves, Gansu Province, China, within the elevation range of 2500 ‒ 3300 m. The results showed that the concentrations of non-structural carbohydrates decreased with increasing elevation in the current-year needles and current-year branches, as well as in the coarse and fine roots. RNA-Seq demonstrated that 483 genes were upregulated and 681 were downregulated in the comparison of 2900 and 2500 m (2900 vs. 2500), 770 were upregulated and 1006 were downregulated in 3300 vs. 2500, and 282 were upregulated and 295 were downregulated in 3300 vs. 2900. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that the differentially expressed genes were highly enriched in photosynthesis-related processes, carbon fixation and metabolism, and nitrogen metabolism. Furthermore, almost all photosynthesis-related genes were downregulated, whereas many genes involved in cuticle lipids and flavonoid biosynthesis were upregulated, contributing to the survival of P. crassifolia under the treeline condition. Thus, our study provided not only molecular evidence for carbon limitation hypothesis in treeline formation, but also a better understanding of the molecular mechanisms of treeline tree survival under adverse conditions. one SDHA (MA_10427069g0010) encoding succinate dehydrogenase (ubiquinone) flavoprotein subunit, and two GAD s (MA_161013g0010, MA_89828g0010) encoding glutamate decarboxylase. These results suggested that the expression of genes involved in photosynthetic C assimilation and energy metabolism was significantly altered by the increasing elevations. two GLYM (serine hydroxymethyltransferase) for serine. The other five DEGs were upregulated, including two ARGI ( arginase ) and SPD (spermidine synthase) for arginine, GAD (glutamate decarboxylase) for γ -aminobutyric acid, and SDC (serine decarboxylase) for serine. This also indicated that high altitude conditions mostly had a negative effect on amino acid metabolism.


Introduction
The treeline, defined here as the highest altitude where groups of upright trees reach at least 3 m of height, is one of the most recognizable bioclimatic boundaries and it has been suggested as a critical indicator of trees' response to climate warming [1,2].As one of the important ecological boundaries, the alpine treeline is globally determined by the soil temperature of 6.5 ± 0.8 • C at 10 cm of depth during the vegetation season [3].Although the decrease in temperature with increasing elevations has been recognized as a major reason for treeline formation, and the mechanism by which low temperature affects the growth of treeline trees is unclear.The concordance of this treeline isotherm suggests a common low-temperature effect on plant metabolism, either related to carbon (C) gain (source activity; photosynthesis) or sink activity (tissues formation; meristems) [4].
The physiological mechanisms driving treeline formation captured the attention of ecologists long ago, and various environmental and biophysiological hypotheses have been formulated [5,6], yet they are not fully understood.The environmental hypothesis stated that the formation of treeline was due to the following factors: climatic stress, mechanical damage, short growing season, and temperature-limited reproduction [7][8][9][10].This hypothesis can only be used to interpret local treeline phenomena, which are related to various local site conditions [11].However, the biophysiological hypothesis focuses on growth limitation and C supply limitation [1,2,12].The growth limitation hypothesis (GLH) considers that colder temperatures at higher elevations limit C sink activity [5], supported by the results of trees in 13 alpine tree line regions [13] and Larix potaninii in the eastern Himalayas treeline [14].This is contrary to the C limitation hypothesis (CLH), which states that reduced C gain with higher elevation (lower temperature) results in treeline formation [5], which is further supported by the finding that the elevational decrease in non-structural carbohydrate (NSC) concentrations in Nothofagus pumilio trees at treeline elevation [15].However, in the birch (Betula ermanii) near the treeline on the Changbai Mountain, Northeast China, the decrease in NSCs with elevation in main storage organs supports the CLH, whereas the increasing or static trends in new developing organs support the GLH [16].Moreover, limitation of nutrients, including P and N, is also found to contribute to the formation of the Arctic treeline, particularly in areas with cold winters [17].So far, there is no conclusive evidence about the formation of altitudinal treelines [11][12][13][14][15][16][17][18].Hence, these conflicting findings prompted us to further explore a new way to explain the physiological mechanism that drives high-altitude treeline formation.
Changes in temperature and water with the increasing altitudes have been shown to have significant effects on the physiological adaptations of alpine trees [6,19], which is a consequence of plant responses at the transcriptional level.Although physiological changes in trees of treelines have been reported [1,5,10,20,21], little information is available about the mechanisms regulating treeline formation at the transcriptomic level, compared with those other vegetation zones.In recent years, RNA sequencing (RNA-Seq) has enabled to detect differentially expressed genes (DEGs) because of its dynamic range of expression levels under abiotic stress, such as salinity [22] and drought [23].
In the present study, in order to verify the GLH or CLH at the molecular level, we conducted a transcriptome analysis of Picea crassifolia growing at different elevations, and combined it with the analysis of physiological changes in trees.Consequently, we examined the expression responses of candidate genes related to carbohydrate fixation and metabolism.In order to assess the adaptive potential of trees in the treeline environment, it is necessary to identify genes that are related to environmental responses of trees forming treelines.

Site Description
This study was carried out in the middle of the Qilian Mountain Reserves in the treeline ecotone of Pailugou Basin, Gansu Province (100 • 17' E, 38 • 24' N), with the elevation ranging from 2500 to 3400 m.The vegetation season in the treeline usually lasts from May to mid-October.The annual mean air temperature is around −0.6 • C to −2 • C, and the mean annual precipitation is 430 mm.The soil in the research site is mountain gray cinnamon soil, and the vegetation in the alpine treeline ecotone is dominated by P. crassifolia (ranging from 2500 to 3300 m), and treelike multi-stem shrubs such as Potentilla fruticosa, Caragana jubata, and Salix gilashanica (ranging from 2500 to 3400 m).

Sample Collection
Different tissues, including current-year needles, current-year branches, coarse roots, and fine roots, were collected during the early growing season (in June) at three elevations across the treeline ecotones on north-facing slopes: closed forest stands at lower elevations (2500 m), middle forest stands between closed forest and treeline (2900 m), and treeline stands in the uppermost limit of P. crassifolia trees (3300 m).Sampling was performed according to the method described by Li et al. [10] and Ma et al. [24] with some modifications.Briefly, nine healthy mature trees, consisting of three individuals each repeat, were selected for repeated sampling at each elevation.The current-year needles and current-year branches were collected from non-shaded leading branches on the upslope canopy side.Simultaneously, coarse roots (0.5-1.0 cm in diameter) and fine roots (<0.5 cm in diameter) were collected.One part of the current-year needle samples was immediately stored in liquid N for RNA extraction, whereas the other part, intended for physiological measurements, was immediately stored in a cool box.After short time treatment in a microwave oven (40 s at 600 W), the sample was dried at 70 • C until it reached a constant weight after two consecutive days.

Non-Structural Carbohydrate Analysis
Soluble sugars: Approximately 0.5 g of dry powdered sample was added into a triangular bottle along with 46 mL of distilled water.After boiling the sample for 2 h in a pressure cooker and cooling it, we set the capacity of sample solution to 50 mL, and filtered the solution.The filtrate was analyzed with the Waters 2695 High Performance Liquid Chromatograph (Waters-Millipore, Milford, MA, USA) using a Sugar-Pak I column.The liquid phase was distilled water, flow rate was 0.6 mL•min −1 , column temperature was 70 • C; a parallax detector was used.All separated components were summed as soluble sugars [25,26].
Starch: Approximately 0.1 g of dry powdered sample was added into a plug tube with 10 mL of distilled water, and 1 mL of 2:1 hydrochloric acid, and incubated at 100 • C for 8 h.After cooling, the pH of the sample was adjusted to neutral with 40% NaOH, and the capacity of sample solution was set to 15 mL.After filtration, the filtrate was used for starch analysis.Starch content was analyzed using the same method as the determination of soluble sugars.The total concentration of NSC was defined as the sum of soluble sugar and starch concentrations [26,27].
Total N: Total N content (% d.m.) was determined in finely ground oven-dried samples with the Automatic Kjeldahl Analyzer (UK 152, Italy), using CuSO 4 , K 2 SO 4 , and H 2 SO 4 for digestion, and NH 3 was determined using the indophenol-blue colorimetric method according to Li et al. [10].

RNA-Seq Library Construction, Sequencing, and Mapping
Nine libraries for RNA-Seq, including three biological replicates for 2500, 2900, and 3300 m, were constructed according to the manuals provided by Illumina and sequenced on the HiSeq2000 platform as described by Zheng et al. [22].Briefly, the total RNA was isolated from the current-year needles using the PLANT easy extraction protocol (BLKW, Beijing, China), and mRNAs were purified, fragmented, and used for complementary DNA (cDNA) synthesis with random hexamer primers and Superscript II reverse transcriptase (Invitrogen, Carlsbad, CA, USA).After the purification of the PCR products by AMPure XP system, (Beckman Coulter, Brea, CA, USA) and assessing the library quality on the Agilent Bioanalyzer 2100 system (Agilent, Palo Alto, CA, USA), paired-end cDNA libraries were sequenced using the Illumina-Solexa sequencing platform.
The poplar reference genome and annotation files were downloaded directly from the National Center for Biotechnology Information (NCBI) website (ftp://ftp.ncbi.nih.gov/genomes/PLANTS/Populus_trichocarpa/, accessed on 13 January 2020).An index of the reference genome was built using Bowtie v2.0.6, and single-end clean reads were aligned with the reference genome using TopHat v2.0.9 [28].Read counts were calculated using HTSeq (www.huber.embl.de/users/anders/HTSeq,accessed on 13 January 2020) and normalized to the expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced (FPKM) [29] in order to obtain the relative expression levels.

Identifying and Clustering of DEGs
The R packages of DESeq [30] were used for pairwise gene expression comparisons between two different elevations.The significance of DEGs was identified with a 'p-value ≤ 0.05' and '|log2Ratio| ≥ 1.0' with at least one FPKM ≥ 1 within two elevations.Hierarchical clustering was performed using pheatmap package of R statistical environment in order to cluster the identified DEGs on the base of the FPKM values.
All DEGs were mapped according to the terms from Gene Ontology (GO), using the bioconductor package GoSeq [31], and Kyoto Encyclopedia of Genes and Genomes (KEGG), according to the methods of Kanehisa [32].Significantly enriched terms were identified in comparison with the genome background.
2.6.Quantitative Real-Time PCR (qRT-PCR) cDNA was synthesized from 1.0 µg of total RNA and digested with DNase I using the PrimeScript RT reagent kit for RT-PCR (Takara, Dalian, China), according to the manufacturer's instructions.The specific primers for poplar (Table S1) were synthesized by the Sangon Biotechnology Company (Shanghai, China).The qRT-PCR was performed using the SYBR Premix ExTaq II kit (Takara) on Roche LightCycler 480 (Roche, Penzberg, Germany).To quantify the relative expression levels of the target genes, the poplar reference gene Actin (GenBank accession no.AAF03692) was used as the internal control [33].Fold changes were calculated by the relative expression of treatments compared with that of the control.Three independent biological replicates were performed for each selected gene (at least two technical repeats per biological replicate).

Statistical Analysis
Quantitative data of NSC concentration, N concentration, and transcription analysis were analyzed by one-way analysis of variance (ANOVA) according to the general linear model.Post-hoc comparison was performed by Tukey's honestly significant difference (HSD) test.Differences were accepted as statistically significant when p < 0.05.All statistical analyses were performed using SPASS 16.0.

Changes of NSC and Total N content with Altitudes
With increasing elevation, soluble sugar content in the current-year needles did not show an obvious change (p > 0.05), although there was a decreasing trend at 3300 m (Figure 1A).However, in the current-year branches, there was an apparent increase (15.5%) in the soluble sugar content at 2900 m, but a significant decrease (11.3%) was recorded at 3300 m (Figure 1D).The starch content showed a sharp decreasing trend (70.0%-85.4%) in both needles and branches at 2900 m (Figure 1B,E), but the starch content at 3300 m was significantly higher than that at 2900 m, indicating that the conversion of starch into soluble sugars could be inhibited at 3300 m.These changes resulted in a significant decrease over 17.9%-21.9%in total NSC in the current-year needles at 2900 and 3300 m, as well as in the current-year branches at 3300 m (18.9%) with the increasing elevation.These results suggested that the conditions at higher elevation limited the supply of C in the current-year needles and branches.Similarly, with increasing elevations, soluble carbohydrates were also negatively affected in both coarse roots and fine roots (Figure 2), which resulted in a decrease in the total NSC content, especially at 3300 m (Figure 2C,F).The total N content increased evidently by 18.8%-38.5% at 2900 m compared with that at 2500 m in both needles and branches, whereas at 3300 m, the total N content showed a decrease, but it was still higher than that at 2500 m (Figure 3A,B).However, the total N content significantly decreased by 12.2%-15.2% in the coarse roots (Figure 3C), whereas it obviously increased by 11.3%-36.8% in the fine roots (Figure 3D).These results demonstrated that N was not a main limiting factor for the growth of P. crassifolia under the treeline condition.The total N content increased evidently by 18.8%-38.5% at 2900 m compared with that at 2500 m in both needles and branches, whereas at 3300 m, the total N content showed a decrease, but it was still higher than that at 2500 m (Figure 3A,B).However, the total N content significantly decreased by 12.2%-15.2% in the coarse roots (Figure 3C), whereas it obviously increased by 11.3%-36.8% in the fine roots (Figure 3D).These results demonstrated that N was not a main limiting factor for the growth of P. crassifolia under the treeline condition.The total N content increased evidently by 18.8%-38.5% at 2900 m compared with that at 2500 m in both needles and branches, whereas at 3300 m, the total N content showed a decrease, but it was still higher than that at 2500 m (Figure 3A,B).However, the total N content significantly decreased by 12.2%-15.2% in the coarse roots (Figure 3C), whereas it obviously increased by 11.3%-36.8% in the fine roots (Figure 3D).These results demonstrated that N was not a main limiting factor for the growth of P. crassifolia under the treeline condition.

Identification of Differentially Expressed Genes in P. crassifolia at Different Elevations by RNA-Seq
A summary of raw data is shown in Table S2, and the quantitative expression of genes (estimated using the RPKM values) is shown in Table S3.A high correlation of R 2 > 0.91 between biological replicates was observed for three treatments at different elevations (Figure S1).

Identification of Differentially Expressed Genes in P. crassifolia at Different Elevations by RNA-Seq
A summary of raw data is shown in Table S2, and the quantitative expression of genes (estimated using the RPKM values) is shown in Table S3.A high correlation of R 2 > 0.91 between biological replicates was observed for three treatments at different elevations (Figure S1).
In order to determine the gene expression changes in P. crassifolia resulting from the variations in elevation, differentially expressed genes (DEGs) were identified in response to the increasing elevation (2500, 2900, and 3300 m).The |log2Ratio|≥ 1 was used as the threshold to screen the DEGs and to obtain more information on genes potentially involved in the response to the elevation variations.We found 483 upregulated genes and 681 downregulated genes by comparison of trees from 2900 m with trees from 2500 m (2900 vs. 2500), 770 upregulated genes and 1006 downregulated genes in 3300 vs. 2500, and 282 upregulated genes and 295 downregulated genes in 3300 vs. 2900 (Figure 4A; Table S4).Moreover, only 13 DEGs were common to 3 pairwise comparisons of different elevations, but 1146 DEGs were specific to the elevation of 3300 m (Figure 4B).Furthermore, gene expression patterns of 12 randomly selected DEGs were verified by qRT-PCR, which demonstrated similar expression trends with RNA-Seq results, indicating the reliability of DEG analysis (Figure 5).This indicated the presence of a large proportion of DEGs in P. crassifolia growing in the treeline environment, and these DEGs could aid the species to cope with the harsh environmental conditions.In order to determine the gene expression changes in P. crassifolia resulting from the variations in elevation, differentially expressed genes (DEGs) were identified in response to the increasing elevation (2500, 2900, and 3300 m).The |log2Ratio|≥ 1 was used as the threshold to screen the DEGs and to obtain more information on genes potentially involved in the response to the elevation variations.We found 483 upregulated genes and 681 downregulated genes by comparison of trees from 2900 m with trees from 2500 m (2900 vs. 2500), 770 upregulated genes and 1006 downregulated genes in 3300 vs. 2500, and 282 upregulated genes and 295 downregulated genes in 3300 vs. 2900 (Figure 4A; Table S4).Moreover, only 13 DEGs were common to 3 pairwise comparisons of different elevations, but 1146 DEGs were specific to the elevation of 3300 m (Figure 4B).Furthermore, gene expression patterns of 12 randomly selected DEGs were verified by qRT-PCR, which demonstrated similar expression trends with RNA-Seq results, indicating the reliability of DEG analysis (Figure 5).This indicated the presence of a large proportion of DEGs in P. crassifolia growing in the treeline environment, and these DEGs could aid the species to cope with the harsh environmental conditions.These DEGs may be involved in different functions.The analysis of GO enrichment and KEGG pathway was applied to gain an insight into their functions.DEGs in 2900 vs. 2500, 3300 vs. 2500, and 3300 vs. 2900 were annotated in GO, and categorized into three major categories (biological process, cellular component, and molecular function) (Figure 6; Table S5).GO annotation suggested that biological processes and molecular functions, including "primary metabolism", "organic substance metabolism", "biological regulation", "organic nitrogen compound biosynthesis", "cellular nitrogen compound biosynthesis", "cellular aromatic compound metabolism", "carbohydrate metabolism", and "lipid metabolism" were related to C and N metabolism at different elevations.In particular, biological processes of C and N metabolism were significantly enriched in 3300 vs. 2500.These DEGs may be involved in different functions.The analysis of GO enrichment and KEGG pathway was applied to gain an insight into their functions.DEGs in 2900 vs. 2500, 3300 vs. 2500, and 3300 vs. 2900 were annotated in GO, and categorized into three major categories (biological process, cellular component, and molecular function) (Figure 6; Table S5).GO annotation suggested that biological processes and molecular functions, including "primary metabolism", "organic substance metabolism", "biological regulation", "organic nitrogen compound biosynthesis", "cellular nitrogen compound biosynthesis", "cellular aromatic compound metabolism", "carbohydrate metabolism", and "lipid metabolism" were related to C and N metabolism at different elevations.In particular, biological processes of C and N metabolism were significantly enriched in 3300 vs. 2500.In the KEGG pathway analysis, 83, 96, and 56 pathways were categorized into 14 major KEGG pathways from the pairwise comparisons between 2900 vs. 2500, 3300 vs. 2500, and 3300 vs. 2900, respectively.The pathways of metabolism, mainly including secondary metabolites, carbohydrates, energy, lipid, amino acids, cofactors, and vitamins, were mostly influenced by the variations in In the KEGG pathway analysis, 83, 96, and 56 pathways were categorized into 14 major KEGG pathways from the pairwise comparisons between 2900 vs. 2500, 3300 vs. 2500, and 3300 vs. 2900, respectively.The pathways of metabolism, mainly including secondary metabolites, carbohydrates, energy, lipid, amino acids, cofactors, and vitamins, were mostly influenced by the variations in elevation, especially for 3300 vs. 2500, rather than being influenced by signal transduction, environmental adaptation, and transport (Figure 7; Table S6).The response of these pathways to increasing elevations could be helpful for the survival of P. crassifolia.Among these responses, the metabolism of carbohydrates and energy was enriched in 23 pathways (Table S6).Compared with those at 2500 m, 39 DEGs were upregulated and 40 were downregulated at 2900 m, whereas 74 were upregulated and 44 were downregulated at 3300 m.This indicated that most of the genes were negatively regulated around the elevation of the treeline compared with the lower elevations.We found that all 40 genes enriched in photosynthesis, antenna proteins of photosynthesis, and C fixation in photosynthetic organisms were downregulated at 3300 m.Similarly, DEGs involved in glyoxylate and dicarboxylate metabolism ( 14), pentose phosphate pathway (9), and fructose and mannose metabolism (9) were all downregulated at 3300 m.Moreover, 26 of the 30 DEGs in the C metabolism were also downregulated at 3300 m and 19 DEGs enriched in cutin, suberine, and wax biosynthesis were upregulated at 3300 m.These results demonstrated that the increasing elevations inhibited C fixation and altered the direction of C fluxes.

Analysis of Gene Expression Pattern
To further explore the gene expression patterns in response to the increasing altitudes, we applied K-means clustering to the identified DEGs.In this analysis, 1466 DEGs with a knownfunction were grouped into six categories (referred to as Sub1, Sub2, Sub3, Sub4, Sub5, and Sub6), comprising 113, 616, 165, 220, 194, and 158 genes, respectively (Figure 8; Table S7).The DEGs in the Sub1 and Sub2 groups were characterized by constant upregulation with the increasing elevations (Table S7).The GO and KEGG analysis showed that these genes were mainly enriched in secondary metabolite metabolism (such as biosynthesis of secondary metabolites; flavonoid biosynthesis; and stilbenoid, diarylheptanoid, and gingerol biosynthesis), carbohydrate metabolism (cutin, suberine, and wax biosynthesis; glycolysis/gluconeogenesis; and pyruvate metabolism), and lipid metabolism (fatty acid elongation, fatty acid biosynthesis, glycerophospholipid metabolism, and fatty acid degradation) (Table S8,S9).The expression of DEGs in Sub3 increased at 2900 m and decreased at 3300 m (Figure 8), and these DEGs were mainly enriched in signal transduction (such as plant hormone signal transduction), environmental adaptation (plant-pathogen interaction), amino acid metabolism (beta-alanine metabolism, tryptophan metabolism, and taurine and hypotaurine metabolism), and secondary metabolite metabolism (terpenoid backbone biosynthesis and carotenoid biosynthesis) (Table S8,S9).However, the DEGs in Sub4 were initially downregulated and then upregulated (Table S7), and they were mainly enriched in folding, sorting, and degradation (protein processing in the endoplasmic reticulum), secondary metabolite metabolism (stilbenoid,

Analysis of Gene Expression Pattern
To further explore the gene expression patterns in response to the increasing altitudes, we applied K-means clustering to the identified DEGs.In this analysis, 1466 DEGs with a known-function were grouped into six categories (referred to as Sub1, Sub2, Sub3, Sub4, Sub5, and Sub6), comprising 113, 616, 165, 220, 194, and 158 genes, respectively (Figure 8; Table S7).The DEGs in the Sub1 and Sub2 groups were characterized by constant upregulation with the increasing elevations (Table S7).The GO and KEGG analysis showed that these genes were mainly enriched in secondary metabolite metabolism (such as biosynthesis of secondary metabolites; flavonoid biosynthesis; and stilbenoid, diarylheptanoid, and gingerol biosynthesis), carbohydrate metabolism (cutin, suberine, and wax biosynthesis; glycolysis/gluconeogenesis; and pyruvate metabolism), and lipid metabolism (fatty acid elongation, fatty acid biosynthesis, glycerophospholipid metabolism, and fatty acid degradation) (Tables S8 and S9).The expression of DEGs in Sub3 increased at 2900 m and decreased at 3300 m (Figure 8), and these DEGs were mainly enriched in signal transduction (such as plant hormone signal transduction), environmental adaptation (plant-pathogen interaction), amino acid metabolism (beta-alanine metabolism, tryptophan metabolism, and taurine and hypotaurine metabolism), and secondary metabolite metabolism (terpenoid backbone biosynthesis and carotenoid biosynthesis) (Tables S8 and S9).However, the DEGs in Sub4 were initially downregulated and then upregulated (Table S7), and they were mainly enriched in folding, sorting, and degradation (protein processing in the endoplasmic reticulum), secondary metabolite metabolism (stilbenoid, diarylheptanoid, and gingerol biosynthesis), and energy metabolism (carbon fixation in photosynthetic organisms) (Tables S8 and S9).The last two groups (Sub5 and Sub6) were continuously downregulated with increasing elevations, and they were mainly enriched in secondary metabolite metabolism (biosynthesis of secondary metabolites; phenylpropanoid bio synthesis; flavonoid biosynthesis; and stilbenoid, diarylheptanoid, and gingerol biosynthesis), carbohydrate metabolism (carbon metabolism, glyoxylate and dicarboxylate metabolism, pentose phosphate pathway, and fructose and mannose metabolism), energy metabolism (photosynthetic antenna proteins, photosynthesis, carbon fixation in photosynthetic organisms, and nitrogen metabolism), and amino acid metabolism (biosynthesis of amino acids, and phenylalanine, tryptophan, glycine, serine, and threonine metabolism) (Tables S8 and S9).This clustering demonstrated that the genes involved in photosynthetic C fixation, carbohydrate metabolism, and N metabolism were significantly influenced by elevation.Hence, we further explored the DEGs involved in C and N metabolism.

DEGs Involved in Carbon Assimilation and Energy Metabolism
In this study, 25 of the 33 DEGs involved in C assimilation and metabolism were downregulated with the increasing elevations (Table 1; Table S10A).Among these 25 DEGs, almost all photosynthesis-related genes were inhibited by higher altitudes (3300 m), whereas CLH2, a gene that encodes chlorophyllase-2 (which degrades chlorophyll), significantly increased with the increasing elevations.Correspondingly, the expression of eight in the nine LHCBs or LHCA (chlorophyll a-b binding proteins) was downregulated, which decreased the ability of capturing and delivering excitation energy to photosystems.Meanwhile, 12 of the 13 photosystem genes (photosystem I complex/PSI and photosystem II complex/PSII), excluding the 4.0-5.7-foldincrease in photosystem I

DEGs Involved in Carbon Assimilation and Energy Metabolism
In this study, 25 of the 33 DEGs involved in C assimilation and metabolism were downregulated with the increasing elevations (Table 1; Table S10A).Among these 25 DEGs, almost all photosynthesis-related genes were inhibited by higher altitudes (3300 m), whereas CLH2, a gene that encodes chlorophyllase-2 (which degrades chlorophyll), significantly increased with the increasing elevations.Correspondingly, the expression of eight in the nine LHCBs or LHCA (chlorophyll a-b binding proteins) was downregulated, which decreased the ability of capturing and delivering excitation energy to photosystems.Meanwhile, 12 of the 13 photosystem genes (photosystem I complex/PSI and photosystem II complex/PSII), excluding the 4.0-5.7-foldincrease in photosystem I reaction center subunit N (PSAN) (MA_3482g0010), were negatively regulated, although ribulose bisphosphate carboxylase small chain (RBCS) (MA_476g0010), which is responsible for the primary event in carbon dioxide fixation, was upregulated.This indicated that photosynthetic capacity might be significantly inhibited at the altitude of 3300 m.Subsequently, central C metabolism was also diversely affected by the increasing elevations.Half of the genes were downregulated, including one F16P (MA_10432295g0010) encoding fructose-1,6-bisphosphatase, one KPRS1 (MA_10431983g0010) encoding ribose-phosphate pyrophosphokinase 1, and two G3PCs encoding glyceraldehyde-3-phosphate dehydrogenase, whereas others were upregulated, including one SDHA (MA_10427069g0010) encoding succinate dehydrogenase (ubiquinone) flavoprotein subunit, and two GADs (MA_161013g0010, MA_89828g0010) encoding glutamate decarboxylase.These results suggested that the expression of genes involved in photosynthetic C assimilation and energy metabolism was significantly altered by the increasing elevations.

DEGs Involved in Non-Struct(ural Carbon Metabolism
In order to examine the effects of elevations on NSC metabolism, three KEGG pathways (soluble sugars, fatty acid metabolism, and flavonoid biosynthesis) were examined (Table 2; Table S10B).With the increasing elevation, seven DEGs involved in soluble sugar metabolism were identified, including AGLU (α-glucosidase ), GOLS1 (galactinol synthase 1), three GAOAs (galactose oxidase), and two RFSs (probable galactinol-sucrose galactosyltransferase).AGLU, as a gene responsible for the breakdown of starch and disaccharides into glucose, was upregulated by more than 2.3-fold.Similarly, the expression of two RFSs, which are involved in the synthesis of raffinose, was also significantly increased at a certain elevation.However, GOLS1, involved in the biosynthesis of raffinose family oligosaccharides (RFOs) that function as osmoprotectants, was downregulated; three GAOAs were also downregulated.Interestingly, with the increasing elevation, we observed 60 DEGs involved in lipid metabolism, and 54 of them were upregulated, which might contribute to the survival in cold conditions at high altitudes.Moreover, flavonoids are important C-based secondary metabolites, which play a crucial role in the survival of trees under high altitude conditions.In this study, 29 genes involved in flavonoid biosynthetic pathway were identified as DEGs.Twelve of these genes were upregulated, including C4H (trans-cinnamate 4-monooxygenase), CAMT (caffeic acid 3-O-methyltransferase), CSE (caffeoylshikimate esterase), three DFRs (dihydroflavonol-4-reductase), F3 5 H (flavonoid 3 ,5 -hydroxylase), two F3 Hs (flavonoid 3 -monooxygenase), UFGT (anthocyanidin 3-O-glucosyltransferase) and BAN (anthocyanidin reductase), and they triggered the biosynthesis of anthocyanidin or condensed tannins, contributing to the resistance of trees to UV and other stresses.

DEGs Involved in Nitrogen Metabolism
With the increasing elevation, eight DEGs were found to be involved in N assimilation (Table 3; Table S10).Three NIAs (nitrate reductases), which are key genes involved in the first step of nitrate assimilation in plants, were downregulated.Meanwhile, one NIR (ferredoxin-nitrite reductase), which catalyzes the six-electron reduction of nitrite into ammonium, was downregulated.Moreover, two NRTs (nitrate transporters) were also downregulated, although one NRT and one AMT (ammonium transporter) were found to be upregulated.These results suggested that high altitude conditions blocked the uptake, assimilation, and transport of N.

Response of NSC and N to the Increasing Elevations
The elevational trend in plant tissue NSC (Figure 1) indicated that the tissue NSC concentrations of trees growing at the alpine treeline were lower than that of trees at lower elevations, in terms of the early growing season.This result suggested that the treeline trees may suffer from a 'C limitation' in the early growing season.In the evergreen trees, branch storage and remobilization were found to be crucial for new shoot growth in spring [34].We speculated that the concentrations of NSC in tissues were decreased in the early growing season, indicating that the mobile carbohydrates stored in tissues were utilized for shoot growth.On the other hand, plants grown under environmental stress conditions have lower tissue NSC concentrations than those grown under normal conditions.Many studies have shown that trees' hardness and ability to withstand cold are positively related to soluble carbohydrate concentrations in the tissues [10,35,36].Alpine trees grown under a cold environment consume large amounts of stored carbohydrates to maintain normal growth and resist to a low temperature stress.Our result is consistent with that work of Dang et al. [11], which showed that the NSC concentrations in tissues were highest in the late growing season and decreased from early to late growing season on both slopes, indicating that the mobile carbohydrates stored in tissues were utilized from bud break and shoot growth in the early growing season [11].
In this study, the increasing trends of total N in the needles and stems suggested that N availability did not constrain sink activity at higher altitudes, which is why the treeline trees had no sufficient nutrient supply.This result was in agreement with those of previous studies, which showed that the tissue N content in treeline trees is higher than that in lower elevation trees [10,20,[37][38][39].Reich et al. [40] also reported that plants at high elevations exhibit N-physiological adaptation strategies to compensate for lower efficiency of physiological processes in low temperature environments.Therefore, these studies indicated that N seemed not to be a limiting factor for tree growth and development around the alpine treelines.Our results also showed that N concentrations in the needles were higher than those in roots at all elevations, supporting the fact that N allocation to storage tissues and remobilization for reuse in trees are seasonally programmed, and that the annual growth of many tree species heavily relies on the internal cycling of N, by which N is stored over winter and subsequently remobilized in spring for new growth [41,42].However, whether N indeed limits the growth of treeline trees and determines the treeline formation needs to be further studied.

Changes of Genes Involved in C Assimilation and Metabolism in P. crassifolia with Increasing Elevations
In this study, many genes were found to be differentially expressed with the changes of elevation, including most of the downregulated genes involved in photosynthesis and carbohydrate metabolism.This result was consistent with the decrease in NSC content with the increasing elevation.All aspects of photosynthesis including the light reactions, the Calvin cycle, and photorespiration can be affected by low temperature stress [43,44], as well as the light-harvesting chlorophyll a/b-binding proteins, which are major constituents of the light harvesting complex in the thylakoid membranes [45].This study showed that the expression of six chlorophyll a-b binding protein genes were significantly downregulated with the increasing elevations, which helped to explain the decreased C fixation, also suggesting that cold soil and ultraviolet-B radiation (UV-B ) radiation leads to the downregulation of several transcripts encoding key photosynthetic proteins.
RFOs, particularly raffinose, play a role in the acquisition of cold tolerance and other stress responses in many plant species [46,47].α-Galactosidase (AGAL) is a key enzyme involved in the breakdown of RFOs [47], and plays a role in producing osmoprotectants that stabilize the membranes and protect cells against oxidative stress [48].Moreover, α-glucosidase (AGLU) digests carbohydrates, including starch and stable sugar [46].Thus, the increase in AGAL and AGLU expression observed in this study might help P. crassifolia to resist the cold stress at the higher elevations.However, the expression of GOLS1 (encoding galactinol synthase 1), which is the first committed enzyme in the RFO synthesis pathway [46], decreased with the increasing elevations.These results indicated that the RFO pathway might play a crucial role in the response of P. crassifolia to elevation changes.
Plant cuticle lipids form outer protective layers to resist environmental stresses, including water loss, UV radiation, and low temperature [49,50], and many of these related genes are linked by their biological function to a specific response to low temperature and strong UV [51,52].For example, 3-ketoacyl-CoA synthases (KCSs) belong to the group of "very long chain fatty acids", which are required for wax synthesis [53,54].Moreover, fatty acyl-coA reductases (FARs), which are reported to be involved in cuticular wax production, i.e., involved in wax synthesis [49], and long-chain acyl-CoA synthetases (LACSs), have important roles in lipid synthesis and storage, synthesis of cutin and wax, and stress responses [55].In this study, six KCSs, three FARs, and three LACSs were identified to be significantly induced with the increasing elevation, which might contribute to the survival of P. crassifoliain in the alpine environment.Furthermore, a significant enhancement in 17 of the 21 GDLs-encoding GDSL esterase/lipases, which play important roles in plant physiological processes, such as development, and stress responses [56]-would benefit the survival of P. crassifolia in the alpine environment.
Flavonoids are a large class of plant secondary metabolites that exert diverse biological and pharmacological effects, including the protection of plants against abiotic stresses [57].With respect to enzymatic coding genes involved in the flavonoid or anthocyanin biosynthetic pathway, in this study, 12 of the 29 DEGs were upregulated at different elevations.For example, F3Hs and F3 5 Hs, encoding flavanone 3-hydroxylases and flavonoid 3 ,5 -hydroxylase, respectively, are involved in an early step in flavonoid metabolism leading to the formation of dihydroflavonols from flavanones [57,58].In grape berries, UV-B increases the expression of these two genes, resulting in increased flavonol and anthocyanin content [59].Our study showed that half of these genes were upregulated, indicating their roles in the resistance to cold stress and UV stress at higher elevations.Meanwhile, the increase of elevations also significantly promoted the transcript of anthocyanin downstream pathway genes, including 3 DFRs (dihydroflavonol-4-reductase), 1 UFGT (anthocyanidin 3-O-glucosyltransferas), and 1 ANR (anthocyanidin reductase), and these would contribute to the resistance of trees to UV or other stresses [57].Therefore, the present study results indicated that flavonoids might play an important role in the survival of P. crassifoliain in the alpine environment.

Conclusions
This study demonstrated that the concentration of NSC significantly decreased in the current-year needles of P. crassifolia with increasing elevations, especially starch, indicating C limitation in trees around the treeline environment in the early growing season.This was further supported by the results of RNA-Seq, which showed almost all the photosynthesis-related genes were downregulated by the increasing elevations.Meanwhile, we found that large amounts of differentially expressed genes were involved not only in the central C metabolism, but also in the C-based secondary metabolism, including wax and flavonoid biosynthesis, contributing to the survival of P. crassifolia in the alpine environment.Thus, our study provided molecular evidence of C limitation hypothesis in treeline formation, as well as provided a better understanding of molecular mechanism of treeline tree survival in severely adverse conditions.

Figure 2 .
Figure 2. Changes of non-structural carbohydrates in the roots of Picea crassifolia with the increasing elevations (m).(A-C), coarse roots with 0.5-1.0cm in diameter; (D-F), fine roots with <0.5 cm in diameter.Different lowercase letters represented the significance levels (p< 0.05) among the different elevations.

Figure 1 .
Figure 1.Changes of non-structural carbohydrates of Picea crassifolia with the increasing elevations (m).(A-C) Current-year needles; (D-F) current-year branches.Different lowercase letters represented the significance levels (p < 0.05) among the different elevations.NSC: non-structural carbohydrate.

Figure 2 .
Figure 2. Changes of non-structural carbohydrates in the roots of Picea crassifolia with the increasing elevations (m).(A-C), coarse roots with 0.5-1.0cm in diameter; (D-F), fine roots with <0.5 cm in diameter.Different lowercase letters represented the significance levels (p< 0.05) among the different elevations.

Figure 2 .
Figure 2. Changes of non-structural carbohydrates in the roots of Picea crassifolia with the increasing elevations (m).(A-C), coarse roots with 0.5-1.0cm in diameter; (D-F), fine roots with <0.5 cm in diameter.Different lowercase letters represented the significance levels (p< 0.05) among the different elevations.

Figure 3 .
Figure 3. Changes of total N in the different tissues of Picea crassifolia with the increasing elevations (m).(A) Current-year needles, (B) current-year branches, (C) coarse roots with 0.5-1.0cm in diameter, (D) fine roots with <0.5 cm in diameter.Different lowercase letters represented the significance levels (p < 0.05) among the different elevations.

Figure 3 .
Figure 3. Changes of total N in the different tissues of Picea crassifolia with the increasing elevations (m).(A) Current-year needles, (B) current-year branches, (C) coarse roots with 0.5-1.0cm in diameter, (D) fine roots with <0.5 cm in diameter.Different lowercase letters represented the significance levels (p < 0.05) among the different elevations.

Figure 7 .
Figure 7. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in Picea crassifolia with elevations.

Figure 8 .
Figure 8. Clustering analysis by K-means clustering of differentially expressed genes (DEGs) in Picea crassifolia with elevations.

Figure 8 .
Figure 8. Clustering analysis by K-means clustering of differentially expressed genes (DEGs) in Picea crassifolia with elevations.

Table 1 .
DEGs involved in carbon assimilation and carbon metabolism.

Table 2 .
DEGs involved in non-structural carbon metabolism.

Table 3 .
DEGs involved in nitrogen metabolism.