Transcriptional Responses of Chilean Quinoa (Chenopodium quinoa Willd.) Under Water Deficit Conditions Uncovers ABA-Independent Expression Patterns

HIGHLIGHTS

• R49 genotype displayed best performance on selected physiological parameters and highest tolerance to drought.

• R49 drought over-represented transcripts has exhibited 19% of genes (306 contigs) that presented no homology to published databases.

• Expression pattern for canonical responses to drought such as ABA biosynthesis and other genes induced in response to drought were assessed by qPCR.

Global freshwater shortage is one of the biggest challenges of our time, often associated to misuse, increased consumption demands and the effects of climate change, paralleled with the desertification of vast areas. Chenopodium quinoa (Willd.) represents a very promising species, due to both nutritional content and cultivation under water constraint. We characterized drought tolerance of three Chilean genotypes and selected Genotype R49 (Salares ecotype) based upon Relative Water Content (RWC), Electrolyte Leakage (EL) and maximum efficiency of photosystem II (F_v/F_m) after drought treatment, when compared to another two genotypes. Exploratory RNA-Seq of R49 was generated by Illumina paired-ends method comparing drought and control irrigation conditions. We obtained 104.8 million reads, with 54 million reads for control condition and 51 million reads for drought condition. Reads were assembled in 150,952 contigs, were 31,523 contigs have a reading frame of at least 300 nucleotides (100 aminoacids). BLAST2GO annotation showed a 15% of genes without homology to NCBI proteins, but increased to 19% (306 contigs) when focused into drought-induced genes. Expression pattern for canonical drought responses such as ABA biosynthesis and other genes induced were assessed by qPCR, suggesting novelty of R49 drought responses.

Introduction

Plants are sessile organisms that need to respond to a wide array of abiotic and biotic stresses. This condition confers strong selective pressures on their local adaptation to different environments. The selective pressure at the abiotic level implies stress responses, that includes numerous and complex environmental conditions, such as light intensity, temperature, salinity and drought. Therefore, understanding abiotic stress responses at the physiological and genomic level is a relevant issue to provide an essential foundation for future breeding and genetic engineering efforts. Indeed, roughly 75% of the world's freshwater supplies are utilized in agriculture, and the increasing climatic variability and the demographic pressures have led to ecosystem degradation and have exacerbated the vulnerability to drought and other abiotic stress factors (Food Agriculture Organization of the United Nations, 2011). Consequently, increasing crop productivity under conditions of limiting water availability is of major importance (Avramova et al., 2015). Research on drought responses in Arabidopsis, maize, tomato and rice among others plants, have determined that a large number of genes as well as signal transduction pathways are involved in drought responses (Shinozaki et al., 2003; Shinozaki and Yamaguchi-Shinozaki, 2007). The identification of genes regulated by drought conditions has a high significance, because it provides a comprehensive understanding of the transcriptional responses and the identification of stress responsive promoters and cis-elements.

Quinoa (Chenopodium quinoa Willd.) is an Andean native crop that belongs to the Amaranthaceae family. It is an allotetraploid plant (2n = 4× = 36) with a genome size estimated of 967 Mbp (Mujica and Jacobsen, 2006; Stevens et al., 2006) and sequenced very recently with a total assembly size of 1.39 gigabases (Gb) (Jarvis et al., 2017), which shows disomic inheritance for most qualitative traits (Simmons, 1971; Risi and Galwey, 1984; Ward, 2000; Maughan et al., 2004). Quinoa was domesticated and has been cultivated in the Andes for the last 7,000 years before present (BP). It's diversity is comprised by five major ecotypes related to different sub-center origins that include: Highlands (Peru and Bolivia), Inter-Andean valleys (Bolivia, Colombia, Ecuador and Peru), Salares (Bolivia, Chile and Argentina), Yungas (Bolivia) and Lowlands (Chile) (Risi and Galwey, 1989a,b; Bertero et al., 2004; Zurita-Silva et al., 2014; Bazile et al., 2015). Among nutritional characteristics it is know that quinoa's seeds have an exceptional balance between oil (4–9%), protein (averaging 16%, with high nutritional relevance due to optimal balance of essential amino acid content) and carbohydrates (64%) (Bhargava et al., 2006; Vega-Gálvez et al., 2010). Moreover, its high starch content (51–61%) enabling flour production (Mastebroek et al., 2000; Repo-Carrasco et al., 2003; Stikic et al., 2012) but with the advantages of gluten absence, which has allowed the development of foods for consumers with celiac disease (i.e., people allergic to gluten) (Jacobsen, 2003). Additionally, quinoa is a good source of vitamins, oil with high linoleate and linolenate content (55–66% of the lipid fraction), natural antioxidants such as α- and γ-tocopherol, and a wide range of minerals (Repo-Carrasco et al., 2003; Vega-Gálvez et al., 2010; Fuentes and Bhargava, 2011; Stikic et al., 2012; Ruiz et al., 2014). Interestingly, quinoa consumption may lead to comparatively lower weight gain, and improved lipid profile and potential antioxidant effects, physiological outcomes that have been linked to bioactive compounds, such as saponins, quinoa proteins, polyphenolic compounds and 20-hydroxyecdysone by yet unknown mechanisms (Simnadis et al., 2015). Considering the attributes and potential to contribute to food security worldwide, the draft genome sequence of an inbred line has been recently published, comprising a free-access Quinoa Genome DataBase (QGDB), which will provide insights into the mechanisms underlying agronomically important traits of quinoa (Yasui et al., 2016).

Quinoa is an interesting abiotic stress tolerant crop that should be adopted as a model because has a good tolerance to high salinity, boron, light intensity and drought (Orsini et al., 2011; Ruiz-Carrasco et al., 2011; Zurita-Silva et al., 2014; Razzaghi et al., 2015; Ruiz et al., 2016). In particular quinoa has a good tolerance to water shortage that is due to its intrinsic lower water requirement, the capability to regain its original level of photosynthesis after a drought period, and both slow growth and smaller leaf area during acclimation (Galwey, 1989; Jensen et al., 2000; Jacobsen et al., 2003, 2009; Sun et al., 2014). There are others factors than can be considered for this tolerance i.e., high instantaneous photosynthetic efficiency in drought conditions (Winkel et al., 2002; Bosque Sanchez et al., 2003); leaf shedding (Jensen et al., 2000); its higher root branching and foraging capacity of root system (Alvarez-Flores et al., 2014), and the presence of leaf vesicles containing calcium oxalate, which could reduce transpiration (Jensen et al., 2000; Siener et al., 2006). It's noteworthy that this tolerance to water deficit has allowed reaching harvests with only 75 mm of rainfall (Martínez et al., 2009). In this work we characterized three Chilean genotypes (Martínez et al., 2015), one corresponding to Salares and two from Lowlands native ecotypes, which frequently experience water deficit during growth season, in terms of drought tolerance. Later, for the most tolerant genotype, a RNA-seq analysis of the transcriptome was performed to explore candidate transcripts under contrasting water conditions i.e., deficit and normal irrigation. A deep bioinformatics analysis is presented that assessed a putative signal transduction pathway involved in quinoa responses to drought. To validate the assembly of the contigs, primer pairs were designed to amplify the sequences of 15 selected transcripts, related to drought response and ABA biosynthesis and to study the expression level of target genes through Real Time -quantitative PCR (RT-qPCR).

Materials and Methods

Plant Material

Three Chilean genotypes of Chenopodium quinoa Willd. were used in this study: R49, from Colchane, Tarapacá Region, 3,850 m altitude (Salares ecotype); PRJ, from Cahuil, O'Higgins Region, 39 m altitude and BO78 from Collipulli, Araucanía Region, 243 m altitude (the last two coastal/lowland ecotypes; Table 1). The seeds were kindly provided by the National Seed Bank of Chile managed by Instituto de Investigaciones Agropecuarias INIA Intihuasi (Vicuña, Chile).

TABLE 1

Table 1. Chilean genotypes studied: ecotype and geographical origin.

Plant Growth and Stress Conditions

Experiment consisted of 12 pots each containing two plants for each genotype that were grown in 1L pots with a soil-sand mixture (1:1) with irrigation every 2 days and temperatures that fluctuated between 12 and 32°C during 4 weeks at greenhouse conditions. To select the genotypes tolerant to drought, at the fourth week since germination half of the plants from each genotype were maintained as control in the same irrigation regimen, and the other half were deprived of water, and distributed in six blocks that were set up according to a randomized block design. The experiments were performed in La Serena, Coquimbo Region, Chile (29°54′; 71°15′) during summer (January).

