De novo sequencing of the transcriptome of Fargesia macclureana (Poaceae) reveals regulators of the floral transition and ecological adaptations to high altitude

Background Fargesia macclureana (Poaceae) is a woody bamboo species found on the Qinghai–Tibet Plateau (QTP) approximately 2,000 ~ 3800 m above sea level. It rarely blossoms in the QTP, but it flowered 20 days after growing in our lab, which is in a low-altitude area outside the QTP. To date, little is known regarding the molecular mechanism of bamboo flowering, and no studies of flowering have been conducted on wild bamboo plants growing in extreme environments. Here, we report the first de novo transcriptome sequence for F. macclureana to investigated the putative mechanisms underlying the flowering time control used by F. macclureana to adapt to its environment. Results Illumina deep sequencing of the F. macclureana transcriptome generated 140.94 Gb of data, assembled into 99,056 unigenes. A comprehensive analysis of the broadly, specifically and differentially expressed unigenes (BEUs, SEUs, and DEUs) and a weighted gene co-expression network analysis (WGCNA) revealed that changes in expressions of unigenes related to the circadian cycle may account for the differences in the floral transition of F. macclureana after being transplanted from the QTP to a laboratory outside. In addition, there were differences in active carbohydrate metabolism and signal transduction between the flowering and non-flowering plants. Moreover, we detected the expression of unigenes related to DNA repair and plant-pathogen interactions, which may be of adaptive importance. Finally, we detected 9,296 simple sequence repeats (SSRs) that may be useful for further molecular marker-assisted breeding. Conclusions F. macclureana may have evolved specific reproductive strategies for flowering-related pathways in response to photoperiodic cues to ensure long vegetation growing period. Our findings will provide new insights to future investigations into the mechanisms of flowering time control and adaptive evolution in plants growing at high altitudes. F-leaves plant

In addition, there were differences in active carbohydrate metabolism and signal transduction between the flowering and non-flowering plants. Moreover, we detected the expression of unigenes related to DNA repair and plant-pathogen interactions, which may be of adaptive importance. Finally, we detected 9,296 simple sequence repeats (SSRs) that may be useful for further molecular marker-assisted breeding. Conclusions F. macclureana may have evolved specific reproductive strategies for flowering-related pathways in response to photoperiodic cues to ensure long vegetation growing period. Our findings will provide new insights to future investigations into the mechanisms of flowering time control and adaptive evolution in plants growing at high altitudes. 3

