Funct Integr Genomics DOI 10.1007/s10142-014-0419-7
ORIGINAL PAPER
Transcriptome analysis reveals diversified adaptation of Stipa purpurea along a drought gradient on the Tibetan Plateau Yunqiang Yang & Xiong Li & Xiangxiang Kong & Lan Ma & Xiangyang Hu & Yongping Yang
Received: 22 April 2014 / Revised: 7 November 2014 / Accepted: 17 November 2014 # Springer-Verlag Berlin Heidelberg 2014
Abstract Natural selection drives species adaptations to biotic and abiotic stresses. Species distributed along a moisture gradient, such as Stipa purpurea, a dominant grass in alpine arid and semi-arid meadows on the Tibetan Plateau, provide an opportunity to evaluate the effects of long-term adaptation to differing degrees of drought stress on gene expression. However, the genetic basis of this divergence remains largely unknown. Next-generation sequencing technologies have provided important genome-wide insights on the evolution of organisms for which genomic information is lacking. To understand how S. purpurea responds to drought stress, we selected five populations distributed along the degressive rainfall line on the northwestern Tibetan Plateau that currently present evolutionary acclimation to localized drought pressure at the physiological and biochemical levels and compared their transcriptome responses. In addition, we performed de novo assembly of the S. purpurea transcriptome using short read sequencing technology and successfully assembled 84,298 unigenes from approximately 51 million sequencing reads. We quantified gene expression level to compare their transcriptome responses using mRNA-Seq and identified differentially expressed transcripts that are involved in primary Electronic supplementary material The online version of this article (doi:10.1007/s10142-014-0419-7) contains supplementary material, which is available to authorized users. Y. Yang : X. Li : X. Kong : L. Ma : X. Hu : Y. Yang (*) Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Science, Kunming 650204, China e-mail: [email protected] Y. Yang : X. Li : X. Kong : L. Ma : X. Hu : Y. Yang The Germplasm Bank of Wild Species, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming, Yunnan 650201, China X. Li : L. Ma University of Chinese Academy of Sciences, Beijing 100049, China
and secondary plant metabolism, plant hormone synthesis, defense responses, and cell wall synthesis. Furthermore, physiological and biochemical evidence supports that abscisic acid (ABA) accumulation and cell wall strengthening derived from the differential transcripts contribute to the tolerance of S. purpurea to drought stress. The mechanisms by which S. purpurea adapts to drought stress provide new insight into how plants ecologically adapt and evolve. Keywords Transcriptome . Stipa purpurea . Gradient drought . Tibetan Plateau
Introduction With a mean altitude of more than 4000 m, the Tibetan Plateau is the highest plateau in the world. It covers approximately 2.5 million km2 and contains enormous glaciers, huge alpine lakes, and impressive waterfalls. There are also plentiful freshwater resources, and this region serves as the headwaters of many of Asia’s rivers (Zou et al. 2005; Wang et al. 2009). However, the distribution of water resources on the Tibetan Plateau is predominantly determined by its topographic configuration and atmospheric circulation patterns. Atmospheric circulation in Asia is currently dominated by a westerly and monsoon wind system, which causes an approximate rainfall gradient from east to west on the Tibetan Plateau (Klein et al. 2004; Shen et al. 2008). Drought is accelerating on the Tibetan Plateau, possibly because of global climate change and human activities. Treeline data indicate that extreme drought conditions have progressed into the northwest of this region over the past 300 years (Gou et al. 2006, 2007); the available soil water, which is largely determined by regional rainfall, affects plant growth, biomass accumulation, and leaf gas exchange rates. Therefore, rainfall is one of the most important factors
Funct Integr Genomics
determining the structure and function of plants and terrestrial ecosystems on the Tibetan Plateau. Additionally, because this area is one of the most sensitive to global climate change, studying the effects of rainfall on the Tibetan Plateau should be very helpful for its ecosystem management, conservation, and development. However, most previous studies have focused on macro-ecological effects, whereas very little is known about the effects of rainfall on the evolution of the ecological community and on adaptations at the physiological and molecular levels. Stress constantly challenges plants to adapt to their environments, and plants have evolved multiple mechanisms to survive harsh environments. Drought stress has drawn much attention because it is a main limiting factor for crop yield, which will be exacerbated by increasing global temperatures (Zhu 2002). Plants endure drought stress mainly by avoiding tissue dehydration, while maintaining an as-high-as-possible tissue potential or by tolerating a low tissue water potential. In fact, plants initiate an integrative mechanism during this process, involving phenotypic plasticity, alterations in biochemical enzyme activity or proteomic dynamics, and stressresponsive changes in gene expression (Chelaifa et al. 2010). Under drought stress, plants close their stomata to reduce water transpiration and photosynthesis; water loss is also minimized by reducing light absorbance by rolling leaves, a dense trichome layer, steep leaf angels, a more rigid cell wall, or small cells. In addition, plants can increase the root investment, such as increasing the root depth. Morphological plasticity and physiological adjustments are involved in the tolerance of low tissue water potential (Chaves et al. 2003). The plant hormone abscisic acid (ABA) signals the closure of stomata (Zhang et al. 2006). Biochemical processes, such as ion influx changes, post-translation modifications (e.g., protein phosphorylation, ubiquitination), and defense protein accumulation (e.g., dehydrin) also greatly contribute to tolerance drought stress (Chaves et al. 2003; Zhang et al. 2006; Nicotra and Davidson 2010). Water shortage can trigger drought-responsive gene expression, which is characterized as changes in signal transduction, enzymatic metabolism, channel protein genes, transcription factors, and other defense responses to stabilize metabolism and physiology (Fordyce 2006). In the past several years, revolutionary advances in methodology, including high-throughput sequencing and integrated tools such as stable isotopes and thermal and fluorescence imaging, have provided new insights into molecular mechanisms (Birren et al. 2011). Among them, the identification of transcriptome divergence associated with adaptive population divergence has attracted much attention because changes in gene expression can drive evolution (Ness et al. 2011). However, plant drought responses are complex biological processes that are controlled by multiple loci, and most studies to date have focused on model and crop plants to the neglect of alpine
plants that survive under harsh environmental conditions and may be highly valuable as resistance gene resources (Marioni et al. 2008; Mortazavi et al. 2008; Ekblom and Galindo 2011). Clearly, many challenges remain for the comprehensive understanding of the mechanism underlying plant tolerance to drought stress (Ozsolak and Milos 2011; Malenke et al. 2013). Stipa purpurea (Poaceae) is an endemic species that is predominantly distributed in the northwestern QinghaiXizang Plateau, Pamirs Plateau, and high mountains in central Asia (Liu et al. 2009; Yue et al. 2011). It is the dominant species in alpine arid and semi-arid meadows and contributes to stabilizing the landscape because of its outstanding tolerance to harsh environments, including cold, drought, and gale. S. purpurea also serves as forage for grassland animals and plays important roles in conserving soil and water, as a windbreak, and/or in sand fixation. Thus, investigating the genetic variation of this plant and uncovering its adaptive mechanisms for tolerating adverse environments are necessary for realizing full utilization potential in harsh, high-elevation areas. However, relatively little research has been conducted on this species. In the past 50 years, most of the work in Tibet has focused on the vegetation ecology of forest and grassland resources (Chang and Gauch 1986; Wang 1988; Zhang et al. 1988). During the 1990s, ecosystem structure and function based on vegetation types of the Tibetan Plateau received more attention (Li and Zhou 1998). The productivity of natural vegetation and the functional ecology of alpine plants have been examined recently (Luo et al. 2002). However, few studies have applied next-generation sequencing to analyze natural populations along an environmental gradient of acclimation and adaptation in Tibet. The aim of this study was to investigate putative genes and molecular functions of S. purpurea populations involved in adaptation to different field conditions and an ecological drought gradient. Therefore, we first applied high-throughput next-generation sequencing to decipher the transcriptomes of several S. purpurea populations that were collected along a rainfall gradient on the northwestern Tibetan Plateau. Our results indicated that many genes related to water metabolism and cell wall structure and synthesis differed dramatically among populations, possibly contributing to this plant’s adaptation to the drought stress gradient. It should aid us in better understanding adaptive diversity in natural populations of S. purpurea on the Tibetan Plateau.
Materials and methods Sampling and meteorological data collection A total of 125 individuals were sampled from five Tibetan populations of S. purpurea located in Laqu (population I or PI 30°37′ N, 91°44′ E), Danxiong (PII, 30°30′ N, 91°03′ E),
Funct Integr Genomics
Rikeza (PIII, 29°24′ N, 85°33′ E), Ali (PIV, 30°57′ N, 81°14′ E), and Ge’er (PV, 30°57′ N, 80°34′ E) in the summer of 2011 (Fig. 1). Fresh leaves of 20–30 individuals, spaced at least 50 m apart, were collected from each population at the flowering stage, immediately frozen in liquid nitrogen, and stored at −80 °C until RNA extraction. RNA preparation, Illumina HiSeq™ 2000 library construction, and sequencing Total RNA of 20 individuals from each population was isolated using RNAiso Plus (Qiagen, Valencia, CA, USA). RNA quality was initially identified via agarose gel electrophoresis and using a NanoDrop ND1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), as previously reported (Bai et al. 2012). Total
RNA of each sample was pooled with equivalent quantity. cDNA library was prepared according to the manufacturer’s instructions (Illumina, San Diego, CA, USA). Poly-A mRNA molecules were purified from 30 μg of total RNA (an equal mixture of RNA from all five populations) using Sera-mag Magnetic Oligo (dT) beads (Illumina) and eluted with 10 mM Tris–HCl. To avoid priming bias during the cDNA synthesis, the purified mRNA was first cleaved into small pieces using RNA fragmentation reagents (Ambion, Austin, TX, USA). The cleaved mRNA fragments were then converted into double-stranded cDNA using random hexamer primers (Illumina) and the SuperScript Double-Stranded cDNA Synthesis kit (Invitrogen, Carlsbad, CA, USA). The resulting cDNAs were purified using the QiaQuick PCR Purification Kit (Qiagen) and then subjected to end-
Fig. 1 Sampling locations of S. purpurea in Tibet (a) and the details of the sites (b)
Funct Integr Genomics
repair and phosphorylation using T4 DNA polymerase, Klenow DNA polymerase, and T4 PNK (NEB, Ipswich, MA, USA). The repaired cDNA fragments were 3′ adenylated using Klenow Exo- (Illumina), producing cDNA fragments with a single “A”’ base overhang at the 3′ ends. Illumina paired-end adapters were ligated to the ends of these 3′ adenylated cDNA fragments, and the products were purified from a 2 % TAE-agarose gel (certified low-range ultra agarose; BioRad, Hercules, CA, USA) to select for downstream enrichment. A range of cDNA fragments (200–500 bp) was excised from the gel and extracted using a QIAquick Gel Extraction kit (Qiagen). Fifteen rounds of PCR amplification were performed using primers complementary to the ends of the adapters (PCR Primer PE 1.0 and PCR Primer PE 2.0; Illumina) with Phusion DNA Polymerase. Finally, after a validation step using an Agilent 2100 Bioanalyzer and Agilent DNA 1000 chip kit (Agilent, Santa Clara, CA, USA), the cDNA library products were sequenced on a paired-end flow cell using an Illumina Genome Analyzer II at the Beijing Genomics Institute in Shenzhen, China. De novo sequence assembly and annotation The raw sequence data were transformed by base calling into raw reads, including dirty reads, such as adapters and unknown (N >5 %) or low-quality bases (raw reads showing more than 10 % of bases with a Phredscaled probability (Q) less than 20). As these data will negatively affect the subsequent bioinformatic analyses, a stringent cDNA sequencefiltering process was employed to select for clean reads and to discard dirty ones. De novo assembly was carried out using the Trinity program (Grabherr et al. 2011) and de Bruijn graphs with default settings except K-mer value (Li et al. 2010); 21-mer provided the best transcriptome assembly and formed longer contigs. Contigs without Ns were obtained by conjoining the K-mers. The reads were then mapped back to the contigs, with the paired-end reads able to detect contigs from the same transcript and the distances between these contigs. Finally, Trinity connected the contigs to provide sequences that could not be extended at either end, i.e., unigenes. Unigene length was a criterion for assembly success. S. purpurea unigenes were further processed by sequence splicing and redundancy removal, using the TGI Clustering tool to acquire non-redundant unigenes (Pertea et al. 2003). The data sets are available at the NCBI Short Read Archive (SRA) with accession number SRR825213 (SAMN02044398 and SAMN02044399). The best Blastx alignments (e-value Swiss-Prot > KEGG > COG (Ye et al. 2006). When a unigene did not align to any of the databases, ESTScan was used to predict its coding regions and direction (Iseli et al. 1999). RNA isolation and library preparation for mRNA-seq analysis To identify genes with different expression levels among samples, high-throughput mRNA sequencing (mRNASeq) analysis of each population (PI–PV) was performed separately. Briefly, total RNA was isolated form each population (20 individuals) using RNAiso Plus (Qiagen) and treated with RNase-free DNase I. Poly (A) mRNA was purified using oligo dT beads. Fragmentation buffer was used to break the mRNA into short fragments (about 200 bp), which were used as templates for first-strand cDNA synthesis using random hexamer primers. Second-strand cDNA was synthesized using RNaseH and DNA polymerase I. Then, single-end and paired-end mRNA-seq libraries were prepared following Illumina’s protocols and sequenced via Illumina HiSeq™ 2000. Analysis of mRNA-Seq reads and mapping to the reference transcriptome Clean reads were obtained as described above and were mapped to our transcriptome reference sequences using SOAPaligner/soap2 (Li et al. 2009) with no more than two mismatched bases (Mortazavi et al. 2008; Tang et al. 2009). Analysis of differential gene expression The gene expression levels were calculated using the RPKM method (Mortazavi et al. 2008), and differentially expressed genes (DEGs) were screened as previously described (Audic and Claverie 1997) to compare gene expression among the five populations. False discovery rate (FDR) was used to determine the threshold P-value in multiple tests and analysis. We used an FDR of ≤0.001 and the absolute value of log2 ratio ≥1 as thresholds to judge the significance of gene expression differences. Cluster (Eisen et al. 1998) and Java Tree View (Saldanha 2004) softwares were used to perform the hierarchical clustering analysis of gene expression patterns. Further investigation of expression profiles using semi-quantitative RT-PCR Total RNA was extracted using the RNAiso Plus (Qiagen, Valencia, CA, USA). DNA-free total RNA (5 μg) from each population (three replications; n=10 for each population) was used individually for first-strand cDNA synthesis in a 20-μL reaction volume containing 2.5 units of avian myeloblastosis virus reverse transcriptase XL (TaKaRa, Tokyo, Japan) and
Funct Integr Genomics
1 μM oligo (dT) primer. PCR was performed using 2 μL of a twofold dilution of the cDNA, 10 pmol of each oligonucleotide primer, and 1 unit of Taq polymerase (TaKaRa) in a 25-μL reaction volume. Primer sequences are listed in Supporting Information Table S1. The relative abundance of alpha tubulin was determined and used as an internal standard. The PCR cycle number was normalized to obtain a clearly visible band for the sample with the highest transcriptional level. Aliquots of the PCR products were separated on 1.5 % agarose gels and visualized using ethidium bromide. Specific amplification products of the expected sizes were observed, and their intensities were analyzed using ImageJ software (http://www.ncbi.nih.gov). ABA extraction and determination ABA quantification was performed as previously described (Jiang et al. 2012). Sample from each population was extracted overnight in 1 ml of ABA extraction buffer (methanol, containing 100 mg/L butylated hydroxyl toluene, 500 mg/L citric acid monohydrate) using the Phytodetek competitive ELISA kit (Agdia, Elkhardt, IN, USA). Samples were then centrifuged 25,000g for 30 min at 4 °C. The supernatant was dried and re-suspended with 100 μL 100 % (w/v) methanol followed by the addition of 900 μl of 1×TRIS-TBS buffer. Analysis was conducted following the Phytodetek protocol. Three independent experiments were done for each population (n=15). Stomatal observation and measurement Scanning electron microscopy was performed to evaluate the stomatal aperture using an FEI QUANTA 200 3D equipped with a large-field secondary electron detector and operated at 20 kV. Fresh leaves from each population were fixed in methanol/acetic acid/ethanol/H2O (5/5/50/40, v/v/v/v, FAA) at 4 °C for 24 h. The leaves were then dehydrated through a graded ethanol series of 70, 85, 95, and 100 % ethanol, each for 45 min. Specimens were critical-point-dried using liquid carbon dioxide and then coated with gold and palladium. With the scanning electron microscopy (SEM) results, five leaves in the same general location were used to record stomatal aperture, and the stomatal aperture of 30 stomata was measured for each sample. Histology Root tissue of 20 individuals from each population was first fixed in 2 % (v/v) glutaraldehyde in PEMT buffer (50 mmol/L PIPES, 2 mmol/L EGTA, 2 mmol/L MgSO4, and 0.05 % [v/v] Triton X-100 [pH 7.2]) at 4 °C overnight. After washing in phosphate buffer (50 mmol/L [pH 7.2]), the samples were post-fixed in 2 % (v/v) OsO4 for 2 h and then dehydrated first
through a gradient series of ethanol concentrations (30, 50, 70, 85, 95, and 100 %) and washed twice with ethanol (100 %). Tissues were finally dehydrated with acetone and embedded in low-viscosity (Spurr’s) resin (Electron Microscopy Sciences). Sections (1 μm) were cut, stained for 1 min with toluidine blue, washed twice with Milli-Q water, and viewed using light microscopy (Carl Zeiss, Thornwood, NY, USA). Proline and glycine content measurements The proline and glycine contents from the different S. purpurea populations were measured using previously reported methods (Bates et al. 1973; Greive and Grattan 1983). In brief, to determine glycine betaine, the leaf samples from different regions were weighted and extracted in 3 % sulfosalicylic acid. An aliquot of each extract (2 mL) was incubated with 2 mL of ninhydrin reagent (2.5 % w/v, ninhydrin, 60 % glacial acetic acid, 40 % 6 M phosphoric acid) and 2 mL of glacical acetic acid at 100 °C for 40 min, and the reaction was determined in an ice bath. Toluene (5 mL) was added, followed by vortexing and incubation at 23 °C for 24 h. The absorbance was measured at 520 nm, and L-Pro was used as reference standard. To determine glycine betaine, 0.5 g of dry plant material was mounted in 20 mL of deionized water and mechanically shaken for 48 h at 25 °C. The samples were then filtered, and the filtrate was diluted with (1:1) 2 N sulphuric acid (H2SO4). Aliquot (0.5 mL) was measured into a test tube and cooled in ice water for 1 h. Cold potassium iodine–iodine reagent (0.2 mL) was added, and the mixture was gently mixed with vortex mixture. The samples were stored at 4 °C for 16 h and then centrifuged at 10,000g for 15 min at 4 °C. The supernatant was carefully aspirated with 1 mL micropipette, and the peroidite crystals were dissolved in 9 mL of 1,2-dichloro ethane for 3 h. The absorbance was measured at 365 nm. The standards of glycine betaine were prepared in 1 N sulphuric acid. Three independent experiments were done for each population (n=15).
Results Adaptation of S. purpurea to an ecological drought gradient To understand how the alpine grass S. purpurea adapts to the ecological drought stress of the northwestern Tibetan Plateau, we collected samples from five populations (I–V) along an ecological gradient of annual precipitation (Fig. 1). The local soil water content gradually decreased from 14.71 % (east) to 3.13 % (west) (Supporting Information Table S2). Because many environmental stress factors affect growth and adaptation, we evaluated their contributions by PCA (Supporting Information Table S2); rainfall, water content in the soil, and
Funct Integr Genomics
air humidity were the main contributing factors, whereas other factors did not substantially affect distribution. Thus, we considered only effects of drought stress on S. purpurea adaptation. Meteorological data were provided by the National Meteorological Information Center (NMIC) of the China Meteorological Administration (CMA) (http:// cdc.cma.gov.cn/home.do). De novo assembly of Illumina short-read sequences from the S. purpurea populations To further investigate the molecular mechanism of S. purpurea adaptation to local drought stress, we first determined the gene expression profiles of the different populations. However, because S. purpurea is not a “model” plant, there was no available public genome information. Therefore, we first de novo sequenced and assembled an S. purpurea transcriptome using a mixture of RNA from the five populations. As shown in Table 1, we obtained approximately 51 million 75-bp paired-end reads after shotgun sequencing and cleaning of the raw data. Read quality was assessed using the basecalling quality scores from Illumina’s base-caller Bustard. In addition, 95.31 % of the clean reads had Phred-like quality scores at the Q20 level (a sequencing error probability of 0.01). Read assembly resulted in 180,242 contigs with mean sizes of 271 bp. Clustering generated 84,298 unigenes (48,834,004 total length, mean size 579 bp) from phases 1 and 2. The length distributions of the contigs and unigenes are Table 1 Statistical summary of the S. purpurea cDNA sequences generated using the Illumina GAII platform Raw reads Total raw reads Total clean reads Total clean nucleotides (nt) Q20 percentage (%) N percentage (%) GC percentage (%) Contig Total number Total length (nt) Mean length (nt) N50 Unigene Total number Total length (nt) Mean length (nt) N50 Total consensus sequences Distinct clusters Distinct singletons
59,189,474 51,185,972 4,606,737,480 95.31 0.00 50.86
listed in Supporting Information Table S3. The data sets are available at the NCBI Short Read Archive (SRA) with accession number SRR825213. Annotation of predicted proteins The unigene sequences were functionally annotated based on the proteins in the nr, Swiss-Prot, KEGG, and COG databases with the highest sequence similarity. Using this approach, 58,966 unigenes (69.95 %) were annotated (Supporting Information Table S4). Because of the lack of genome information for S. purpurea, 25,332 unigenes (30.05 %) could not be matched to known genes. As shown in Fig. 2, longer assembled sequences were more likely to be matched; match efficiency was 99.50 % for sequences longer than 2,000 bp but was only 57.28 % for sequences 100–500 bp in length. A total of 55.46 % of the mapped sequences had strong homology (evalue
ORIGINAL PAPER
Transcriptome analysis reveals diversified adaptation of Stipa purpurea along a drought gradient on the Tibetan Plateau Yunqiang Yang & Xiong Li & Xiangxiang Kong & Lan Ma & Xiangyang Hu & Yongping Yang
Received: 22 April 2014 / Revised: 7 November 2014 / Accepted: 17 November 2014 # Springer-Verlag Berlin Heidelberg 2014
Abstract Natural selection drives species adaptations to biotic and abiotic stresses. Species distributed along a moisture gradient, such as Stipa purpurea, a dominant grass in alpine arid and semi-arid meadows on the Tibetan Plateau, provide an opportunity to evaluate the effects of long-term adaptation to differing degrees of drought stress on gene expression. However, the genetic basis of this divergence remains largely unknown. Next-generation sequencing technologies have provided important genome-wide insights on the evolution of organisms for which genomic information is lacking. To understand how S. purpurea responds to drought stress, we selected five populations distributed along the degressive rainfall line on the northwestern Tibetan Plateau that currently present evolutionary acclimation to localized drought pressure at the physiological and biochemical levels and compared their transcriptome responses. In addition, we performed de novo assembly of the S. purpurea transcriptome using short read sequencing technology and successfully assembled 84,298 unigenes from approximately 51 million sequencing reads. We quantified gene expression level to compare their transcriptome responses using mRNA-Seq and identified differentially expressed transcripts that are involved in primary Electronic supplementary material The online version of this article (doi:10.1007/s10142-014-0419-7) contains supplementary material, which is available to authorized users. Y. Yang : X. Li : X. Kong : L. Ma : X. Hu : Y. Yang (*) Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Science, Kunming 650204, China e-mail: [email protected] Y. Yang : X. Li : X. Kong : L. Ma : X. Hu : Y. Yang The Germplasm Bank of Wild Species, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming, Yunnan 650201, China X. Li : L. Ma University of Chinese Academy of Sciences, Beijing 100049, China
and secondary plant metabolism, plant hormone synthesis, defense responses, and cell wall synthesis. Furthermore, physiological and biochemical evidence supports that abscisic acid (ABA) accumulation and cell wall strengthening derived from the differential transcripts contribute to the tolerance of S. purpurea to drought stress. The mechanisms by which S. purpurea adapts to drought stress provide new insight into how plants ecologically adapt and evolve. Keywords Transcriptome . Stipa purpurea . Gradient drought . Tibetan Plateau
Introduction With a mean altitude of more than 4000 m, the Tibetan Plateau is the highest plateau in the world. It covers approximately 2.5 million km2 and contains enormous glaciers, huge alpine lakes, and impressive waterfalls. There are also plentiful freshwater resources, and this region serves as the headwaters of many of Asia’s rivers (Zou et al. 2005; Wang et al. 2009). However, the distribution of water resources on the Tibetan Plateau is predominantly determined by its topographic configuration and atmospheric circulation patterns. Atmospheric circulation in Asia is currently dominated by a westerly and monsoon wind system, which causes an approximate rainfall gradient from east to west on the Tibetan Plateau (Klein et al. 2004; Shen et al. 2008). Drought is accelerating on the Tibetan Plateau, possibly because of global climate change and human activities. Treeline data indicate that extreme drought conditions have progressed into the northwest of this region over the past 300 years (Gou et al. 2006, 2007); the available soil water, which is largely determined by regional rainfall, affects plant growth, biomass accumulation, and leaf gas exchange rates. Therefore, rainfall is one of the most important factors
Funct Integr Genomics
determining the structure and function of plants and terrestrial ecosystems on the Tibetan Plateau. Additionally, because this area is one of the most sensitive to global climate change, studying the effects of rainfall on the Tibetan Plateau should be very helpful for its ecosystem management, conservation, and development. However, most previous studies have focused on macro-ecological effects, whereas very little is known about the effects of rainfall on the evolution of the ecological community and on adaptations at the physiological and molecular levels. Stress constantly challenges plants to adapt to their environments, and plants have evolved multiple mechanisms to survive harsh environments. Drought stress has drawn much attention because it is a main limiting factor for crop yield, which will be exacerbated by increasing global temperatures (Zhu 2002). Plants endure drought stress mainly by avoiding tissue dehydration, while maintaining an as-high-as-possible tissue potential or by tolerating a low tissue water potential. In fact, plants initiate an integrative mechanism during this process, involving phenotypic plasticity, alterations in biochemical enzyme activity or proteomic dynamics, and stressresponsive changes in gene expression (Chelaifa et al. 2010). Under drought stress, plants close their stomata to reduce water transpiration and photosynthesis; water loss is also minimized by reducing light absorbance by rolling leaves, a dense trichome layer, steep leaf angels, a more rigid cell wall, or small cells. In addition, plants can increase the root investment, such as increasing the root depth. Morphological plasticity and physiological adjustments are involved in the tolerance of low tissue water potential (Chaves et al. 2003). The plant hormone abscisic acid (ABA) signals the closure of stomata (Zhang et al. 2006). Biochemical processes, such as ion influx changes, post-translation modifications (e.g., protein phosphorylation, ubiquitination), and defense protein accumulation (e.g., dehydrin) also greatly contribute to tolerance drought stress (Chaves et al. 2003; Zhang et al. 2006; Nicotra and Davidson 2010). Water shortage can trigger drought-responsive gene expression, which is characterized as changes in signal transduction, enzymatic metabolism, channel protein genes, transcription factors, and other defense responses to stabilize metabolism and physiology (Fordyce 2006). In the past several years, revolutionary advances in methodology, including high-throughput sequencing and integrated tools such as stable isotopes and thermal and fluorescence imaging, have provided new insights into molecular mechanisms (Birren et al. 2011). Among them, the identification of transcriptome divergence associated with adaptive population divergence has attracted much attention because changes in gene expression can drive evolution (Ness et al. 2011). However, plant drought responses are complex biological processes that are controlled by multiple loci, and most studies to date have focused on model and crop plants to the neglect of alpine
plants that survive under harsh environmental conditions and may be highly valuable as resistance gene resources (Marioni et al. 2008; Mortazavi et al. 2008; Ekblom and Galindo 2011). Clearly, many challenges remain for the comprehensive understanding of the mechanism underlying plant tolerance to drought stress (Ozsolak and Milos 2011; Malenke et al. 2013). Stipa purpurea (Poaceae) is an endemic species that is predominantly distributed in the northwestern QinghaiXizang Plateau, Pamirs Plateau, and high mountains in central Asia (Liu et al. 2009; Yue et al. 2011). It is the dominant species in alpine arid and semi-arid meadows and contributes to stabilizing the landscape because of its outstanding tolerance to harsh environments, including cold, drought, and gale. S. purpurea also serves as forage for grassland animals and plays important roles in conserving soil and water, as a windbreak, and/or in sand fixation. Thus, investigating the genetic variation of this plant and uncovering its adaptive mechanisms for tolerating adverse environments are necessary for realizing full utilization potential in harsh, high-elevation areas. However, relatively little research has been conducted on this species. In the past 50 years, most of the work in Tibet has focused on the vegetation ecology of forest and grassland resources (Chang and Gauch 1986; Wang 1988; Zhang et al. 1988). During the 1990s, ecosystem structure and function based on vegetation types of the Tibetan Plateau received more attention (Li and Zhou 1998). The productivity of natural vegetation and the functional ecology of alpine plants have been examined recently (Luo et al. 2002). However, few studies have applied next-generation sequencing to analyze natural populations along an environmental gradient of acclimation and adaptation in Tibet. The aim of this study was to investigate putative genes and molecular functions of S. purpurea populations involved in adaptation to different field conditions and an ecological drought gradient. Therefore, we first applied high-throughput next-generation sequencing to decipher the transcriptomes of several S. purpurea populations that were collected along a rainfall gradient on the northwestern Tibetan Plateau. Our results indicated that many genes related to water metabolism and cell wall structure and synthesis differed dramatically among populations, possibly contributing to this plant’s adaptation to the drought stress gradient. It should aid us in better understanding adaptive diversity in natural populations of S. purpurea on the Tibetan Plateau.
Materials and methods Sampling and meteorological data collection A total of 125 individuals were sampled from five Tibetan populations of S. purpurea located in Laqu (population I or PI 30°37′ N, 91°44′ E), Danxiong (PII, 30°30′ N, 91°03′ E),
Funct Integr Genomics
Rikeza (PIII, 29°24′ N, 85°33′ E), Ali (PIV, 30°57′ N, 81°14′ E), and Ge’er (PV, 30°57′ N, 80°34′ E) in the summer of 2011 (Fig. 1). Fresh leaves of 20–30 individuals, spaced at least 50 m apart, were collected from each population at the flowering stage, immediately frozen in liquid nitrogen, and stored at −80 °C until RNA extraction. RNA preparation, Illumina HiSeq™ 2000 library construction, and sequencing Total RNA of 20 individuals from each population was isolated using RNAiso Plus (Qiagen, Valencia, CA, USA). RNA quality was initially identified via agarose gel electrophoresis and using a NanoDrop ND1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), as previously reported (Bai et al. 2012). Total
RNA of each sample was pooled with equivalent quantity. cDNA library was prepared according to the manufacturer’s instructions (Illumina, San Diego, CA, USA). Poly-A mRNA molecules were purified from 30 μg of total RNA (an equal mixture of RNA from all five populations) using Sera-mag Magnetic Oligo (dT) beads (Illumina) and eluted with 10 mM Tris–HCl. To avoid priming bias during the cDNA synthesis, the purified mRNA was first cleaved into small pieces using RNA fragmentation reagents (Ambion, Austin, TX, USA). The cleaved mRNA fragments were then converted into double-stranded cDNA using random hexamer primers (Illumina) and the SuperScript Double-Stranded cDNA Synthesis kit (Invitrogen, Carlsbad, CA, USA). The resulting cDNAs were purified using the QiaQuick PCR Purification Kit (Qiagen) and then subjected to end-
Fig. 1 Sampling locations of S. purpurea in Tibet (a) and the details of the sites (b)
Funct Integr Genomics
repair and phosphorylation using T4 DNA polymerase, Klenow DNA polymerase, and T4 PNK (NEB, Ipswich, MA, USA). The repaired cDNA fragments were 3′ adenylated using Klenow Exo- (Illumina), producing cDNA fragments with a single “A”’ base overhang at the 3′ ends. Illumina paired-end adapters were ligated to the ends of these 3′ adenylated cDNA fragments, and the products were purified from a 2 % TAE-agarose gel (certified low-range ultra agarose; BioRad, Hercules, CA, USA) to select for downstream enrichment. A range of cDNA fragments (200–500 bp) was excised from the gel and extracted using a QIAquick Gel Extraction kit (Qiagen). Fifteen rounds of PCR amplification were performed using primers complementary to the ends of the adapters (PCR Primer PE 1.0 and PCR Primer PE 2.0; Illumina) with Phusion DNA Polymerase. Finally, after a validation step using an Agilent 2100 Bioanalyzer and Agilent DNA 1000 chip kit (Agilent, Santa Clara, CA, USA), the cDNA library products were sequenced on a paired-end flow cell using an Illumina Genome Analyzer II at the Beijing Genomics Institute in Shenzhen, China. De novo sequence assembly and annotation The raw sequence data were transformed by base calling into raw reads, including dirty reads, such as adapters and unknown (N >5 %) or low-quality bases (raw reads showing more than 10 % of bases with a Phredscaled probability (Q) less than 20). As these data will negatively affect the subsequent bioinformatic analyses, a stringent cDNA sequencefiltering process was employed to select for clean reads and to discard dirty ones. De novo assembly was carried out using the Trinity program (Grabherr et al. 2011) and de Bruijn graphs with default settings except K-mer value (Li et al. 2010); 21-mer provided the best transcriptome assembly and formed longer contigs. Contigs without Ns were obtained by conjoining the K-mers. The reads were then mapped back to the contigs, with the paired-end reads able to detect contigs from the same transcript and the distances between these contigs. Finally, Trinity connected the contigs to provide sequences that could not be extended at either end, i.e., unigenes. Unigene length was a criterion for assembly success. S. purpurea unigenes were further processed by sequence splicing and redundancy removal, using the TGI Clustering tool to acquire non-redundant unigenes (Pertea et al. 2003). The data sets are available at the NCBI Short Read Archive (SRA) with accession number SRR825213 (SAMN02044398 and SAMN02044399). The best Blastx alignments (e-value Swiss-Prot > KEGG > COG (Ye et al. 2006). When a unigene did not align to any of the databases, ESTScan was used to predict its coding regions and direction (Iseli et al. 1999). RNA isolation and library preparation for mRNA-seq analysis To identify genes with different expression levels among samples, high-throughput mRNA sequencing (mRNASeq) analysis of each population (PI–PV) was performed separately. Briefly, total RNA was isolated form each population (20 individuals) using RNAiso Plus (Qiagen) and treated with RNase-free DNase I. Poly (A) mRNA was purified using oligo dT beads. Fragmentation buffer was used to break the mRNA into short fragments (about 200 bp), which were used as templates for first-strand cDNA synthesis using random hexamer primers. Second-strand cDNA was synthesized using RNaseH and DNA polymerase I. Then, single-end and paired-end mRNA-seq libraries were prepared following Illumina’s protocols and sequenced via Illumina HiSeq™ 2000. Analysis of mRNA-Seq reads and mapping to the reference transcriptome Clean reads were obtained as described above and were mapped to our transcriptome reference sequences using SOAPaligner/soap2 (Li et al. 2009) with no more than two mismatched bases (Mortazavi et al. 2008; Tang et al. 2009). Analysis of differential gene expression The gene expression levels were calculated using the RPKM method (Mortazavi et al. 2008), and differentially expressed genes (DEGs) were screened as previously described (Audic and Claverie 1997) to compare gene expression among the five populations. False discovery rate (FDR) was used to determine the threshold P-value in multiple tests and analysis. We used an FDR of ≤0.001 and the absolute value of log2 ratio ≥1 as thresholds to judge the significance of gene expression differences. Cluster (Eisen et al. 1998) and Java Tree View (Saldanha 2004) softwares were used to perform the hierarchical clustering analysis of gene expression patterns. Further investigation of expression profiles using semi-quantitative RT-PCR Total RNA was extracted using the RNAiso Plus (Qiagen, Valencia, CA, USA). DNA-free total RNA (5 μg) from each population (three replications; n=10 for each population) was used individually for first-strand cDNA synthesis in a 20-μL reaction volume containing 2.5 units of avian myeloblastosis virus reverse transcriptase XL (TaKaRa, Tokyo, Japan) and
Funct Integr Genomics
1 μM oligo (dT) primer. PCR was performed using 2 μL of a twofold dilution of the cDNA, 10 pmol of each oligonucleotide primer, and 1 unit of Taq polymerase (TaKaRa) in a 25-μL reaction volume. Primer sequences are listed in Supporting Information Table S1. The relative abundance of alpha tubulin was determined and used as an internal standard. The PCR cycle number was normalized to obtain a clearly visible band for the sample with the highest transcriptional level. Aliquots of the PCR products were separated on 1.5 % agarose gels and visualized using ethidium bromide. Specific amplification products of the expected sizes were observed, and their intensities were analyzed using ImageJ software (http://www.ncbi.nih.gov). ABA extraction and determination ABA quantification was performed as previously described (Jiang et al. 2012). Sample from each population was extracted overnight in 1 ml of ABA extraction buffer (methanol, containing 100 mg/L butylated hydroxyl toluene, 500 mg/L citric acid monohydrate) using the Phytodetek competitive ELISA kit (Agdia, Elkhardt, IN, USA). Samples were then centrifuged 25,000g for 30 min at 4 °C. The supernatant was dried and re-suspended with 100 μL 100 % (w/v) methanol followed by the addition of 900 μl of 1×TRIS-TBS buffer. Analysis was conducted following the Phytodetek protocol. Three independent experiments were done for each population (n=15). Stomatal observation and measurement Scanning electron microscopy was performed to evaluate the stomatal aperture using an FEI QUANTA 200 3D equipped with a large-field secondary electron detector and operated at 20 kV. Fresh leaves from each population were fixed in methanol/acetic acid/ethanol/H2O (5/5/50/40, v/v/v/v, FAA) at 4 °C for 24 h. The leaves were then dehydrated through a graded ethanol series of 70, 85, 95, and 100 % ethanol, each for 45 min. Specimens were critical-point-dried using liquid carbon dioxide and then coated with gold and palladium. With the scanning electron microscopy (SEM) results, five leaves in the same general location were used to record stomatal aperture, and the stomatal aperture of 30 stomata was measured for each sample. Histology Root tissue of 20 individuals from each population was first fixed in 2 % (v/v) glutaraldehyde in PEMT buffer (50 mmol/L PIPES, 2 mmol/L EGTA, 2 mmol/L MgSO4, and 0.05 % [v/v] Triton X-100 [pH 7.2]) at 4 °C overnight. After washing in phosphate buffer (50 mmol/L [pH 7.2]), the samples were post-fixed in 2 % (v/v) OsO4 for 2 h and then dehydrated first
through a gradient series of ethanol concentrations (30, 50, 70, 85, 95, and 100 %) and washed twice with ethanol (100 %). Tissues were finally dehydrated with acetone and embedded in low-viscosity (Spurr’s) resin (Electron Microscopy Sciences). Sections (1 μm) were cut, stained for 1 min with toluidine blue, washed twice with Milli-Q water, and viewed using light microscopy (Carl Zeiss, Thornwood, NY, USA). Proline and glycine content measurements The proline and glycine contents from the different S. purpurea populations were measured using previously reported methods (Bates et al. 1973; Greive and Grattan 1983). In brief, to determine glycine betaine, the leaf samples from different regions were weighted and extracted in 3 % sulfosalicylic acid. An aliquot of each extract (2 mL) was incubated with 2 mL of ninhydrin reagent (2.5 % w/v, ninhydrin, 60 % glacial acetic acid, 40 % 6 M phosphoric acid) and 2 mL of glacical acetic acid at 100 °C for 40 min, and the reaction was determined in an ice bath. Toluene (5 mL) was added, followed by vortexing and incubation at 23 °C for 24 h. The absorbance was measured at 520 nm, and L-Pro was used as reference standard. To determine glycine betaine, 0.5 g of dry plant material was mounted in 20 mL of deionized water and mechanically shaken for 48 h at 25 °C. The samples were then filtered, and the filtrate was diluted with (1:1) 2 N sulphuric acid (H2SO4). Aliquot (0.5 mL) was measured into a test tube and cooled in ice water for 1 h. Cold potassium iodine–iodine reagent (0.2 mL) was added, and the mixture was gently mixed with vortex mixture. The samples were stored at 4 °C for 16 h and then centrifuged at 10,000g for 15 min at 4 °C. The supernatant was carefully aspirated with 1 mL micropipette, and the peroidite crystals were dissolved in 9 mL of 1,2-dichloro ethane for 3 h. The absorbance was measured at 365 nm. The standards of glycine betaine were prepared in 1 N sulphuric acid. Three independent experiments were done for each population (n=15).
Results Adaptation of S. purpurea to an ecological drought gradient To understand how the alpine grass S. purpurea adapts to the ecological drought stress of the northwestern Tibetan Plateau, we collected samples from five populations (I–V) along an ecological gradient of annual precipitation (Fig. 1). The local soil water content gradually decreased from 14.71 % (east) to 3.13 % (west) (Supporting Information Table S2). Because many environmental stress factors affect growth and adaptation, we evaluated their contributions by PCA (Supporting Information Table S2); rainfall, water content in the soil, and
Funct Integr Genomics
air humidity were the main contributing factors, whereas other factors did not substantially affect distribution. Thus, we considered only effects of drought stress on S. purpurea adaptation. Meteorological data were provided by the National Meteorological Information Center (NMIC) of the China Meteorological Administration (CMA) (http:// cdc.cma.gov.cn/home.do). De novo assembly of Illumina short-read sequences from the S. purpurea populations To further investigate the molecular mechanism of S. purpurea adaptation to local drought stress, we first determined the gene expression profiles of the different populations. However, because S. purpurea is not a “model” plant, there was no available public genome information. Therefore, we first de novo sequenced and assembled an S. purpurea transcriptome using a mixture of RNA from the five populations. As shown in Table 1, we obtained approximately 51 million 75-bp paired-end reads after shotgun sequencing and cleaning of the raw data. Read quality was assessed using the basecalling quality scores from Illumina’s base-caller Bustard. In addition, 95.31 % of the clean reads had Phred-like quality scores at the Q20 level (a sequencing error probability of 0.01). Read assembly resulted in 180,242 contigs with mean sizes of 271 bp. Clustering generated 84,298 unigenes (48,834,004 total length, mean size 579 bp) from phases 1 and 2. The length distributions of the contigs and unigenes are Table 1 Statistical summary of the S. purpurea cDNA sequences generated using the Illumina GAII platform Raw reads Total raw reads Total clean reads Total clean nucleotides (nt) Q20 percentage (%) N percentage (%) GC percentage (%) Contig Total number Total length (nt) Mean length (nt) N50 Unigene Total number Total length (nt) Mean length (nt) N50 Total consensus sequences Distinct clusters Distinct singletons
59,189,474 51,185,972 4,606,737,480 95.31 0.00 50.86
listed in Supporting Information Table S3. The data sets are available at the NCBI Short Read Archive (SRA) with accession number SRR825213. Annotation of predicted proteins The unigene sequences were functionally annotated based on the proteins in the nr, Swiss-Prot, KEGG, and COG databases with the highest sequence similarity. Using this approach, 58,966 unigenes (69.95 %) were annotated (Supporting Information Table S4). Because of the lack of genome information for S. purpurea, 25,332 unigenes (30.05 %) could not be matched to known genes. As shown in Fig. 2, longer assembled sequences were more likely to be matched; match efficiency was 99.50 % for sequences longer than 2,000 bp but was only 57.28 % for sequences 100–500 bp in length. A total of 55.46 % of the mapped sequences had strong homology (evalue