To determinate the most drought tolerant genotype, we analyzed the following physiological parameters: (1) Relative Water Content (RWC) as described previously by Barrs and Weatherley (1962), (2) Electrolyte Leakage (EL) as described by Pinhero et al. (1997), and F_v/F_m determination (Woo et al., 2008). The experiment was done using 1-month-old plants under a water shortage condition (dry season). The starting of the experiment was ceasing the plants irrigation and the measurement of the physiological parameters previously mentioned. For F_v/F_m determination, three different leaves by treatment were measured from dark-adapted plants during 30 min, and then were subjected to an initial saturating pulse of 900 μmol photons m⁻² s⁻¹, followed by a 40″ delay in darkness and subsequently 10' of actinic illumination with saturating flashes at 20″ intervals, using a foliar chamber 6400-40 coupled to LI-6400 fluorometer (LICOR), and were calculated following reported methods (Woo et al., 2008).

For the transcriptome approach, a second experiment was performed using the most drought tolerant genotype. The growth conditions were similar to the previous experiment and plant tissues were collected at the following stages: (a) Stage 1 (D1), where relative water content (RWC) was near 80% (similar to control) and soil potential was significant different between control and drought (−0.6 MPa and −1.3 MPa, respectively); (b) Stage 2 (D2), where leaf RWC was near 50–60% and the soil potential was −1.3 MPa; and (c) Stage 3 (D3), where leaf RWC was near 30% and soil potential was −2.7 MPa. For the control samples, the same stages were used (C1, C2, C3), in which the analyzed parameters were similar to stage 1 control (80% RWC and −0.6 MPa soil potential) (Table S1). The collected tissues (whole plants) were frozen in liquid nitrogen and then placed at −80°C until further use.

RNA Extraction

Total RNA was obtained from plant tissues using the Trizol reagent (Invitrogen Corp., Carlsbad, CA, USA) following the manufacturer's protocol. We extracted the RNA individually from three different plants for each stage (D1, D2, and D3; C1, C2, and C3), and tissues were separated between root and shoot. After quality analysis (Meisel et al., 2005) we composed two different pools (drought and control) using the same quantity from each plant tissue (root and shoot, 1.6 μg each) and sent to Macrogen (http://www.macrogen.com) for Illumina 100 bp paired-ends sequencing procedure.

Library Construction, Deep Sequencing and De novo Transcriptome Assembly

Library construction and deep sequencing for each sample were performed at Macrogen (Inc. Seoul, South Korea) using Solexa HiSeq2000 platform with the previous construction of a Truseq mRNA library for paired end application according to Macrogen's protocol with an insertion length of 550 bp. The sequence reads were quality trimmed and assembled using the CLC Genome Workbench version 4.8 (CLC Bio: CLC genomics workbench [http://www.clcbio.com]). Trimming was done following parameters: Q ≥ 20; no more than 2 ambiguities; final read length ≥ 50 bp. Reads assembly was done from a pooling of all the paired end short-read data (hybrid assembly) using the following parameters: similarity = 0.95; length fraction = 0.7; insertion/deletion cost = 3; mismatch cost = 3; automatic bubble and word size; minimum contig length of 200 bp to avoid singlets. Paired end range distance was 217 to 442 bp for the control sample and 211 to 380 bp for the drought treatment sample. To compare contigs assembly we used both CLC and Oases 0.1.22 (Schulz et al., 2012) [http://www.ebi.ac.uk/~zerbino/oases]. Contigs were scanned for full length coding sequence (full length CDS) using the tool GETORF from EMBOSS version 6.3.1 (Rice et al., 2000), looking for open reading frame regions between START and STOP codons and a minimum of 300 nucleotides. Considering that drought related proteins described in other models have more than 100 aa, and incomplete genes can't be resolved without a deep genome exploration that is out of the scope of this work, we chose this cutoff to avoid the presence of incomplete genes. Results were filtered through Perl script to obtain the longest predicted CDS by contig.

Exploratory Data Analysis

Contigs with predicted full length CDS obtained from the de novo hybrid assembly were used as a reference set of transcripts for RNA-seq analysis. Short-read sequence data from control and drought samples were separately mapped against the reference set of assembled transcripts using the CLC Genome Workbench RNA-seq function using the following parameters: similarity = 0.95; length fraction = 0.7; maximum mismatches = 2; unspecific match limit = 10. Paired reads were counted as two, and paired end distances were set as described previously for the assembly.

For exploratory expression analysis, we selected contigs with at least 5 total reads mapped to each hybrid contig to increase the results confidence (Tittarelli et al., 2009). The relative “over/under-represented” gene levels were defined as the number of reads mapped uniquely to the gene. The representation levels were compared using a Z-Test (Kal et al., 1999) with the control sample as reference. This test compares counts by considering the proportions that make up the total sum of counts in each sample, correcting the data for sample size. For visual inspection, original over/under-represented values were transformed by Log₁₀ method and then normalized by the Quantile method that was the best to fits the results (Bolstad et al., 2003).

Functional Annotation