Background
The flowering time is of crucial importance to ensure the reproductive success of flowering plants. Previous results have indicated that the floral transition is orchestrated by several parallel and interactive genetic pathways that are regulated by a variety of environmental and endogenous signals [1]. Many key genes and regulatory networks have been identified in herbaceous annual plants such as Arabidopsis [2,3], rice [4], gourds [5], potato [6] and sorghum [7]. However, much less is known about such regulation in perennial plants.
Despite the increasing attention on perennial dicotyledonous woody plants such as poplar [8,9], eucalyptus [10] and citrus [11] species, to date, the molecular mechanism underlying floral regulation in monocotyledonous woody plants remains elusive.
Furthermore, previous studies investigated flowering mainly by artificially altering the external signals (e.g. photoperiod and light intensity) and did not assess the impact of the original environment on the adaptive evolution of species-specific reproductive strategies.
Bamboo plants are an important group in the Bambusoideae subfamily of the monocotyledonous Poaceae. They exhibit a wide degree of variation in the timing (1-120 years) and nature (sporadic vs. gregarious) of flowering among species [12]. Sporadic flowering involves flowering in only a few isolated clumps, which set little or no seed and usually remain alive afterward [13]. In contrast, gregarious flowering involves all individuals of a species regardless of age and/or location within and among the populations at the same time, which is usually followed by death and seed setting [14].
And the simultaneous death of many individuals triggers serious ecological consequences, including changes in the population dynamics of neighboring plants, differences in soil properties, various effects on endangered animals that depend on bamboo [15], and the knock-on effects on human economies in many parts of the world [16]. Therefore, dissecting the regulators that control the unique life history of bamboo may be of use for 4 plant ecology and human society. However, to date, little is known regarding the molecular mechanisms of bamboo flowering, in part because of the sporadic occurrence of these flowering episodes and the long intervals between events.
Many genes have been identified as regulators of reproductive development in different bamboo species, including the MADS-box transcription factors [17][18][19], CONSTANS ( CO) [20] and FLOWERING LOCUS T (FT) [21], among others. In addition, studies of sequenced transcriptomes have identified microRNAs related to floral development [22][23][24]. However, samples collected in these analyses were limited to mature spikelets or to different spikelets at different development stages. Thus, it is likely that dynamic changes in genes occurring at different development stages may be missing. In addition, the specific response of particular tissues to internal and external cues and how plants integrate these signals to regulate different phases of reproductive development (including the floral transition, florigen transport, and floral organ specification) has not yet been elucidated in bamboo. Furthermore, no studies of flowering have been conducted on wild bamboo plants growing in extreme environments.
Here, we took advantage of an unexpected flowering event in highland arrow bamboo, Fargesia macclureana [25], and performed the first de novo transcriptome analysis. This transcriptome includes data from six different tissues collected at different development stages, including inflorescences in the initial and peak flower stage (I-and P-spikelets), branchlets, and leaves from both flowering and non-flowering bamboo plants (F/NFbranchlets and F/NF -leaves). F. macclureana is a woody bamboo species found in areas 2,000 ~ 3,800 m above sea level on the Qinghai-Tibet Plateau (QTP) (Fig. 1), which is the highest and largest plateau in the world. The growth environment of the QTP is characterized by low temperature and low oxygen availability, reduced pathogen incidence, and intense radiation [26]. F. macclureana rarely blossoms in the QTP, but it 5 flowered 20 days after growing in our lab, which is in a low-altitude area outside the QTP.
Our goal is to use the transcriptomic data to gain a deeper understanding of the mechanisms underlying the control of flowering time and the adaptation of F. macclureana to the complex extreme conditions of the QTP. On one hand, we expect to detect regulatory hubs involved in the flowering mechanisms. On the other hand, we aim to discover signs of the adaptive evolutionary changes in F. macclureana in response to the harsh environmental conditions in the QTP, which may, in turn, provide a broader insight into the adaptive mechanisms for plants that grow at high altitudes.

De novo transcriptome assembly yielded 99,056 unigenes
Illumina deep sequencing of the F. macclureana transcriptome generated 140.94 Gb of data, including 471,537,304 clean reads in 18 unique samples (Additional file 1: Table S1).
The average Q20 (sequencing error rate less than 1%) and Q30 (sequencing error rate less than 0.1%) percentages were 100.00% and 89.95% respectively. The GC content of all samples ranged from 53.78% to 55.86%, with an average of 54.81%. Sample data were assembled into 289,122 transcript scaffolds, with an N50 and average length of 1,765 bp and 1,183 bp, respectively. The final de novo assembly included 99,056 unigenes, with an N50 and average length of 1,587 bp and 926 bp, respectively. Among these unigenes, 71.02% (70,354) were shorter than 1,000 bp and 12.06% (11,950) were longer than 2,000 bp (Table 1).

Most unigenes were functionally annotated and classified
A total of 47,306 unigenes were annotated (Additional file 2: Table S2). Of these, 45,516 (96.22%) unigenes were found to encode products that showed significant similarity to characterized proteins in the non-redundant protein sequence database (Nr) at an E-value 6 threshold of 10 -5 (Table 2). We also found that 7,027 (15.45 %) unigenes showed similarity to genes found in rice, 11.33% were similar to those found in Brachypodium distachyon, and we also found a significant proportion of the unigenes that were similar to those found in Setaria italica, Oryza brachyantha, and Zea mays (Fig. 2a). We identified 24,847 (52.52%), 28,317 (59.86%) and 43,909 (92.82%) unigenes that showed significant matches to entries in the Swiss-Prot, Pfam, and eggnog databases, respectively ( Table 2).
Many unigenes expressed in the F. macclureana transcriptome were functionally annotated as regulators of plant responses to evolutionarily important phenotypes, including membrane stabilization, heat stress response and pathogen defense (Additional file 2: Table S2).

