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).

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).

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.

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.