Functional annotation was performed on assembled contigs with predicted full length CDS by BLAST2GO software (Conesa and Gotz, 2008). We used as input BLASTX results of our contigs against NCBI RefSeqProt database (Pruitt et al., 2007; Mestanza et al., 2015) and Arabidopsis TAIR10 database (The Arabidopsis Information Resource, Genome Release 10 [http://www.arabidopsis.org]) with an e-value cut-off of 1e⁻⁶. Also, it was performed an INTERPROSCAN analysis (Hunter et al., 2009) with BLAST2GO default parameters. We used BLAST2GO program defaults in all mapping and annotation steps and the False Discovery Rate (FDR) cut-off was set to 0.05% probability level. The data from INTERPROSCAN terms, enzyme classification codes (EC), and metabolic pathways (KEGG, Kyoto Encyclopedia of Genes and Genomes) were merged with GO terms for a wide functional range cover in annotation (Conesa and Gotz, 2008).

Quantitative Real-Time PCR Analysis

For each sample, 10 μg total RNA was treated with DNAse RNAse-free (Fermentas), 5 μg of which was reverse transcribed in a 20 μL volume using Affinity Script Multiple Temperature cDNA Synthesis Kit (Agilent) primed with oligo dT. The resulting cDNA was diluted to 200 μL with distillate water. Gene-specific primers were designed to span the selected genes using Primer3 software (http://frodo.wi.mit.edu/primer3/). qPCR was carried out on 1 μL diluted cDNA by triplicate using the MaxPro3000P Stratagene Sequence Detection System, Brilliant III Ultra Fast SYBR Green QPCR master mix (Agilent) and primers at a final concentration between 250 and 450 nM. The primers list is shown in (Table S2). qPCR analysis was performed and data (±S.D.) represents three biological replicates following previous reports (Ruiz-Carrasco et al., 2011; Mestanza et al., 2015).

Statistical Analysis

A two-way analysis of variance (ANOVA), followed by Tukey's post-test, was performed to evaluate the drought effects on relative water content, electrolyte leakage and F_v/F_m determination to compare genotypes among each other, after testing for normality and homogeneity of variances using the Shapiro-Wilkes test and variation coefficient, respectively.

Results

R49 Exhibits the Most Drought Tolerant Phenotype

Based on the altitudinal and latitudinal distribution of quinoa in Chile we analyzed three genotypes representing the biogeographical areas where this staple crop is actually grown. The genotype R49 is representative of the Chilean highlands (Salares ecotype), where 200 mm/year of rainfall are concentrated in only 1 month during the summer rainy season at an altitude of 3,800 m. The PRJ genotype is grown in central coastal Chile with a rainfall of 600 mm/year that is also concentrated during the fall season, and the genotype BO78 is from the southern Chile were the rainfall reaches 1,300 mm and the rain period occur throughout the whole year, the last two genotypes corresponding to coastal/lowland ecotype.

Figure 1 shows the genotype responses for relative water content, electrolytic leakage and maximum efficiency of photosystem II (F_v/F_m) parameter. The analysis of the relative water content (Figure 1A) showed a gradual decrease for all genotypes from 8 days from the treatment onward thus decreasing their RWC below 30%; the genotype R49 showed a slight difference compared to the others genotypes after 21 days of drought. We could observe a significant difference [F_(1.63) = 129.01; p < 0.0001] since 18 days post treatment compared to PRJ. On the other hand, BO78 showed an intermediate behavior between both genotypes but not significantly different. When measuring electrolytic leakage (Figure 1B), a significant increase [F_(1.63) = 34.37; p < 0.0001] for the central and south genotypes was observed at 18 days, increasing over 70% EL; nevertheless, the Salares genotype did not display major differences through the whole experiment similar to control of all genotypes. The parameter of maximum efficiency of photosystem II (F_v/F_m; Figure 1C) significantly decreased [F_(1.63) = 6.47; p < 0.0001] as response to drought treatment for genotypes PRJ and BO78 only from 18 days of drought treatment, in contrast to a stable F_v/F_m ratio in genotype R49 throughout the whole experiment, showing that photosynthetic machinery in R49 leaves remained functional despite prolonged water deficit condition.

FIGURE 1

Figure 1. Physiological performance of quinoa genotypes in days under drought conditions. (A) Percentage of Relative water content (RWC). (B) Percentage of Electrolytic leakage (EL). (C) Maximum efficiency of photosystem II (F_v/F_m). X-axis represents days under drought. RWC was measured based in the plant total weight that corresponded to water; EL represents the percentage of ions that were released compared to total amount present in the plant; F_v/F_m represents the capacity for photon energy absorbed by photosystem II (PSII) to be utilized in photochemistry under dark- and light-adapted conditions. Values and Bars represent averages and standard deviation (n = 3); ** and different letters denote significant differences (p < 0.01 and p < 0.05 respectively).

RNA-Sequencing, Reference Transcriptome Assembly and Gene Space of Quinoa under Drought Conditions

Two different samples were sequenced to obtain the R49 quinoa transcriptome. One of the samples corresponded to plants with the drought stress treatment and the other one to control plants as described in materials and methods. After library construction, Illumina paired -end sequencing was performed. 54 million and 51 million reads were obtained from the control and drought conditions respectively. Both reads from control and drought conditions were used to build a reference transcriptome due to the lack of an available quinoa genome at that stage.

The sequence assembly was performed with the CLC Genome Workbench software (version 4.8). The 104.8 millions of reads were assembled into 150,952 contigs, with an average length of 538 bp, their size distribution is showed in (Figure S1). Out of the assembled contigs, 31,523 with predicted full length CDS were used for further analysis as our reference transcriptome. The size distribution of this subset of sequences is showed in Figure 2A and the list of sequences is available in Additional File 1. Reads of both conditions (drought and control) were individually mapped back to our reference transcriptome and results were used in an RNA-seq assay to do a qualitative analysis (e.g., present/absent of a gene) (Mortazavi et al., 2008). Gene data is available in Additional File 2.

FIGURE 2

Figure 2. Identity analysis of quinoa genes. (A) Contig size distribution (bp). (B) Blast hit sequence similarity distribution. (C) Blast hit sequence E-value distribution. Red line represents the accumulation rate from low value to high whereas green line represents the accumulation ratio from high value to low.

Functional Annotation of Reference Transcriptome

To obtain a functional annotation, the sequences with predicted full length CDS (n = 31,523) were aligned by BLAST against the NCBI NR database and Arabidopsis TAIR10 database using a cut-off E-value of 10⁻⁶. The number of quinoa genes that exhibited homology with the sequences described in the NCBI NR database was 29,334 (93.1%); on the other hand, 2,189 genes (6.9%) did not match the database. These results were different to others results previously reported (Riggins et al., 2010; He et al., 2012; Huang et al., 2012; Zhang et al., 2012; Raney et al., 2014), where the genes without match represented more than 40%. In this case contigs with predicted full length CDS where used but the cited published works included the whole set of contigs. We also cannot rule out that unknown sequences may be novel genes in quinoa, which might be possibly derived from chimerical sequences (assembly's errors), or could be non-conserved regions of proteins. The similarity distribution showed that 41% of the query sequences have a similarity higher than 80%, while 91% of the hits have a similarity higher than 50%, (Figure 2B). The E-value distribution of the top hits in the NR database showed that 59.4% of the mapped sequences have strong homology (<1.0E⁻⁶³), whereas 28.6% of the homolog sequences ranged between 1.0E⁻⁷ and 1.0E⁻⁶³ and only 7% was >1.0E⁻³, which included sequences with no hit (Figure 2C). Those quinoa genes displayed a significant range of identity to sustain that a putative gene function was correlated to a functional gene homology. Indeed, with the very recent availability of a reference quinoa genome (Yasui et al., 2016), out of 31,523 contigs with predicted full CDS, 30,445 mapped to 90,438 genome predicted genes (e < 1e⁻¹⁰; score > = 500; coverage >= 50%), meaning that our built-in reference transcriptome covers a 40% of the genome transcripts (data not shown).

Annotation and determination of gene ontology of transcripts with predicted full length CDS was done by BLAST2GO platform. BLASTX, Interproscan, KEGG and Classification Codes enzymes (EC) were combined to perform a classification with a bigger GO coverage. Additional File 3 contains all the results associated to functional annotation. The results of the GO annotation were used to compare GO categories that changed in response to drought in the expression qualitative subset (2,456 contigs). Figure 3 presents the results for the categories: biological process and molecular function, level 3, and cellular component, level 8. As expected, it was possible to identify a higher number of genes over-represented by drought (85) in the biological process category “response to stress” (Figure 3A) and also 60 genes that were down-represented. In the category “cellular response to stimulus” we found the induction of 18 genes and in the category “response to abiotic stimulus” 57 genes were over-represented and 12 down-represented by drought treatment. Analysis of molecular functions that were affected by drought (Figure 3B), showed a majority of genes related to “organic cyclic compound binding and heterocyclic compound binding,” where 214 of them were over-represented and 105 were down-represented, whereas in “small molecule binding” 123 genes and 64 genes respectively. On the other hand, in the most general category “protein binding” 95 genes were over-represented and 80 genes were down-represented. Also “transferase activity” (with 143 and 113 respectively), “transcription factor activity, sequence-specific DNA binding” (18 and 16), and “hydrolase activity” (133 and 118 genes) in response to drought treatment.

FIGURE 3

Figure 3. GO classification of genes with differential expression under drought conditions. Functional categories (A) Biological process; (B) Molecular function; (C) Cellular component. Red bars: over-represented genes; Blue bars: down-represented genes.

Regarding the “cellular compartment” (Figure 3C), level 8 of Gene Ontology was used to gain insights, and the results showed that the number of genes related to “plastids” and “mitochondrion” that modified their representation was remarkable: 44 genes were over-represented and 21 genes were down-represented for “plastids,” whereas 35 genes and 17 genes were modified for “mitochondrion.” This was followed by the category “nucleus” with 48 genes and 27 genes respectively. Also we found 2 categories that were only present in drought subset: “vacuole” with 19 genes and “cytosol” with 17 genes.

Gene Ontology (GO) Enrichment Analysis

With each selected gene set we performed a gene ontology enrichment analysis with focus on biological processes terms and found that only GO terms “fruit ripening” and “reproduction” were enriched in drought over-represented genes. Furthermore, “pollen pistil interaction,” “abscission” and “carbohydrate metabolic process” were enriched in drought down-represented genes (Table 2). Indeed, we used 4 X fold change as a cutoff taking in consideration that enrichment analysis did not produce any enriched biological process at 2 X of fold of change.

TABLE 2

Table 2. Enriched Biological Processes in drought over-represented gene set or down-represented gene set.

RT-qPCR Analysis of Differentially-Expressed Genes

To confirm both the assembly and the exploratory approach to select target genes, expression analysis by qPCR amplification was performed (target genes and their representation are included in Additional File 1). We also performed Sanger sequencing of all genes used in the qPCR analysis of this work finding a perfect match between the assembly and the sequenced fragment in 90% of the cases (see Additional Files 4, 5). On the other hand, when mapping the predicted contigs with full CDS against the Quinoa genome (Yasui et al., 2016) we found that 90% had a structure according to the genome (e < 1e⁻¹⁰; score > = 500; coverage > = 50%). The gene selection was made based on two criteria. First, genes that were induced under drought conditions in other plant models were identified through literature mining. In addition, we attempted to reconstruct a canonical pathway associated with drought stress, the Abscisic Acid (ABA) pathway. The ABA biosynthesis pathway was used as a reference, corresponding to the described by Seo and Koshiba (2011). This pathway consists of 5 genes, three of which are located in the plastids: ABA1, ABA4 and NCED3 and two in the cytosol: ABA2 and ABA3. In addition to ABA biosynthesis pathway genes, two important genes involved in ABA transport were included, ABCG25 and ABCG40. Others genes involved in stress response were selected based in the changes of representation reads detected by RNA-Seq in quinoa: CqHSP20 (putative chaperones hsp20-protein superfamily), CqCAP160 (cold acclimation protein 160), CqLEA (late embryogenesis abundant protein family protein), CqAP2/ERF (integrase-type DNA-binding protein superfamily), CqPP2C (protein phosphatase protein family 2c), CqHSP83 (chaperone protein, protein family HTPG), and CqP5CS (delta 1-pyrroline-5-carboxylate synthase 2). To address the RNA-Seq changes, expression analysis by qPCR of selected 15 unigenes were performed. The qPCR analysis showed a similar pattern to the in silico analysis (Figure 4). We determined that CqNCED3a and CqNCDE3b were the only genes involved in the ABA biosynthesis pathway that were up-regulated by drought in quinoa. For CqABA1 and CqABA3, we determined slight down-regulation by drought (≤2 Fold Change), whereas CqABA2, and CqABA4 were down-regulated (>2 Fold Change). Similar slight down-regulation were determined for CqABCG25 and CqABCG40, involved in ABA transport (Figure 4A). In addition, the genes CqHSP20, CqLEA, CqCAP160, CqAP2/ERF, CqPP2C, CqHSP83, and CqP5CS were up-regulated at variable magnitudes in response to drought conditions as Figure 4B shown, which is in agreement with previous reports in different plant models (Kaye et al., 1998; Wang et al., 2004; Ali-Benali et al., 2005; Nakano et al., 2006; Swindell et al., 2007; Umezawa et al., 2009; Sharma and Verslues, 2010; Merewitz et al., 2011; Candat et al., 2014). Interestingly, CqHSP20 and CqLEA experienced a shift over 140-fold expression level, whereas other group of genes (CqCAP160, CqAP2/ERF, CqPP2C, CqHSP83, and CqP5CS) the expression level change was between 24 and 2-fold respectively.

FIGURE 4

Figure 4. Gene expression levels measured by qPCR. (A) ABA transport and biosynthesis genes expression levels. (B) Genes that respond to drought. The expression levels are relative to the normalizer gene (Pre-mRNA splicing PRP18-interacting factor) that was identified from the in silico expression analysis. Bars represent the standard deviation of three replicates.

Discussion

Quinoa or quinua (Chenopodium quinoa Willd.) is vastly distributed throughout the Andean region covering a range of 12,000 km from Colombia to Chile, also defining five ecotypes: inter- Andean valleys, Highlands, Yungas, Salares (salt flats), and coastal/lowlands. Due to the geographical diversity of our country, several genotypes can be found in the only two ecotypes present: Salares, which is distributed in the Tarapacá and Antofagasta regions (18–25°S), with elevations over 3,000 m high and precipitation fluctuating between 100 and 200 mm per year falling during the southern hemisphere summer, whereas coastal/lowlands, the latter being the only temperate latitude ecotype is distributed in central Chile and more southern latitudes (43°S) (O'Higgins to Lakes' regions) and are rain-fed with variable altitudes between sea level and 1,000 m height. A remarkable difference is that compared to the extremely dry conditions where the “Salares” quinoa is grown in northern Chile, rainfall in the central and southern zones of Chile occurs during the southern hemisphere winter (June–August), with rainfall fluctuations between 500 and 2,000 mm per year. This rainfall increases steadily across 34–40°S (Martínez et al., 2015). Among this diversity it is possible to distinguish two groups through genetic distance, most probably due to the lack of seed exchanges: “coastal/lowlands” and “Salares” (Fuentes et al., 2008; Martínez et al., 2015). Therefore, it was necessary to perform a selection among the Chilean genotypes that could exhibit distinctive tolerance to drought. Three representative Chilean genotypes were assessed to select the most drought tolerant: R49 (Salares), PRJ (coastal/lowlands), and BO78 (coastal/lowlands). These genotypes represent three bio geographical areas or growing zones (north, central and south) previously described (Fuentes et al., 2008, 2012). After the analysis of relative water content, electrolytic leakage and maximum efficiency of photosystem II, we determined that the R49 genotype was the one with the best performance on physiological parameters selected and the highest tolerance to drought, considering that decline in photosynthetic parameters occurred concurrently with the appearance of physical symptoms of drought stress, including water loss and membrane stability among others.

The quinoa R49 transcriptome was sequenced by the Illumina pair-ends method, using total RNA extracted from quinoa plants under stress and from total RNA extracted out of control plants. We obtained 104.8 millions of reads, which were assembled into 150,952 contigs with an average size of 538 bp, where 18,124 contigs had a higher or equal length of 1,000 bp and 132,826 contigs were less than 1,000 bp in length. The high percentage of contigs less than 1 Kb is expected on this type of approach (de novo assembly) were there were no reference transcriptome to date of experiments and also due to the polyploidy of this specie (Crawford et al., 2010; Mizrachi et al., 2010; Zhang et al., 2012). Taking in consideration that a high number of these contigs could be fragments of genes, we decided to select those contigs that has been predicted as full length CDS (n = 31,523 contigs). To determine the identity of those predicted full length CDS quinoa contigs, BLAST, Interproscan, KEGG and EC tools were used. The results indicated that 6.9% of the analyzed genes had no identity with any other gene in the public databases used in this study. This is a small number of unknown genes compared to other similar reports but the difference may lay in the fact that we used only a subset of contigs with predicted full length CDS to ensure more confidence avoiding assembly problems related to polyploidy. Among the contigs with known identity, we observed a small percentage having a low homology; therefore we assumed that the function of the genes was indeed the one predicted for the orthologous gene.

We indeed found six processes that could be reduced to two main biological processes: “fruit ripening” and “reproduction” under a 4 X fold change criterion and FDR p < 0.05, 2,456 differentially represented genes were identified, from which 1,579 were over-represented and 877 were down-represented under drought conditions. By analyzing the genes identity, we found 76 genes with unknown identity that were down-represented under drought, equivalent to 9% of the total of genes that were down represented, whereas 306 genes with unknown identity were over-represented under this stress, equivalent to 19% of the total of genes over-represented in drought. When we searched for the number of drought-induced genes in other models such Arabidopsis (Seki et al., 2002), the number of genes with known identity was one fifth to the ones induced in quinoa under the same stress (277 vs. 1,273 respectively). By taking a sample of the genes population, it was expected that the percentage of genes without identity would be conserved from what is observed in the entire population, in our case around 6.9%. We finally determined that among drought over-represented transcripts, genes with unknown identity were 19%, well above the expected 6.9%. These results allow us hypothesized that the drought response in quinoa might present several unknown paths. Recently, Raney et al. (2014) published the transcriptome sequencing of two quinoa genotypes in response to drought. These authors obtained a lower number of genes in comparison to those reported in our work, i.e., 462 differentially expressed contigs identified upon water stress treatments, 251 which presented sequences that could be annotated with a functional gene ontology and assigned to a GO category, despite the use of a combination of methodological approaches, including 454 and ESTs sequencing, and different genotypes analyzed (Raney et al., 2014). Indeed, these authors identified a group of 27 differentially expressed genes that can be classified as having regulatory functions, including one whose inferred polypeptide product showed high amino acid sequence homology to Naringenin, 2-oxogluturate 3-dioxygenase, which is intermediate in the biosynthesis of flavonoids in plants (Raney et al., 2014). It is interesting that the two GO terms enriched on the drought over-represented genes, the processes of “fruit ripening” and “reproduction,” has been associated to drought stress in grapevine, cereals and other plants on studies related to climate change (Castellarin et al., 2007; Barnabas et al., 2008; Chaves et al., 2010; Fleury et al., 2010) where they have been related to an escape strategy before the onset of severe drought stress. This strategy involves several biological processes, most of them directed to dehydration avoidance like minimized water loss (e.g., caused by stomatal closure, trichomes, reduced leaf area, senescence of older leaves, etc.) or maximized water uptake (e.g., by increased root growth). These changes must be equilibrated energetically and we found that, under drought condition quinoa also displayed reduction in pollen-pistil interactions, possibly related to CLE genes in response to environmental stimuli such as heat stress (Wang et al., 2016), reduction in abscission processes and reduction in carbohydrate metabolic processes suggesting different ways of metabolic plasticity under the drought condition.

The results obtained from our in silico analysis served as reference to further identify target genes for in vivo experiments, where the expression patterns of 15 genes were assessed by qPCR, which were chosen for canonical responses to drought such as ABA biosynthesis, since ABA regulation seemed to be one of the mechanisms utilized by quinoa when facing drought induced decrease of turgor of stomata guard cells (Jacobsen et al., 2009), and another group of genes directly induced in response to drought. We could determine that CqABA1, CqABA2, and CqABA4 from the ABA biosynthesis pathway were repressed over two times. One hypothesis to explain why these two genes were repressed might be that in quinoa a family of corresponding genes might be present as tetraploid species, and the designed primers did not differentiate between the members of the family. A previous study concluded that ABA apparently plays a minor role under drought conditions in quinoa, and the authors suggested that quinoa might produce other antitranspirants compounds than ABA in the xylem sap (Jacobsen et al., 2009). Also, hormonal stress signals may exist and may play an important role in quinoa, suggesting that both cytokinin and ethylene reactions should be further dissected in quinoa.

The relative abundance of studied target genes was highly variable in response to stressful conditions. There were genes that had very low expression levels, as ABA-pathway transcripts described above. In the other hand, the highest up-regulation (over 200-fold) was exhibited by CqHSP20, might indicates that drought treatment affected the thermal regulation via water shortage, since Hsp20 genes have been induced to a larger extent in tomato plants under various abiotic stresses, including heat, salt, and drought treatments (Yu et al., 2016), which may suggest that plants were inducing ABA-independent responses as well. Small HSPs (sHSPs) are involved in folding and assembling protein, keeping protein stabilization, activating protein, and degrading protein in many normal cellular processes and under stress conditions. Indeed, most of the HSPs (HSP18.1, HSP18.3, HSP20, HSP21.7, HSP25.3, HSP26.26, and HSP30) were significantly up-regulated during the heat treatment in spinach (Spinacia oleracea L.) leaves under heat stress (Yan et al., 2016). In the case of CqLEA gene that was induced 140-fold in treated plants, consistent with the fact that the ABA-responsive element was first described for a group 2 LEA gene from rice, there are genes for this group of proteins whose expression in response to stress is mediated by ABA. Moreover, there are examples of dual regulation; that is, their response to stress is mediated by more than one pathway, one of which may be ABA dependent (Welling et al., 2004; Battaglia et al., 2008). A different case was CqAP2/ERF, whose up-regulation might be pointing to a pivotal role in coordination of regulatory networks underlying abiotic stress tolerance (Golldack et al., 2014). Furthermore, functional analysis of berry transcriptomic responses to higher temperatures revealed the induction of heat shock protein (HSP) chaperones coincident with up-regulation of ERF subfamily transcription factors and increased ABA levels, suggesting their participation in the maintenance of the acclimation response (Carbonell-Bejerano et al., 2013).

Differentially abundant genes both over-represented and down-represented in drought treatment were compared among them using Blast2GO. The categories that shown over-represented genes that are commonly associated to drought were “response to stress” as expected with 85 genes and “response to abiotic stimulus” with 57 genes. Among the biological processes that were found in both subsets but over-represented in drought we highlight: “biosynthetic process,” “single-organism metabolic process,” “nitrogen compound metabolic process,” “establishment of localization,” “single-multicellular organism process” and “single-organism developmental process.” Also there were four processes only detectable in the over-represented subset: “cellular response to stimulus” (18 genes), “single organism signaling” (18 genes), “developmental process involved in reproduction” (13 genes) and “single organism reproductive process” (13 genes). In general, these results coincide with reports for Thellungiella, wheat, poplar, rice, cassava, Arabidopsis and chickpea (Wong et al., 2005; Mochida et al., 2006; Street et al., 2006; Gorantla et al., 2007; Lokko et al., 2007; Huang et al., 2008; Varshney et al., 2009). These results indicated that in quinoa most responses to drought are conserved, as for most species where this abiotic stress has been studied. The difference in tolerance levels might rely in the population of genes without identity. When looking at categories of genes that were regulated with stress, we highlighted the plasma membrane related genes, because these variations can be linked to cell turgor maintenance (Razzaghi et al., 2015) or membrane stability observed by measuring electrolyte leakage. Indeed, we found an up-representation of transporter activity genes such as ABC, ERD6-like, MATE, SWEET-like among others, in parallel with down-representation of cell wall modifying genes. These variations can be attributed to the ability to withstand stress levels to which the plant was submitted, remobilization of assimilates and might be linked to ROS (reactive oxygen species) detoxification, since they cause severe cellular damage by peroxidation and de-esterification of membrane-lipids (Golldack et al., 2014; Raney et al., 2014).

Conclusions

This work has determined that quinoa transcriptome from a tolerant genotype R49 has exhibited a 15% of genes (382 contigs) that not presented homology to the published databases from 2,456 identified differentially represented transcripts under drought conditions. The over-represented genes were higher (1,579) than down-represented genes (877) by drought treatment, and 19% of over-represented genes (306 contigs) were unknown. In Arabidopsis, approximately 40% of enzyme- and transporter-encoding genes have credible functional annotations, and this number is even lower in non-model plants. The slight up-regulation of ABA genes in response to drought stress in quinoa might indicate that ABA-independent mechanisms are committed to coordinate responses to acclimate to hydric deficit. Also functional characterization of unknown genes remains a challenge, but various databases and homologs cross-kingdom comparative genomics could be mined to provide clues (Niehaus et al., 2015). Therefore, integrative efforts are still necessary to unravel how this resilient species is able to withstand in adverse environments where others species fail to succeed. Indeed, the information represents a very useful tool for selecting drought tolerant parentals or lines with active tolerance mechanisms for breeding purposes, being useful to explore the differentially expressed gene space and valuable for loci identification in ongoing quinoa breeding efforts.

Author Contributions

AM performed drought treatments, prepared RNA samples and interpreted the results. JM performed assembling, expression analysis and interpreted the results. AM, HS, and AZ designed the experiment and provided guidance of the study. AZ and HS provided assistance in the data analysis, management data interpretation and wrote the manuscript. All authors have read and approved the final manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

HS was supported by Iniciativa Científica Milenio, MNPCB P06-065-F and CONICYT, FONDECYT/Regular No. 1120261; AZ was supported by Innova Chile (BioTecZA 06FC01IBC-71), FIC-R Atacama BIP-30432772-0, CONICYT FONDECYT/Regular No. 1140039 2014/INIA and Regional Research Center INIA-Intihuasi; AM was supported by a CONICYT Fellowship D-24110186, Universidad Andres Bello Fellowship DI-22-11//I and was a member of the Biotechnology Doctoral Program at Universidad Andres Bello. Authors want to thank technical assistance of Msc María Alejandra Montoya.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017.00216/full#supplementary-material

Figure S1. Distribution of contigs size obtained from the genome reference (hybrid assembly).

Table S1. Soil water potential for quinoa plants ecotype R49.

Table S2. qPCR primers list for selected genes in this study.

Additional File 1. FASTA file with predicted full length CDS contigs used for further analysis as our reference transcriptome.

Additional File 2. Excel spreadsheet with differentially expressed genes (DEG analysis) when compared reads from drought and control conditions. Proportions fold of change of values minus infinite and plus infinite was set to −1,000,000 and 1,000,000 respectively.

Additional File 3. Excel spreadsheet with the annotation and gene ontology of transcripts with predicted full length CDS.

Additional File 4. Excel spreadsheet with the validation of qPCR amplicons by BLAST alignment against quinoa transcriptome. Amplicons were sequenced using Sanger technology.

Additional File 5. Excel spreadsheet with the structure validation of selected quinoa genes by BLAST alignment against the genome predicted transcripts (Yasui et al., 2016).

Availability of Data and Materials

The nucleotide sequences of raw reads from this study were submitted to NCBI's Sequence Read Archive through the BioProject ID: PRJNA305752. Contigs with predicted full length CDS are available in Additional File 3.

Abbreviations

RWC, Relative Water Content; EL, Electrolyte Leakage; qPCR, Quantitative Real-Time PCR; ABA, Abscisic acid; BP, Biological process; FC, Fold-change; GO, Gene ontology.

References

Ali-Benali, M., Alary, R., Joudrier, P., and Gautier, M. (2005). Comparative expression of five Lea Genes during wheat seed development and in response to abiotic stresses by real-time quantitative RT-PCR. Biochim. Biophys. Acta 1730, 56–65. doi: 10.1016/j.bbaexp.2005.05.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Alvarez-Flores, R., Winkel, T., Nguyen-Thi-Truc, A., and Joffre, R. (2014). Root foraging capacity depends on root system architecture and ontogeny in seedlings of three Andean Chenopodium species. Plant Soil 380, 415–428. doi: 10.1007/s11104-014-2105-x

CrossRef Full Text | Google Scholar

Avramova, V., AbdElgawad, H., Zhang, Z., Fotschki, B., Casadevall, R., Vergauwen, L., et al. (2015). Drought induces distinct growth response, protection and recovery mechanisms in the maize leaf growth zone. Plant Physiol. 169, 1382–1396. doi: 10.1104/pp.15.00276

PubMed Abstract | CrossRef Full Text | Google Scholar

Barnabas, B., Jager, K., and Feher, A. (2008). The effect of drought and heat stress on reproductive processes in cereals. Plant Cell Environ. 31, 11–38. doi: 10.1111/j.1365-3040.2007.01727.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrs, H. D., and Weatherley, P. E. (1962). A re-examination of the relative turgidity technique for estimating water deficits in leaves. Aust. J. Biol. Sci. 15, 413–428. doi: 10.1071/BI9620413

CrossRef Full Text | Google Scholar

Battaglia, M., Olvera-Carrillo, Y., Garciarrubio, A., Campos, F., and Covarrubias, A. A. (2008). The enigmatic LEA proteins and other hydrophilins. Plant Physiol. 148, 6–24. doi: 10.1104/pp.108.120725

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazile, D., Bertero, H. D., and Nieto, C. (eds.) (2015). State of the Art Report on Quinoa Around the World in 2013. Rome: FAO, CIRAD, 603. Available online at: http://www.fao.org/3/a-i4042e.pdf

Google Scholar

Bertero, H. D., de la Vega, A. J., Correa, G., Jacobsen, S. E., and Mujica, A. (2004). Genotype and genotype-by-environment interaction effects for grain yield and grain size of quinoa (Chenopodium quinoa Willd.) as revealed by pattern analysis of international multi-environment trials. Field Crops Res. 89, 299–318. doi: 10.1016/j.fcr.2004.02.006

CrossRef Full Text | Google Scholar

Bhargava, A., Shukla, S., and Ohri, D. (2006). Chenopodium quinoa. An Indian perspective. Indust. Crops Prod. 23, 73–87. doi: 10.1016/j.indcrop.2005.04.002

CrossRef Full Text | Google Scholar

Bolstad, B. M., Irizarry, R. A., Åstrand, M., and Speed, T. P. (2003). A comparison of normalization methods for high-density oligonucleotide array data based on variance and bias. Bioinformatics 19, 185–193. doi: 10.1093/bioinformatics/19.2.185

PubMed Abstract | CrossRef Full Text | Google Scholar

Bosque Sanchez, H., Lemeur, R., Van Damme, P., and Jacobsen, S. E. (2003). Ecophysiological analysis of drought and salinity stress of quinoa (Chenopodium quinoa Willd.). Food Rev. Int. 19, 111–119. doi: 10.1081/FRI-120018874

CrossRef Full Text | Google Scholar

Candat, A., Paszkiewicz, G., Neveu, M., Gautier, R., Logan, D. C., Avelange-Macherel, M. H., et al. (2014). The ubiquitous distribution of late embryogenesis abundant proteins across cell compartments in Arabidopsis offers tailored protection against abiotic stress. Plant Cell 26, 3148–3166. doi: 10.1105/tpc.114.127316

PubMed Abstract | CrossRef Full Text | Google Scholar

Carbonell-Bejerano, P., Santa María, E., Torres-Pérez, R., Royo, C., Lijavetzky, D., Bravo, G., et al. (2013). Thermotolerance Responses in Ripening Berries of Vitis vinifera L. cv Muscat Hamburg. Plant Cell Physiol. 54, 1200–1216. doi: 10.1093/pcp/pct071

PubMed Abstract | CrossRef Full Text | Google Scholar

Castellarin, S. D., Pfeiffer, A., Sivilotti, P., Degan, M., Peterlunger, E., and Gaspero, D. I. (2007). Transcriptional regulation of anthocyanin biosynthesis in ripening fruits of grapevine under seasonal water deficit. Plant Cell Environ. 30, 1381–1399. doi: 10.1111/j.1365-3040.2007.01716.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaves, M. M., Zarrouk, O., Francisco, R., Costa, J. M., Santos, T., Regalado, A. P., et al. (2010). Grapevine under deficit irrigation: hints from physiological and molecular data. Ann. Bot. 105, 661–676. doi: 10.1093/aob/mcq030

PubMed Abstract | CrossRef Full Text | Google Scholar

Conesa, A., and Gotz, S. (2008). Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int. J. Plant Genomics 2008:619832. doi: 10.1155/2008/619832

PubMed Abstract | CrossRef Full Text | Google Scholar

Crawford, J., Guelbeogo, W., Sanou, A., Traoré, A., Vernick, K., Sagnon, N., et al. (2010). De novo transcriptome sequencing in Anopheles funestus using illumina RNA-Seq technology. PLoS ONE 5:e14202. doi: 10.1371/journal.pone.0014202

PubMed Abstract | CrossRef Full Text | Google Scholar

Food and Agriculture Organization of the United Nations (2011). Climate Change, Water and Food Security, Vol. 36. FAO Water Reports, Rome.

Fleury, D., Jefferies, S., Kuchel, H., and Langridge, P. (2010). Genetic and genomic tools to improve drought tolerance in wheat. J. Exp. Bot. 61, 3211–3222. doi: 10.1093/jxb/erq152

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuentes, F., and Bhargava, A. (2011). Morphological analysis of Quinoa germplasm grown under lowland desert conditions. J. Agron. Crop Sci. 197, 124–134. doi: 10.1111/j.1439-037X.2010.00445.x

CrossRef Full Text | Google Scholar

Fuentes, F. F., Bazile, D., Bhargava, A., and Martínez, E. A. (2012). Implications of farmers' seed exchanges for on-farm conservation of quinoa, as revealed by its genetic diversity in Chile. J. Agric. Sci. 150, 702–716. doi: 10.1017/S0021859612000056

CrossRef Full Text | Google Scholar

Fuentes, F. F., Martínez, E. A., Hinrichsen, P. V., Jellen, E. N., and Maughan, P. J. (2008). Assessment of genetic diversity patterns in Chilean quinoa (Chenopodium quinoa Willd.) germplasm using multiplex fluorescent microsatellite markers. Conserv. Genet. 10, 369–377. doi: 10.1007/s10592-008-9604-3

CrossRef Full Text | Google Scholar

Galwey, N. W. (1989). Quinoa. Biologist 36:5.

Google Scholar

Golldack, D., Li, C., Mohan, H., and Probst, N. (2014). Tolerance to drought and salt stress in plants: unraveling the signaling networks. Front. Plant Sci. 5:151. doi: 10.3389/fpls.2014.00151

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorantla, M., Babu, P. R., Lachagari, V. B., Reddy, A. M., Wusirika, R., Bennetzen, J. L., et al. (2007). Identification of stress-responsive genes in an indica rice (Oryza sativa L.) using ESTs generated from drought-stressed seedlings. J. Exp. Bot. 58, 253–265. doi: 10.1093/jxb/erl213

PubMed Abstract | CrossRef Full Text | Google Scholar

He, M., Wang, Y., Hua, W., Zhang, Y., and Wang, Z. (2012). De novo sequencing of Hypericum perforatum transcriptome to identify potential genes involved in the biosynthesis of active metabolites. PLoS ONE 7:e42081. doi: 10.1371/journal.pone.0042081

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, D., Wu, W., Abrams, S. R., and Cutler, A. J. (2008). The relationship of drought-related gene expression in Arabidopsis thaliana to hormonal and environmental factors. J. Exp. Bot. 59, 2991–3007. doi: 10.1093/jxb/ern155

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J., Lu, X., Yan, H., Chen, S., Zhang, W., Huang, R., et al. (2012). Transcriptome characterization and sequencing-based identification of salt-responsive genes in Millettia pinnata, a semi-mangrove plant. DNA Res. 19, 195–207. doi: 10.1093/dnares/dss004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hunter, S., Apweiler, R., Attwood, T. K., Bairoch, A., Bateman, A., Binns, D., et al. (2009). InterPro: the integrative protein signature database. Nucleic Acids Res. 37, D211–D215. doi: 10.1093/nar/gkn785

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobsen, S. E. (2003). The worldwide potential for quinoa (Chenopodium quinoa Willd.). Food Rev. Int. 19, 167–177. doi: 10.1081/FRI-120018883

CrossRef Full Text | Google Scholar

Jacobsen, S. E., Liu, F., and Jensen, C. R. (2009). Does root-sourced ABA play a role for regulation of stomata under drought in quinoa (Chenopodium quinoa Willd.). Sci. Hortic. 122, 281–287. doi: 10.1016/j.scienta.2009.05.019

CrossRef Full Text | Google Scholar

Jacobsen, S. E., Mujica, A., and Jensen, C. R. (2003). The resistance of quinoa (Chenopodium quinoa Willd.) to adverse abiotic factors. Food Rev. Int. 19, 99–109. doi: 10.1081/FRI-120018872

CrossRef Full Text | Google Scholar

Jarvis, D. E., Ho, Y. S., Lightfoot, D. J., Schmöckel, S. M., Li, B., Borm, T. J. A., et al. (2017). The genome of Chenopodium quinoa. Nature 542, 307–312. doi: 10.1038/nature21370

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, C. R., Jacobsen, S. E., Andersen, M. N., Nunez, N., Andersen, S. D., Rasmussen, L., et al. (2000). Leaf gas exchange and water relation characteristics of field quinoa (Chenopodium quinoa Willd.) during soil drying. Eur. J. Agron. 13, 11–25. doi: 10.1016/S1161-0301(00)00055-1

CrossRef Full Text | Google Scholar

Kal, A. J., van Zonneveld, A. J., Benes, V., van den Berg, M., Koerkamp, M. G., Albermann, K., et al. (1999). Dynamics of gene expression revealed by comparison of serial analysis of gene expression transcript profiles from yeast grown on two different carbon sources. Mol. Biol. Cell 10, 1859–1872. doi: 10.1091/mbc.10.6.1859

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaye, C., Neven, L., Hofig, A., Li, Q. B., Haskell, D., and Guy, C. (1998). Characterization of a gene for spinach CAP160 and expression of two spinach cold-acclimation proteins in tobacco. Plant Physiol. 116, 1367–1377. doi: 10.1104/pp.116.4.1367

PubMed Abstract | CrossRef Full Text | Google Scholar

Lokko, Y., Anderson, J. V., Rudd, S., Raji, A., Horvath, D., Mikel, M. A., et al. (2007). Characterization of an 18,166 EST dataset for cassava (Manihot esculenta Crantz) enriched for drought-responsive genes. Plant Cell Rep. 26, 1605–1618. doi: 10.1007/s00299-007-0378-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Martínez, E. A., Fuentes, F., and Bazile, D. (2015). “History of quinoa: its origin, domestication, diversification and cultivation with particular reference to the Chilean context,” in Quinoa: Improvement and Sustainable Production, World Agriculture Series, eds K. Murphy and J. Matanguihan (Hoboken, NJ: Wiley-Blackwell,) 19–24.

Google Scholar

Martínez, E. A., Veas, E., Jorquera, C., San Martin, R., and Jara, P. (2009). Re-introduction of quinoa into arid Chile: cultivation of two lowland races under extremely low irrigation. J. Agron. Crop Sci. 195, 1–10. doi: 10.1111/j.1439-037X.2008.00332.x

CrossRef Full Text | Google Scholar

Mastebroek, H., Limburg, H., Gilles, T., and Marvin, H. (2000). Occurrence of sapogenins in leaves and seeds of quinoa (Chenopodium quinoa Willd). J. Sci. Food Agric. 80, 152–156. doi: 10.1002/(SICI)1097-0010(20000101)80:1<152::AID-JSFA503>3.0.CO;2-P

CrossRef Full Text | Google Scholar

Maughan, P. J., Bonifacio, A., Jellen, E., Stevens, M., Coleman, C., Ricks, M., et al. (2004). A genetic linkage map of quinoa (Chenopodium quinoa) based on AFLP, RAPD and SSR markers. Theor. Appl. Genet. 109, 1188–1189. doi: 10.1007/s00122-004-1730-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Meisel, L., Fonseca, B., Gonzalez, S., Baeza-Yates, R., Cambiazo, V., Campos, R., et al. (2005). A rapid and efficient method for purifying high quality total RNA from peaches (Prunus persica) for functional genomics analyses. Biol. Res. 38, 83–88. doi: 10.4067/s0716-97602005000100010

PubMed Abstract | CrossRef Full Text | Google Scholar

Merewitz, E. B., Gianfagna, T., and Huang, B. (2011). Protein accumulation in leaves and roots associated with improved drought tolerance in creeping bentgrass expressing an ipt gene for cytokinin synthesis. J. Exp. Bot. 62, 5311–5333. doi: 10.1093/jxb/err166

PubMed Abstract | CrossRef Full Text | Google Scholar

Mestanza, C., Riegel, R., Silva, H., and Vásquez, S. (2015). Characterization of the acetohydroxyacid synthase multigene family in the tetraploide plant Chenopodium quinoa. Electron. J. Biotechnol. 18, 393–398. doi: 10.1016/j.ejbt.2015.07.003

CrossRef Full Text | Google Scholar

Mizrachi, E., Hefer, C., Ranik, M., Joubert, F., and Myburg, A. (2010). De novo assembled expressed gene catalog of a fast-growing Eucalyptus tree produced by Illumina mRNA-Seq. BMC Genomics 11:681. doi: 10.1186/1471-2164-11-681

PubMed Abstract | CrossRef Full Text | Google Scholar

Mochida, K., Kawaura, K., Shimosaka, E., Kawakami, N., Shin, I. T., Kohara, Y., et al. (2006). Tissue expression map of a large number of expressed sequence tags and its application to in silico screening of stress response genes in common wheat. Mol. Genet. Genomics 276, 304–312. doi: 10.1007/s00438-006-0120-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Mortazavi, A., Williams, B. A., McCue, K., Schaeffer, L., and Wold, B. (2008). Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. doi: 10.1038/nmeth.1226

PubMed Abstract | CrossRef Full Text | Google Scholar

Mujica, A., and Jacobsen, S. (2006). La quinua (Chenopodium quinoa Willd.) y Sus Parientes Silvestres. La Paz: Botánica Económica de los Andes Centrales, Universidad Mayor de San Andrés.

Google Scholar

Nakano, T., Suzuki, K., Fujimura, T., and Shinshi, H. (2006). Genome-wide analysis of the ERF gene family in Arabidopsis and rice. Plant Physiol. 140, 411–432. doi: 10.1104/pp.105.073783

PubMed Abstract | CrossRef Full Text | Google Scholar

Niehaus, T. D., Thamm, A. M. K., de Crécy-Lagard, V., and Hanson, A. D. (2015). Proteins of unknown biochemical function: a persistent problem and a roadmap to help overcome it. Plant Physiol. 169, 1436–1442. doi: 10.1104/pp.15.00959

PubMed Abstract | CrossRef Full Text | Google Scholar

Orsini, F., Accorsi, M., Gianquinto, G., Dinelli, G., Antognoni, F., Ruiz-Carrasco, K., et al. (2011). Beyond the ionic and osmotic response to salinity in Chenopodium quinoa: functional elements of successful halophytism. Funct. Plant Biol. 38, 818–831. doi: 10.1071/FP11088

CrossRef Full Text | Google Scholar

Pinhero, R., Rao, M., Paliyath, G., Murr, D., and Fletcher, R. (1997). Changes in activities of antioxidant enzymes and their relationship to genetic and paclobutrazol-induced chilling tolerance of maize seedlings. Plant Physiol. 114, 695–704. doi: 10.1104/pp.114.2.695

PubMed Abstract | CrossRef Full Text | Google Scholar

Pruitt, K. D., Tatusova, T., and Maglott, D. R. (2007). NCBI reference sequences (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res. 35, D61–D65. doi: 10.1093/nar/gkl842

PubMed Abstract | CrossRef Full Text | Google Scholar

Raney, J., Reynolds, D. J., Elzinga, D. B., Page, J., Udall, J. A., Jellen, E. N., et al. (2014). Transcriptome analysis of drought induced stress in Chenopodium quinoa. Am. J. Plant Sci. 5, 338–357. doi: 10.4236/ajps.2014.53047

CrossRef Full Text | Google Scholar

Razzaghi, F., Jacobsen, S.-E., Jensen, C. R., and Andersen, M. N. (2015). Ionic and photosynthetic homeostasis in quinoa challenged by salinity and drought – mechanisms of tolerance. Funct. Plant Biol. 42, 136–148. doi: 10.1071/FP14132

CrossRef Full Text | Google Scholar

Repo-Carrasco, R., Espinoza, C., and Jacobsen, S. E. (2003). Nutritional value and use of the Andean crops Quinoa (Chenopodium quinoa) and Ka-iwa (Chenopodium pallidicaule). Food Rev. Int. 19, 179–189. doi: 10.1081/FRI-120018884

CrossRef Full Text | Google Scholar

Rice, P., Longden, I., and Bleasby, A. (2000). EMBOSS: the European Molecular Biology Open Software Suite. Trends Genet. 16, 276–277. doi: 10.1016/S0168-9525(00)02024-2

PubMed Abstract | CrossRef Full Text

Riggins, C. W., Peng, Y., Stewart, C. N. Jr., and Tranel, P. J. (2010). Characterization of de novo transcriptome for waterhemp (Amaranthus tuberculatus) using GS-FLX 454 pyrosequencing and its application for studies of herbicide target-site genes. Pest Manag. Sci. 66, 1042–1052. doi: 10.1002/ps.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Risi, J. C., and Galwey, N. W. (1989a). The pattern of genetic diversity in the Andean grain crop quinoa (Chenopodium quinoa Willd.). I. Associations between characteristics. Euphytica 41, 147–162. doi: 10.1007/BF00022424

CrossRef Full Text | Google Scholar

Risi, J. C., and Galwey, N. W. (1989b). The pattern of genetic diversity in the Andean grain crop quinoa (Chenopodium quinoa Willd.). II. Multivariate methods. Euphytica 41, 135–145. doi: 10.1007/BF00022423

CrossRef Full Text | Google Scholar

Risi, J., and Galwey, N. W. (1984). The Chenopodium grains of the Andes: Inca crops for modern agriculture. Adv. Appl. Biol. 10, 145–216.

Google Scholar

Ruiz, K. B., Aloisi, I., Del Duca, S., Canelo, V., Torrigiani, P., Silva, H., et al. (2016). Salares versus coastal ecotypes of quinoa: salinity responses in Chilean landraces from contrasting habitats. Plant Physiol. Biochem. 101, 1–13. doi: 10.1016/j.plaphy.2016.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruiz, K. B., Biondi, S., Oses, R., Acuña-Rodríguez, I. S., Antognoni, F., Martinez-Mosqueira, E. A., et al. (2014). Quinoa biodiversity and sustainability for food security under climate change. A review. Agron. Sustain. Dev. 34, 349–359. doi: 10.1007/s13593-013-0195-0

CrossRef Full Text | Google Scholar

Ruiz-Carrasco, K., Antognoni, F., Coulibaly, A. K., Lizardi, S., Covarrubias, A., Martínez, E. A., et al. (2011). Variation in salinity tolerance of four lowland genotypes of quinoa (Chenopodium quinoa Willd.) as assessed by growth, physiological traits, and sodium transporter gene expression. Plant Physiol. Biochem. 49, 1333–1341. doi: 10.1016/j.plaphy.2011.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Schulz, M. H., Zerbino, D. R., Vingron, M., and Birney, E. (2012). Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics 28, 1086–1092. doi: 10.1093/bioinformatics/bts094

PubMed Abstract | CrossRef Full Text | Google Scholar

Seki, M., Narusaka, M., Ishida, J., Nanjo, T., Fujita, M., Oono, Y., et al. (2002). Monitoring the expression profiles of 7000 Arabidopsis genes under drought, cold and high-salinity stresses using a full-length cDNA microarray. Plant J. 31, 279–292. doi: 10.1046/j.1365-313X.2002.01359.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Seo, M., and Koshiba, T. (2011). Transport of ABA from the site of biosynthesis to the site of action. J. Plant Res. 124, 501–507. doi: 10.1007/s10265-011-0411-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, S., and Verslues, P. E. (2010). Mechanisms independent of abscisic acid (ABA) or proline feedback have a predominant role in transcriptional regulation of proline metabolism during low water potential and stress recovery. Plant Cell Environ. 33, 1838–1851. doi: 10.1111/j.1365-3040.2010.02188.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shinozaki, K., and Yamaguchi-Shinozaki, K. (2007). Gene networks involved in drought stress response and tolerance. J. Exp. Bot. 58, 221–227. doi: 10.1093/jxb/erl164

PubMed Abstract | CrossRef Full Text | Google Scholar

Shinozaki, K., Yamaguchi-Shinozaki, K., and Seki, M. (2003). Regulatory network of gene expression in the drought and cold stress responses. Curr. Opin. Plant Biol. 6, 410–417. doi: 10.1016/S1369-5266(03)00092-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Siener, R., Honow, R., Seidler, A., Voss, S., and Hesse, A. (2006). Oxalate contents of species of the Polygonaceae, Amaranthaceae and Chenopodiaceae families. Food Chem. 98, 220–224. doi: 10.1016/j.foodchem.2005.05.059

CrossRef Full Text | Google Scholar

Simmons, N. W. (1971). The breeding system of Chenopodium quinoa. I. Male Sterility. Heredity 27, 73–82. doi: 10.1038/hdy.1971.72

CrossRef Full Text | Google Scholar

Simnadis, T. G., Tapsell, L. C., and Beck, E. J. (2015). Physiological effects associated with Quinoa consumption and implications for research involving humans: a review. Plant Foods Hum. Nutr. 70, 238–249. doi: 10.1007/s11130-015-0506-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Stevens, M. R., Coleman, C. E., Parkinson, S. E., Maughan, P. J., Zhang, H. B., Balzotti, M. R., et al. (2006). Construction of a quinoa (Chenopodium quinoa Willd.) BAC library and its use in identifying genes encoding seed storage proteins. Theor. Appl. Genet. 112, 1593–1600. doi: 10.1007/s00122-006-0266-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Stikic, R., Glamoclija, D., Demin, M., Vucelic-Radovic, B., Jovanovic, Z., Milojkovic-Opsenica, D., et al. (2012). Agronomical and nutritional evaluation of quinoa seeds (Chenopodium quinoa Willd.) as an ingredient in bread formulations. J. Cereal Sci. 55, 132–138. doi: 10.1016/j.jcs.2011.10.010

CrossRef Full Text | Google Scholar

Street, N. R., Skogström, O., Sjödin, A., Tuckere, J., Rodriguez-Costa, M., Nilsson, P., et al. (2006). The genetics and genomics of the drought response in Populus. Plant J. 48, 321–341. doi: 10.1111/j.1365-313X.2006.02864.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Y., Liu, F., Bendevis, M., Shabala, S., and Jacobsen, S. E. (2014). Sensitivity of two Quinoa (Chenopodium quinoa Willd.) varieties to progressive drought stress. J. Agron. Crop Sci. 200, 12–23. doi: 10.1111/jac.12042

CrossRef Full Text | Google Scholar

Swindell, W. R., Huebner, M., and Weber, A. P. (2007). Transcriptional profiling of Arabidopsis heat shock proteins and transcription factors reveals extensive overlap between heat and non-heat stress response pathways. BMC Genomics 8:125. doi: 10.1186/1471-2164-8-125

PubMed Abstract | CrossRef Full Text | Google Scholar

Tittarelli, A., Santiago, M., Morales, A., Meisel, L., and Silva, H. (2009). Isolation and functional characterization of cold-regulated promoters, by digitally identifying peach fruit cold-induced genes from a large EST dataset. BMC Plant Biol. 9:121. doi: 10.1186/1471-2229-9-121

PubMed Abstract | CrossRef Full Text | Google Scholar

Umezawa, T., Sugiyama, N., Mizoguchi, M., Hayashi, S., Myouga, F., Yamaguchi-Shinozaki, K., et al. (2009). Type 2C protein phosphatases directly regulate abscisic acid-activated protein kinases in Arabidopsis. Proc. Natl. Acad. Sci. U.S.A. 106, 17588–17593. doi: 10.1073/pnas.0907095106

PubMed Abstract | CrossRef Full Text | Google Scholar

Varshney, R. K., Hiremath, P. J., Lekha, P., Kashiwagi, J., Balaji, J., Deokar, A. A., et al. (2009). A comprehensive resource of drought- and salinity responsive ESTs for gene discovery and marker development in chickpea (Cicer arietinum L.). BMC Genomics 10:523. doi: 10.1186/1471-2164-10-523

PubMed Abstract | CrossRef Full Text | Google Scholar

Vega-Gálvez, A., Miranda, M., Vergara, J., Uribe, E., Puente, L., and Martinez, E. A. (2010). Nutrition facts and functional potential of quinoa (Chenopodium quinoa Willd.), an ancient Andean grain: a review. J. Sci. Food Agric. 90, 2541–2547. doi: 10.1002/jsfa.4158

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, G., Zhang, G., and Wu, M. (2016). CLE Peptide signaling and crosstalk with phytohormones and environmental stimuli. Front. Plant Sci. 6:1211. doi: 10.3389/fpls.2015.01211

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Vinocur, B., Shoseyov, O., and Altman, A. (2004). Role of plant heatshock proteins and molecular chaperones in the abiotic stress response. Trends Plant Sci. 9, 244–252. doi: 10.1016/j.tplants.2004.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Ward, S. M. (2000). Allotetraploid segregation for single-gene morphological characters in quinoa (Chenopodium quinoa Willd.). Euphytica 116, 11–16. doi: 10.1023/A:1004070517808

CrossRef Full Text | Google Scholar

Welling, A., Rinne, P., Vihera-Aarnio, A., Kontunen-Soppela, S., Heino, P., and Palva, E. T. (2004). Photoperiod and temperature differentially regulate the expression of two dehydrin genes during overwintering of birch (Betula pubescens Ehrh.). J. Exp. Bot. 55, 507–516. doi: 10.1093/jxb/erh045

PubMed Abstract | CrossRef Full Text | Google Scholar

Winkel, T., Méthy, M., and Thénot, F. (2002). Radiation-use efficiency, chlorophyll fluorescence and reflectance indices, associated with ontogenic changes in water-limited Chenopodium quinoa leaves. Photosynthetica 40, 227–232. doi: 10.1023/A:1021345724248

CrossRef Full Text | Google Scholar

Wong, C. E., Li, Y., Whitty, B. R., Díaz-Camino, C., Akhter, S. R., Brandle, J. E., et al. (2005). Expressed sequence tags from the Yukon ecotype of Thellungiella reveal that gene expression in response to cold, drought and salinity shows little overlap. Plant Mol. Biol. 58, 561–574. doi: 10.1007/s11103-005-6163-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Woo, N. S., Badger, M. R., and Pogson, B. J. (2008). A rapid, non-invasive procedure for quantitative assessment of drought survival using chlorophyll fluorescence. Plant Methods 4:27. doi: 10.1186/1746-4811-4-27

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, J., Yu, L., Xuan, J., Lu, Y., Lu, S., and Zhu, W. (2016). De novo transcriptome sequencing and gene expression profiling of spinach (Spinacia oleracea L.) leaves under heat stress. Sci. Rep. 6, 1–10. doi: 10.1038/srep19473

PubMed Abstract | CrossRef Full Text | Google Scholar

Yasui, Y., Hirakawa, H., Oikawa, T., Toyoshima, M., Matsuzaki, C., Ueno, M., et al. (2016). Draft genome sequence of an inbred line of Chenopodium quinoa, an allotetraploid crop with great environmental adaptability and outstanding nutritional properties. DNA Res. 23, 535–546. doi: 10.1093/dnares/dsw037

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, J., Cheng, Y., Feng, K., Ruan, M., Ye, Q., Wang, R., et al. (2016). Genome-wide identification and expression profiling of tomato Hsp20 gene family in response to biotic and abiotic stresses. Front. Plant Sci. 7:1215. doi: 10.3389/fpls.2016.01215

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Pei, X., Zhang, C., Lu, Z., Wang, Z., Jia, S., et al. (2012). De novo foliar transcriptome of Chenopodium amaranticolor and analysis of its gene expression during virus-induced hypersensitive response. PLoS ONE 7:e45953. doi: 10.1371/journal.pone.0045953

PubMed Abstract | CrossRef Full Text | Google Scholar

Zurita-Silva, A., Fuentes, F., Zamora, P., Jacobsen, S. E., and Schwember, A. R. (2014). Breeding quinoa (Chenopodium quinoa Willd.): potential and perspectives. Mol. Breed. 34, 13–30. doi: 10.1007/s11032-014-0023-5

CrossRef Full Text | Google Scholar