Functional annotation indicated that many unigenes were involved in metabolism and genetic information processing
We were able to annotate 13,128 unigenes (27.75% of the total) in 25 different categories of the COG (clusters of orthologous groups) classification database (Fig. 2b). Of these, the cluster for "General function prediction only" (3,277, representing 24.96% of the 13,128 unigenes annotated by this database) was the largest group, followed by "Replication, recombination and repair" (2,202, 16.77%), "Transcription" (1,571, 11.97%), and "Translation, ribosomal structure and biogenesis" (1,429, 10.88%). The "Signal transduction mechanisms", "posttranslational modification, protein turnover, chaperones", "carbohydrate and amino acid transport and metabolism" and "transport and metabolism" categories also contained a significant proportion of the annotated unigenes.
GO enrichment analysis indicated that these predicted unigenes were categorized into three main categories-i.e. biological process (BP), cellular component (CC), and molecular function (MF). As shown in Fig. 2c, for unigenes that were enriched in the BP category, they were mainly involved in biological processes related to reproduction, 7 posttranslational modification and signal transduction; as for those in the CC category, they were mainly involved in cellular components related to membrane, ubiquitin ligase complex, mitochondrion, chloroplast and etc.; while for those in the MF category, they were mainly involved in molecular functions related to signaling transduction (e.g. "ATP binding", "zinc ion binding", "protein kinase activity", and etc.) (Additional file 3: Table   S3).
We also mapped 14,307 unigenes (representing 30.24% of the total) to six different KEGG subsystems, including metabolism, genetic information processing, environmental information processing, cellular processes, and organismal systems. As shown in Fig. 3, the majority of these unigenes (7,922, representing 66.17% of the 14,307 unigenes classified using KEGG annotations) were assigned to metabolic pathways, including carbohydrate metabolism, energy metabolism, and others. In addition, 4,024 unigenes (28.13%) were assigned to genetic information processing, including transcription, translation, and folding, and 474 unigenes (3.31%) were found to be related to membrane transport and signal transduction. We also found 707 genes (4.94%) that were related to transport and catabolism and 377 genes (2.64%) related to environmental adaptation.

Most BEUs were involved in genetic information processing, environmental adaptation and signal transduction
As shown in the Venn diagram (Fig. 4a), we found nearly equal numbers of unigenes that were broadly and specifically expressed in I-spikelets, P-spikelets, F-branchlets, and Fleaves. COG analysis indicated that most BEUs were clustered in signal transduction mechanisms (T), replication, recombination and repair (L), and transcription (K), besides general function prediction only (R). GO enrichment analysis for these BEUs indicated that they were also mainly involved in reproduction, environmental adaptation and signal transduction, which was largely similar with that for all predicted unigenes (Additional file 8 4: Table S4-a).
KEGG enrichment analysis also indicated that these BEUs were mainly enriched in pathways related to environmental adaptation (including circadian rhythm, endocytosis, and plant-pathogen interactions), signal transduction (including plant hormone signal transduction, phosphatidylinositol signaling system, and inositol phosphate metabolism) and genetic information processing (including spliceosome, mRNA surveillance, and RNA transport and degradation; Additional file 4: Table S4-b).
The SEUs were mostly involved in carbohydrate metabolism, energy metabolism,

and environmental adaptation
As shown in Fig. 4a, we identified 10,653 unigenes that were specifically expressed in spikelets, including 5,528 and 5,025 unigenes in I-and P-spikelets, respectively. We also found 9,067 and 7,437 unigenes that were specifically expressed in F-branchlets and Fleaves, respectively. COG annotation indicated that the distribution patterns of SEUs among the 26 terms were similar, with the number of SEUs within each term varying among the three tissues (Fig. 4b).
The GO enrichment analysis indicated that these SEUs not only shared some common GO terms, but also had some particular ones. As shown in Fig. 4c and Additional file 4: Table   S4-c, for those SEUs that were enriched in the BP category, they were broadly involved in several important biological processes, including "protein phosphorylation", "regulation of flower development", "protein ubiquitination", "regulation of transcription, DNAtemplated", "reciprocal meiotic recombination" and "meiotic chromosome segregation". In addition, SEUs in I-and P-spikelets were also involved in some processes related to reproduction; and those in F-branchlets were mainly involved in processes related to posttranslational modification; while those in F-leaves were mainly involved in processes related to plant-pathogen interaction. As for those in the CC category, they were broadly 9 involved in several important cellular components, including "mitochondrion", "plasma membrane" and "plastid". In addition, SEUs in I-and P-spikelets were also involved in ribosome and mitochondria; and those in F-branchlets were mainly involved in endoplasmic reticulum and proteasome; and those in F-leaves were mostly involved in chloroplast. As for those in the MF category, they were broadly involved in several molecular functions, including "ATP binding", "ubiquitin-protein transferase activity" and "protein tyrosine kinase activity". In addition, SEUs in I-and P-spikelets were also involved in DNA and microtubule binding; those in F-branchlets were also enriched in oxidoreductases activities; and those in F-leaves were also enriched in enzymes involved in carbohydrate metabolism.
As shown in Additional file 5: Fig. S1, KEGG pathway analysis indicated that SEUs in I-and F-spikelets mainly mapped to the ribosome pathway, with those in F-branchlets mainly mapped to the ribosome, amino acid biosynthesis, and carbon metabolism pathways, and those in F-leaves mainly mapped to KEGG pathways related to energy metabolism (including oxidative phosphorylation, fatty acid metabolism, and photosynthesis), environmental adaptation (e.g. proteasomes), genetic information processing, and various unrelated metabolic pathways (e.g. tryptophan metabolism, beta-alanine metabolism, and N-glycan biosynthesis).

DEUs were mostly involved in carbohydrate and energy metabolism, signal transition and environmental adaptation
As shown in Table 3, many unigenes showed differential expressions across all 15 groups sampled. The number of DEUs in each sample pair ranged from 970 between I-vs Pspikelets to 13,577 in NF-leaves vs I-spikelets. For most pairwise comparisons, the number of up-and down-regulated DEUs was approximately the same, except for four groups, including I-vs P-spikelets, F-branchlets vs both I-and P-spikelets, and F-leaves vs P-spikelets.
The Venn diagram of DEU sets shows that 5,494 unigenes were differentially expressed in F-branchlets/F-leaves vs I-and P-spikelets. For those DEUs that were up-regulated in spikelets, they are mainly mapped to KEGG pathways related to carbohydrate metabolism, plant-pathogen interactions and DNA repair (Fig. 5a). Notably, among the 970 DEUs identified between I-and P-spikelets, 916 up-regulated DEUs were mapped to KEGG pathways related to metabolic activity (Additional file 6: Table S5).
A total of 5,494 unigenes were differentially expressed in the DEU sets of spikelets/Fleaves vs F-branchlets. Upregulated DEUs in F-branchlets were mapped to KEGG pathways including phenylalanine metabolism, phenylpropanoid biosynthesis, ABC transporters, and flavone and flavonol biosynthesis (Fig. 5b). Those that were upregulated in F-and NFleaves vs F-branchlets were mainly mapped to plant hormone signal transduction, homologous recombination, base excision repair, and mismatch repair (Additional file 6: Table S5). Notably, 3,275 (50.20% of the total) DEUs found between NF-and F-branchlets were upregulated; these were mainly mapped to KEGG pathways related to replication and recombination (Additional file 6: Table S5). Those that were downregulated were mainly mapped to carbon fixation and photosynthesis (Additional file 6: Table S5).
We also found that 6,966 (43.69% of the total) DEUs found in spikelets/F-branchlets vs Fleaves were up-regulated, and were mainly mapped to KEGG pathways related to carbohydrate metabolism (Fig. 5c). 2,492 (49.52%) DEUs in NF-vs F-leaves were upregulated, and these were mainly mapped to starch and sucrose metabolism (Additional file 6: Table S5). In contrast, downregulated DEUs were mainly mapped to KEGG pathways related to photosynthesis (Additional file 6: Table S5).

WGCNA results identified gene modules related to specific tissues
WGCNA results showed that unigenes expressed in the six different tissues of flowering and nonflowering plants tested here clustered into 18 branches representing 18 different genetic modules (Additional file 11: Fig. S5a). Unigenes within each module were highly co-expressed, while those in different modules were co-expressed to a lower degree (Additional file 11: Fig. S5b). In six of the samples collected, we identified nine significant gene modules including 1,344 unigenes. Here, correlation coefficient of a module with a related trait > 0.7 was used as a threshold of significance (Additional file 11: Fig. S5c).
Notably, these six tissues were more strongly divided into clades according to whether they were flowered or not rather than by the differences among tissues (Additional file 11: S5d).
In addition, the unigenes in gene modules relating to I -and P-spikelets were most strongly enriched in KEGG pathways related to carbohydrate metabolism, genetic information processing, and environmental information processing. In contrast, those related to F-and NF-branchlets were mostly enriched in KEGG pathways related to metabolism, plant hormone signal transduction, and genetic information processing. The gene modules related to F-leaves were enriched in pathways related to plant hormone signal transduction and protein processing, while the gene modules related to NF-leaves were enriched in KEGG pathways related to oxidative phosphorylation (Additional file 12: Table S7).

Activated Hd3a expression probably accelerates flowering in F. macclureana
FT is a key floral regulator that controls the timing of flowering and seasonal growth cessation in response to light and the circadian clock in many plant species [8,10,27]. In this study, FmHd3a, a bamboo ortholog of rice FT, was significantly expressed only in tissue samples collected from flowering plants. In rice, Hd3a functions as a major photoperiodic flowering regulator and participates in the OsGI-Hd1-Hd3a module, which is similar to the GI-CO-FT module in Arabidopsis [28]. Hd1 activates and suppresses Hd3a expression by promoting heading under the short day (SD) and long day (LD) conditions, respectively [29][30]. As F. macclureana rarely blossoms on the QTP, which experiences a long photoperiod with a low ratio of red to far-red light, it may have evolved specific reproductive strategies involving flowering-related pathways in response to photoperiodic cues to ensure long vegetation growing period. It is probably that the weak light intensity with a low proportion of blue light might activate Hd3a expression even in the LD conditions, thereby accelerating flowering. Notably, reproduction pathways play an important role in the mechanisms of plant adaptation to extreme environments. Previous 13 studies showed that the phytochrome and flowering time regulatory protein 1 (PFT1) from Crucihimalaya himalaica, a close relative of Arabidopsis and Capsella, grows on the QTP, showed signs of positive selection for adaptive divergence [31]. We presume that the function of Hd3a in promoting flowering is likely to be conserved between bamboo and rice, because both of them belong to the Poaceae.
In F-leaves, the expression of FmFT and the photoreceptor gene FmCRY (c105898.graph_c2) were both significantly downregulated, while the expression of another FT ortholog, FmHd3a, was significantly upregulated. Photoreceptors mediate light input pathways to synchronize the circadian clock [1,5]. In A. thaliana, CRY activates FT transcription in response to blue light [3,32]. UV-B radiation causes a multitude of lowand high-fluence responses similar to the phytochrome responses [33][34]. Thus, downregulated FmFT expression is likely due to down-regulated FmCRY, which is, in turn, a response to a lower ratio of blue light or reduced light intensity (both in the laboratory or in the QTP). The upregulated expression may indicate that FmHd3a can function in a CRYindependent manner or be negatively regulated by FmFT. We suspect that the floral transition of F. macclureana is regulated by a complex regulatory network in which at least two unique FT orthologs interact with the circadian clock pathways. However, how these circadian clock pathways mediate the activation of FmHd3a and FmFT in response to light signalling remains to be elucidated by future research.
Notably, we detected FmFT expressions in all four tissues collected from the flowering plants, but not in NF-leaves or NF-branchlets. Given that all plants were grown in the same conditions, we suspect the physiological states of the plant itself may be the possible cause. Additionally, the DEUs upregulated in F-leaves relative to NF-leaves were mainly enriched in KEGG pathways related to starch, sucrose, and galactose metabolism, and corresponding down-regulated DEUs were mapped to the light and carbon fixation, plant circadian rhythm, and photosynthesis pathways. Therefore, we speculate that bamboo FT orthologs might be regulated by regulators involved in other pathways. However, the nature of the mechanism responsible for this cross-regulation needs further experimental verification.

Carbohydrate metabolism and signal transport may be major factors in floral transition, organogenesis, and death after flowering
Bamboo exhibits excellent flexibility and fracture toughness, and so far, the presence of fibers within the bamboo culm was thought responsible for these remarkable mechanical properties [35]. And the development of plant fibers is accompanied by the of carbohydrate metabolism [36]. Perhaps this is the reason why many DEUs were involved in carbohydrate metabolism pathway. Interestingly, our results indicated that starch and sucrose metabolism was a major enriched KEGG pathway for the DEUs from several combination pairs, including NF-vs F-leaves, branchlets/leaves vs I-& P-spikelets.
Transcripts and metabolic signatures of maize leaves have shown that the balance between transitory starch and sucrose is associated with the autonomous floral transition [37]. And in Lilium, carbohydrates have been found to be transported from the vascular bundles to floral organs during reproduction [38]. Yang et al., (2017) also reported that the deficiency in the resources in male flowers reduced pollen viability in Tapiscia sinensis due to biased carbohydrate transport toward the female flowers [39]. Therefore, we suspect there may be a correlation between DEUs related to starch and sucrose metabolism and arrow bamboo floral organ development.
In rice, excessive uridine 5´-diphosphoglucose-glucose (UDPG) can result in programmed cell death, accumulation of reactive oxygen species and an increase in the caspase-like activity [40] and inactivate starch synthase disrupted normal male reproduction by delaying programmed cell death in cotton [41], suggestive of a correlation between starch and sucrose metabolism and death. Bamboo flowering, especially in masting species, often causes plants to wilt and die after setting seed. It is possible that increased starch and sucrose metabolism might trigger the excessive accumulation of reactive oxygen species and result in the altered activity of key enzymes in important biological pathways.
In the present study, unigenes related to the signal transduction pathways were significantly upregulated in the tissues of flowering arrow bamboo plants. The transcriptomic profiles of Posidonia oceanica also showed a strong metabolic activation of hormones in the heat stress-induced flowering plants [42]. Signal transduction-related genes were also found to have undergone both significant positive selection and expansion events on the adaptive evolution of Crucihimalaya himalaica [31] and cyanobacterium Trichormus sp. NMC-1 [43] on the QTP. In the present study, unigenes related to the signal transduction pathways were significantly upregulated in the tissues of flowering arrow bamboo plants. We suspect that this may be due to the long distance transport of the FT proteins, which ensures floral promotion at the shoot apex [44], or the phytohormone signaling and calcium signaling, which play diverse roles in the specification of flower organs during arrow bamboo reproductive development [45][46].

F. macclureana has presumably evolved an integrated mechanism to adapt to the harsh environment of the QTP
We detected the broad expressions of unigenes encoding putative Hsp proteins, such as heat shock protein 70 ( Hsp70) and heat shock protein 90 ( Hsp90). Both Hsp70 and Hsp90 are important for maintaining cellular protein homeostasis under stress conditions and they function by activating other targets [47][48][49]. It is likely that the sudden exposure to the higher temperature of the lab (i.e. outside the QTP) triggered the expressions of these unigenes. Warmer temperatures can greatly reduce flowering synchrony among individuals from 72 woody and herbaceous species [50]. It was also reported that P. oceanica, a highly clonal and long-lived species, massively bloomed after a simulated heatwave [42]. Given the cold temperatures present at the high-altitude regions of the QTP [26], it is reasonable to presume that F. macclureana has developed into a heatsensitive but not heat-tolerant bamboo species and flowering is probably a stress-induced response to the higher temperature in the lab.
DNA repair and disease-resistance pathways have been found to play a crucial role in the highland adaptation of Tibetan highland plants [43,51]. In the present study, we detected significant expressions of many unigenes related to pathogen response that contained either a nucleotide binding (NB)-ARC domain or a leucine-rich repeat (LRR) domain, which were present in most resistance (R) proteins [52][53][54]. We also identified many differentially expressed unigenes that were significantly enriched in the DNA repair pathways. Since relatively fewer species of pathogenic microorganisms and intense UV radiation exist on the QTP [55], it is reasonable to presume that F. macclureana has evolved a relatively narrow range of pathogen specificity and specific DNA-repair mechanisms. Sudden exposure to the lab environment, which contains a heavier load of pathogens and weak light intensity, may have induced an innate defensive response of F. macclureana, and those that were enriched in the plant-pathogen interaction and DNA repair response pathways may be important for F. macclureana to cope with the new environment present in the lab. Although further studies are needed to investigate the molecular mechanisms responsible for the putative adaptive evolutionary changes, this study provides insights into how plants adapt to harsh and extreme environments.

Conclusions
In the present study, we constructed a novel de novo transcriptome analysis for F.

Tissues collection
The studied plant species is highland arrow bamboo (Fargesia macclureana), and it grows mainly as a underbrush of coniferous forest or coniferous and broad-leaved mixed forest, and sometimes forms a pure population in the QTP at an altitude of approximately 2,000 3 800 m above sea level (Fig. 9).

RNA extraction, quantification, and qualification
Total RNA was extracted from each of the six unique tissues mentioned above using a RNeasy plant RNA extraction kit (Qiagen, Dusseldorf, Germany), and the extraction procedure was performed according to the manufacturer's instructions. RNA degradation and contamination were monitored using 1% agarose gels. RNA purity was checked using a NanoPhotometer ® spectrophotometer (Implen GmbH, Munich, Germany). RNA concentration was measured using a Qubit ® RNA Assay Kit and a Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA). RNA integrity was assessed using an RNA Nano 6000 Assay Kit run on an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library preparation for transcriptome sequencing
Library construction and RNA-Seq were performed by the Biomarker Biotechnology Corporation (Beijing, China). A total of 3 μg RNA per sample was used for RNA preparation.
Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads, followed by fragmentation carried out using divalent cations at elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase (RNase H-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H.
The remaining overhangs were converted into blunt ends via the exonuclease and polymerase activities. Next, the 3´ ends of the DNA fragments were adenylated and ligated to the NEBNext adaptors with hairpin loop structures to prepare samples for hybridization, this was to select cDNA fragments that are 150-200 bp in length. Library fragments were then purified using an Agencourt AMPure XP system (Beckman Coulter, Brea, CA, USA), and 3 μl USER enzyme (New England Biolabs, Ipswich, MA, USA) was added to the size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by five min at 95°C before PCR. PCR was performed using Phusion High-Fidelity DNA polymerase (Thermo Fisher, Waltham, MA, USA), universal PCR primers, and the Index (X) Primer. Finally, the PCR products were purified using the AMPure XP system and library quality was assessed on the Agilent Bioanalyzer 2100.

Clustering and sequencing
The clustering of index-coded samples was performed using a cBot Cluster Generation System and a TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA), and all experimental procedures were performed according to the manufacturer's instructions.
After that, library preparations were sequenced on an Illumina Hiseq 2000 platform and paired-end reads were generated.

De novo transcriptome assembly
Raw data in fastq format were first processed using in-house perl scripts. Clean data were obtained by removing low-quality reads and reads that contain the adapters or poly-N sequences. Meanwhile, we checked the quality of our unassembled read dataset by 20 examining various measures including Q20, Q30, GC-content, and sequence duplication.
All the downstream analyses were performed using high-quality clean data.
The transcriptome was assembled using clean reads from all libraries and samples. The assembly was produced using Trinity [57] with min_kmer_cov set to 2 and all other parameters set to their respective default values.

Quantification of gene expression levels
Gene expression levels were estimated using RSEM [58] for each sample: clean data were mapped back onto the assembled transcriptome, followed by a read count for each gene.
The expression levels of unigenes were calculated and normalized using FPKM (fragments per kb per million fragments) [59].

Expression analysis of broadly and specifically expressed unigenes
For all unigenes, those that were expressed in all six tissues were defined as broadly expressed unigenes (BEUs). Similarly, unigenes that were specifically expressed in only one tissue were defined as specifically expressed unigenes (SEUs). The identification of BEUs and SEUs was conducted by using tools on the BMKCloud platform (http://www.biocloud.net).

Weighted gene co-expression network analysis (WGCNA)
WGCNA was performed on all unigenes identified using the WCGNA R package. We calculated the adjacency matrices and performed the topological overlap measures 21 (TOMs), which show the degree of overlap in shared neighbors between pairs of genes in the network to define gene clusters in our transcriptome dataset. 1 − TOM was used as a dissimilarity measure for hierarchical clustering and module detection. Modules of the clustered genes were then selected using the Dynamic Tree Cut algorithm as implemented by WGCNA. To identify modules that are significantly related to particular tissues, expression profiles of each module were summarized by a module eigengene (ME) and the correlations between the modules and corresponding tissues were calculated.

Expression analysis of differently expressed unigenes (DEUs)
Before analysis, we conducted a principal component analysis (PCA) and removed one replicate that showed an inconsistent expression pattern in the NF-branchlets and NFleaves to ensure consistency in the expression patterns of unigenes between replicates (Additional file 14: Figure S6).
Expression analysis of the DEUs between pairs of tissues/groups was performed using the DESeq package (1.10.1) in R. DESeq provides statistical routines for identifying differential expression in the digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini-Hochberg method for controlling the false discovery rate [60]. Here, uni-transcripts with an absolute value of log 2 ratio ≥ 2, an FDR significance score < 0.01, and an adjusted P-value < 0.05 were deemed to be differentially expressed.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis
To understand the higher-level functions of the observed unigenes, we performed GO term annotation and KEGG pathway enrichment analysis using BMKCloud (http://www.biocloud.net/; [61]). We used KOBAS 2.0 [62] to test the statistical enrichment of differentially expressed genes in KEGG pathways. Pathways with P values < 0.05 were 22 considered significantly enriched.

Protein-protein interactions (PPIs)
The DEU and SEU sequences were queried using BLASTX against the related species to predict PPIs that the DEUs and SEUs may be involved in. This search procedure was capable of identifying PPIs that may be similar to any others found in the STRING database (http://string-db.org/). These PPIs were then visualized using Cytoscape [63].

Detection of SSRs
Picard-tools version 1.41 and samtools version 0.1.18 were used to sort data, remove duplicated reads, and merge the bam alignment results of each sample. SSRs were identified using MISA (https://webblast.ipk-gatersleben.de/misa/).

Validation of FmHd3a transcript levels by qRT-PCR
To verify the expression pattern of the FmHd3a, we used RT-qPCR to assess the expressions of FmHd3a in six distinct tissues. First-strand cDNA was synthesized from total RNA extracted by using a reverse transcription system (Promega, Madison, WI, USA) following the manufacturer's instructions. Each RT-qPCR amplification was performed at least three times, and NTB and TIP41 were used as internal controls [64]. Primers for these genes are listed in Additional file 15: Table S9. The relative expression levels of FmHd3a in different tissues were calculated using the 2 -ΔΔCT method [65]. The statistical significance of differences in the mean levels of expression was tested using a one-way ANOVA. Significant differences in transcript abundance between different tissues were then compared using Duncan's multiple range tests as implemented by SPSS version 17.0 (IBM SPSS, Chicago, USA). We considered mean differences at P < 0.05 and P < 0.01 to be statistically significant and highly statistically significant, respectively.

24
The RNA sequencing dataset generated during the current study have been submitted to NCBI Sequence Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra) with the accession number PRJNA544133.

Ethics approval and consent to participate
Not applicable.