Transcriptome analysis revealed gene regulatory network involved in PEG-induced drought stress in Tartary buckwheat (Fagopyrum Tararicum)

Tartary buckwheat is a nutritious pseudo-cereal crop that is resistant to abiotic stresses, such as drought. However, the buckwheat’s mechanisms for responding to drought stress remains unknown. We investigated the changes in physiology and gene expression under drought stress, which was simulated by treatment with polyethylene glycol (PEG). Five physiological indexes, namely MDA content, H2O2 content, CAT activity, SOD activity, and POD activity, were measured over time after 20% PEG treatment. All indexes showed dramatic changes in response to drought stress. A total of 1,190 differentially expressed genes (DEGs) were identified using RNA-seq and the most predominant were related to a number of stress-response genes and late embryogenesis abundant (LEA) proteins. DEGs were gathered into six clusters and were found to be involved in the ABA biosynthesis and signal pathway based on hierarchical clustering and GO and KEGG pathway enrichment. Transcription factors, such as NAC and bZIP, also took part in the response to drought stress. We determined an ABA-dependent and ABA-independent pathway in the regulation of drought stress in Tartary buckwheat. To the best of our knowledge, this is the first transcriptome analysis of drought stress in Tartary buckwheat, and our results provide a comprehensive gene regulatory network of this crop in response to drought stress.


INTRODUCTION
Tartary buckwheat (Fagopyrum Tararicum), also called bitter buckwheat, is a pseudocereal crop belonging to the genus Fagopyrum Mill, Polygonaceae (Ohnishi, 1998). Tartary buckwheat has become popular for its rich nutritional composition that includes a high content of flavonoids, resistant starches, crude fibers, proteins, and vitamins, which are all shown to have health benefits (Frias et al., 2011;Qin et al., 2010;Zhu, 2016). Tartary buckwheat is also highly adaptable to adverse soil and climatic conditions and shows a very strong tolerance or resistance to the adverse environment and abiotic stresses, including drought, low temperatures, and acid soils (Zhang et al., 2017).
Drought is a major meteorological disaster in agriculture of China and leads to the reduction of crop yield and quality (Xu et al., 2015). A series of physiological responses occur in Tartary buckwheat under water stress (Chen et al., 2008;Xiang et al., 2013), including a decrease of chlorophyll a, chlorophyll b, and total chlorophyll content and a significant increase in proline accumulation (Xiang et al., 2013). Photosynthesis, transpiration, stomatal conductance, and yield decrease but intercellular CO 2 concentration increases (Xiang et al., 2013). Peroxidase (POD) activity, superoxide dismutase (SOD) activity, catalase (CAT) activity, hydrogen peroxide (H 2 O 2 ) content, and soluble protein content, and proline content increase significantly and relative water content decreases rapidly in Tartary buckwheat under drought stress (Chen et al., 2008;Xiang et al., 2013).
Previous studies on Tartary buckwheat focused on changes in gene regulation under drought and identified some drought-inducible TFs, including three MYB family genes (FtMYB9,FtMYB10,and FtMYB13), eight NAC family genes (from FtNAC2 to FtNAC9), two bZIP family genes (FtbZIP5 and FtbZIP83), and a HLH family gene (FtbHLH3) (Deng et al., 2019;Gao et al., 2016;Gao et al., 2017;Huang et al., 2018;Li et al., 2019;Li et al., 2020;Sun et al., 2019;Yao et al., 2017). The transcriptome analysis of Tartary buckwheat under drought stress has not been reported despite the reports of other transcriptome analyses related to seed development (Huang et al., 2017;Liu et al., 2018). Many drought-inducible genes are unknown and the regulation mechanism of Tartary buckwheat under drought stress is still under investigation.
We examined the dynamic changes of five physiological indexes under drought stress, which was simulated by 20% polyethylene glycol (PEG-6000) treatment. We then performed a comprehensive transcriptome analysis using high throughput RNA-seq. The identification of differentially expressed genes (DEGs), the hierarchical cluster, GO enrichment, and KEGG enrichment was then analyzed. Genes in ABA biosynthesis and signal transduction were significantly enriched; thereafter ABA content after drought treatment was measured. The expression patterns of 18 important DEGs were verified using qRT-PCR. We proposed an ABA-dependent and ABA-independent pathway in the regulation of drought stress for Tartary buckwheat. This is the first known transcriptome analysis under drought stress for Tartary buckwheat, and provides a comprehensive gene regulatory network of Tartary buckwheat in response to drought stress.

Plant materials and growth conditions
A widely cultivated Tartary buckwheat variety, Jinqiao No. 2, was selected for its high stability and adaptability (Li et al., 2011). Seeds were surface sterilized with a 10% H 2 O 2 solution, rinsed three times with double distilled water (ddH 2 O), and placed on wet filter papers for two days to accelerate germination. The uniformly sprouted seeds were moved to rolls of papers soaked with 0%, 10%, 20%, and 30% PEG-6000 solutions, respectively, to evaluate the genotypes after drought stress. After 48 h, the root length of the seedlings was compared across the four treatments. Three biological replicates were performed, and at least 7 seedlings were included in each biological replicate.
The uniformly sprouted seeds were moved to rolls of papers soaked with double distilled water when seedlings grew two true leaves, which occurred around day seven, in order to measure the physiological indexes and transcriptome sequencing. The seedlings with consistent growth were treated with 20% PEG based on the literature (Abdel-Ghany et al., 2020). Samples were taken at 0 h, 1 h, 3 h, and 6 h, frozen in liquid nitrogen, and stored at −80 • C. Three biological replicates were performed, and 40 seedlings were included in each biological replicate.

Physiological indexes measurement
A total of 5 physiological indexes, namely malondialdehyde (MDA) content, H 2 O 2 content, CAT activity, SOD activity, and POD activity were measured, using determination kits (Sino Best Bio-Technology Co. Ltd, China) following the manufacturer's instructions. Three biological replicates were included for each treatment, and 2 technical replicates were included for each biological replicate.

Total RNA isolation and transcriptome sequencing
Total RNA was isolated using the plant RNA purification kit (TianGen Biotech Co. LTD, China), according to the manufacturer's instructions. The purified RNA samples were treated with Dnase I for 20 min to digest the genomic DNA. The quality and quantity of RNA were determined using the NanoDrop 2000 micro spectrophotometer, the Agilent 2100 Bioanalyzer, and the Agilent RNA 6000 Nano Kit. A total of 12 RNA samples were used to construct the library, and high-throughput sequencing was performed using the Illumina 4000 System (Illumina Inc., USA) with a read length of 150 bp and paired-end method.

Measurement of lycopene, zeaxanthin, and ABA content
High performance liquid chromatography was used to measure the content of lycopene, zeaxanthin, and ABA. For the measurement of lycopene and zeaxanthin, approximately 0.300 g well-ground samples were precisely weighed and placed in a 5 mL brown volumetric flask. Then 5 mL 0.1% bht-ethanol solution was added and oscillated for 5 min. The mixture was oscillated at 200 r/min for 4 h at room temperature in the dark. A 0.1% bht-ethanol solution was added to keep the volume of mixture at 10 mL. The mixture was subsequently centrifuged at 4,000 r/min for 10 min and 1 mL of the supernatant was filtered using a 0.22 µm Millipore filter, and the solution was collected in a 1.5 mL brown sample bottle. Lycopene and zeaxanthin were tested using Agilent HPLC-1100 with a DAD detector and Thermopylae C18 chromatographic column. The test conditions were as follows: column temperature, 25 • C; injection volume, 20 ul; flow rate, 1.0 ml/min; mobile phase, acetonitrile: methanol = 65:35 (V:V) isoelution.
Samples were well-ground and weighed to approximately 0.300 g to measure for ABA. Three mL were precooled and 80% methanol was added and mixed by oscillation. The mixture was sealed and stored overnight at 4 • C, then centrifuged at 5,000 r/min at 4 • C for 10 min. The supernatant was taken, and the residue was extracted ultrasonically with precooled 80% methanol at 4 • C for 30 min. This procedure was repeated twice and the supernatants were combined. The combined supernatant was blown to the aqueous phase with nitrogen at 4 • C. Three mL petroleum ether was added three times for decolorization, the aqueous phase was extracted with ethyl acetate three times then combined with the ethyl acetate phase and blown dry at 4 • C with nitrogen. The acetic acid solution (pH = 3.5) was added and purified in a SEP-PakC18 column. The eluent was eluted with methanol at room temperature and reduced to dry. ABA was tested by Agilent HPLC-1100 using the VWD detector and Agilent C18 chromatographic column (250*4.6 mm; 5 µL). The testing conditions were as follows: column temperature, 25 • C; injection volume, 10 ul; wavelength: 254 nm; flow rate, 1.0 ml/min; mobile phase, methanol: aqueous acetic acid solution (pH = 3.6) isoelution.

Quantitative RT-PCR (qRT-PCR) analysis
The transcriptome results were verified by qRT-PCR. A total of 31 DEGs were selected. Actin was used as the inner reference gene. Primer3Plus (http://www.primer3plus.com/cgibin/dev/primer3plus.cgi) was used to select gene-specific primers (Table S1). qRT-PCR was performed using the SYBR R Premix Ex Taq TM II kit (Takara Biomedical Technology (Beijing) Co., Ltd., China) on an ABI 7500 Fast Real-Time PCR system (Applied Biosystems, USA) following the manufacturer's instructions, with three technical replicates. qRT-PCR results were calculated using the 2 − Ct method.

Investigation of drought tolerance of Tartary buckwheat
To determine the drought tolerance of Tartary buckwheat, we measured the root length of Tartary buckwheat seedlings after treatment with 0%, 10%, 20%, and 30% PEG. This method has been widely used to induce a water deficit in plants ( Abdel-Ghany et al., 2020). The results are shown in Fig. 1. The seedlings treated with 10% PEG grew less root hair compared with those treated with 0% PEG; however, the root length showed no significant difference between the two treatments. The roots were dramatically shortened after treatment with 20% PEG compared to those treated with 0% PEG and 10% PEG. The roots did not grow well after treatment with 30% PEG, leading to the loss of approximately 50% of the seedlings. Seedlings treated with 20% PEG showed significant reductions in growth and root length, therefore, we used a 20% PEG treatment for the drought stress experiment.

Physiological changes of Tartary buckwheat seedlings under drought stress
We measured five physiological indexes involved in the drought stress response, namely MDA content, H 2 O 2 content, CAT activity, SOD activity, and POD activity. These indexes were measured at four time points after being treated with 20% PEG and they showed significant alterations (Fig. 2). The H 2 O 2 content was significantly increased after 1 h treatment, but remained nearly unchanged after 3 h and 6 h treatments. The MDA content was significantly increased after 1 h treatment and was significantly decreased afterwards. CAT activity showed an opposite trend to that of the H 2 O 2 content, which was significantly decreased after 1 h treatment, but remained almost unchanged after 3 h and 6 h treatments. The SOD and POD activity showed an opposite trend to that of the MDA content, with a significant decrease after 1 h treatment, and a significant increase afterwards.

High-throughput RNA-Seq of Tartary buckwheat seedlings under drought stress
We performed high-throughput RNA-Seq after 1, 3, and 6 h PEG treatments to obtain the transcriptome dynamics of Tartary buckwheat seedlings after drought tolerance. The 0 h treatment was the control. Three biological replicates were included for each treatment. We performed Pearson's rank correlation analysis to evaluate the repeatability and reproducibility of the transcriptome data. The values of Pearson R between the two samples from the same biological replicates were higher than those of different biological replicates (Fig. S1), indicating that the biological replicates had good repeatability. We obtained 44,303,640 to 62,033,328 raw reads and 42,564,428 to 59,376,546 clean reads for each library (Table 1). The quality of sequencing was high with the Q30 ranging from 92.53% to 93.29%. The genome data of Tartary buckwheat (Zhang et al., 2017) was used as the reference data for mapping, and we successfully mapped 74.49% to 75.45% of the clean reads to the predicted coding sequences (Table 1). A total of 27,490 genes were identified, with 24,433,24,462,24,508,24,271,24,680,24,531,24,617,24,531,24,398, 24,694, 24,553, and 24,263 genes identified in the libraries of PEG0h-1, PEG0h-2, PEG0h-3, PEG1h-1, PEG1h-2, PEG1h-3, PEG3h-1, PEG3h-2, PEG3h-3, PEG6h-1, PEG6h-2, and PEG6h-3, respectively (Table 1). We mapped 25,895, 25,982, 25,957, and 26,033 genes in the treatment of PEG0 h, PEG1 h, PEG3 h, and PEG6 h, respectively ( Fig. 3A). We also performed a principal component analysis (PCA) and the results showed that samples were similar from the same time point after PEG treatment; whereas samples from the time point after the PEG treatment were not (Fig. 3B).

Analysis of DEGs of Tartary buckwheat seedlings under drought stress
We   The remaining genes were differentially regulated at one time point or showed no obvious expression patterns after PEG treatment (Fig. 4). It was worth to mention that circadian affect gene regulation. Previous study has suggested that more than one third expected DEGs were classified as clock-controlled genes comparing Arabidopsis sampled at time 0, 0.5, and 1 h (Hsu & Harmer, 2012). Thus, some fraction of the DEGs identified in our study might not be drought inducible genes, but circadian regulated genes. The top 20 URGs and DRGs were analyzed after being exposed to drought stress for 1 h, 3 h, or 6 h (Tables 2, 3 and 4), among which many stress responsive genes were identified, including FtPinG0005419000.01 (annotated to dehydrin Rab18) and FtPinG0009412200.01 (annotated to carotenoid oxygenase). Interestingly, 5 LEA proteins (FtPinG0002083100.01, FtPinG0005679700.01, FtPinG0004425400.01, FtPinG0001202200.01, and FtPinG0000702400.01) were up-regulated at the third time point of the PEG treatments, suggesting the LEA proteins played crucial roles in response to drought stress. We found that that genes related to the oxidationreduction process and reactive oxygen species biosynthesis were suppressed at the 1 h PEG treatment. These genes included two genes encoding peroxidase (FtPinG0007824100.01 and FtPinG0003282600.01) and one gene encoding allene oxide synthase (FtPinG0002376300.01), suggesting the biosynthesis of the protective enzyme, POD, was suppressed at a molecular level.

Functional category of DEGs of Tartary buckwheat under drought stress
We performed a hierarchical cluster analysis based on the gene expression pattern, with gene function enrichment based on GO annotation, to determine the biological function of  Table S2). C1 included 352 genes that were mildly down-regulated at 1 h and 3 h PEG, and were seriously down-regulated at 6 h PEG. A total of 11 GO biological processes were significantly enriched in this cluster, most of which were related to H 2 O 2 metabolism and catabolism, reactive oxygen species metabolism, and oxidative stress responses. C2 included 189 genes that were up-regulated after PEG treatment; however, there was no significantly enriched GO biological process in this cluster. C3 included 87 genes that were down-regulated at 1 h and 3 h PEG, but up-regulated at 6 h PEG. A total of 14 GO biological processes were significantly enriched in this cluster, most of which were related to cell wall biogenesis, including xyloglucan metabolism, cell wall polysaccharide metabolism, and hemicellulose metabolism. C4 included 291 genes that were up-regulated at 1 h and 3 h PEG, but down-regulated at 6 h PEG, which showed an opposite expression pattern compared to C3. A total of four GO biological processes were significantly enriched in this cluster, and were related to lipid localization and phospholipid transport. C5 included 214 genes that were up-regulated at 1 h and 3 h PEG, but down-regulated at 6 h PEG, which showed a similar expression pattern to C4. A total of 25 GO biological processes were significantly enriched in this cluster, most of which were related to the regulation of candidate processes, such as regulation of gene expression and regulation of biosynthetic process. C6 included 57 genes that were down-regulated at 1 h PEG, but slightly up-regulated at 3 h and 6 h PEG  We also analyzed the metabolic pathways based on KEGG annotation (Fig. 6). At 1 h PEG treatment, 200 of the 213 DEGS were annotated to 30 KEGG pathway, in which two pathways, ''Carotenoid biosynthesis'' and ''Plant hormone signal transduction'', were significantly enriched (Qvalue <0.05, Fig. 6A). At 3 h PEG treatment, 702 of 734 DEGS were annotated to 81 KEGG pathways, in which six pathways were significantly enriched, namely ''Glutathione metabolism'', ''Protein processing in endoplasmic reticulum'', ''alpha-Linolenic acid metabolism'', ''Carotenoid biosynthesis'', ''Biosynthesis of secondary metabolites'', and ''Phenylpropanoid biosynthesis'' (Qvalue <0.05, Fig. 6B). At 6 h PEG treatment, 634 of the 697 DEGS were annotated to 86 KEGG pathways, in which nine pathways were significantly enriched. The pathways were: ''Biosynthesis of secondary metabolites'', ''Circadian rhythm -plant'', ''alpha-Linolenic acid metabolism'', ''Phenylpropanoid biosynthesis'', ''Metabolic pathways'', ''Flavonoid biosynthesis'', ''Protein processing in endoplasmic reticulum'', ''Photosynthesis -antenna proteins'', and ''Carotenoid biosynthesis'' (Qvalue <0.05, Fig. 6C).

The involvement of ABA in response to drought stress in Tartary buckwheat
The carotenoid biosynthesis pathway was significantly enriched at all time points of PEG treatment, based on the KEGG results. Three DEGs were included in this pathway, namely two PSY genes (FtPinG0005737800.01 and FtPinG0004637900.01) and one BCH gene (FtPinG0006960700.01) (Fig. 7). PSY encodes phytoene synthase that catalyzes geranylgeranyl diphosphate to form phytoene, which is the rate-limiting enzyme in the carotenoid biosynthetic pathway. BCH encodes β-carotene hydroxylase of the P-450 monooxygenase family that converts β-carotene to zeaxanthin by a two-step reaction (Ruiz-Sola & Rodriguez-Concepcion, 2012). We determined the contents of phytoene and zeaxanthin to examine whether the carotenoids took a part in a drought response (Fig.  S2). The result showed that neither phytoene nor zeaxanthin content were significantly altered after drought stress, indicating that the carotenoids may not be involved in drought response.
In higher plants, ABA is synthesized from the cleavage of carotenoid precursors (9cis-violaxanthin and 9 -cis-neoxanthin), which is likely the key regulatory step in the ABA biosynthetic pathway. The cleavage reaction is catalyzed by 9-cis-epoxycarotenoid dioxygenase (NCED), and produces xanthoxin, which can be converted into ABA via ABA-aldehyde (Chernys & Zeevaart, 2000). It is possible that ABA was involved in drought response rather than carotenoids. Interestingly, the plant hormone signal transduction    KEGG pathway was one of the most enriched pathways in our results. Fourteen genes related to ABA biosynthesis and signal transduction were differentially expressed in our dataset (Fig. 7)   We measured the ABA content of the Tartary buckwheat seedlings after drought stress and the result of the physiological changes was in accordance with the molecular changes; the ABA content was significantly increased after PEG treatment (Fig. 7).

Transcription Factors (TFs) in response to drought stress in Tartary buckwheat
Two hundred fourteen DEGs in C5 were related to the regulation processes, such as regulation of gene expression and biosynthetic process, based on the hierarchical cluster and the enriched GO biological processes, which is indicative of TFs' role in the drought stress response. We analyzed the differentially expressed TFs in the dataset. We found that a total of 174 TFs belonging to 36 TF families were identified as DEGs, among which the most abundant TF families were WRKY (27), NAC (14), MYB (13), C3H (12), bZIP (10), AP2-EREBP (8), DBP (7), HSF (7), Orphans (7), bHLH (6), C2H2 (6), and MYB-related (6) (Fig. 8). We then analyzed the expression patterns of the top 6 families and the genes in each TF family had a different expression pattern. Overall, most of the TFs were clustered to C5 (59 TFs), followed by C1 (38 TFs), C4 (32 TFs), and C2 (30 TFs) (Table S2), indicating their potential role on the regulation of definite biological processes. For example, among the 27 WRKY TFs, 10 TFs were clustered to C5, suggesting that they may function in the regulation of biosynthetic process; seven TFs were clustered to C4, suggesting they may function in lipid localization and phospholipid transport; three TFs were clustered to C1, suggesting that they may function in H 2 O 2 metabolism and catabolism, reactive oxygen species metabolism, and oxidative stress response; and three TFs were clustered to C3, suggesting they may function in cell wall biogenesis (Figs. 5 and 8). In addition, some TFs, such as FtPinG0002173200.01, FtbZIP83, FtbZIP5, FtPinG0007618600.01, and FtPinG0008274300.01, were homologous to the reported genes within the regulatory network of the drought stress response (Fujita et al., 2004;Lee et al., 2010;Li et al., 2019;Li et al., 2020;Trivedi, Gill & Tuteja, 2016).

DISCUSSION
Physiological changes were highly consistent with molecular changes of Tartary buckwheat seedlings under drought stress ROS accumulate in high amounts when a plant survives oxidative damage or various stresses, such as drought, leading to the damage of the membrane protein and lipids. H 2 O 2 and MDA content increase rapidly under these conditions and the ROS scavenger enzymatic system, including the activities of SOD, APX, and CAT, are induced as an universal response (Mohammadkhani & Heidari, 2007). We measured five widely used physiological indexes in drought stress response, namely MDA content, H 2 O 2 content, CAT activity, SOD activity, and POD activity, after 20% PEG treatment. The results showed that the physiological changes were highly consistent with the molecular changes based on the hierarchical cluster and GO enrichment (biological process) analyses (Figs. 2 and 5, and Table S2). The H 2 O 2 and MDA content were significantly increased in the 1 h PEG treatment, whereas the activity of the protective mechanisms was decreased. Genes related to H 2 O 2 metabolism and catabolism, and ROS metabolism, were clustered to C1 and were down-regulated at this time point, suggesting that ROS accumulation and membrane damage immediately occurred after drought stress in Tartary buckwheat. In the 3 h and 6 h PEG treatment, the activity of the protective enzymes significantly increased, but the MDA content decreased, indicating that the ROS scavenger system was induced against drought stress. The related genes were clustered to C3 and C4, which were associated with cell wall biogenesis and metabolism and lipid localization and phospholipid transport. This suggests that the repair of the ROS and damage to the membrane occurred after the damage to the ROS and membrane.

LEA proteins were involved in drought stress response in Tartary buckwheat
LEA proteins are a type of low molecular weight protein induced by various abiotic stresses, such as drought, high temperature, and cold (Hong-Bo, Zong-Suo & Ming-An, 2005). A number of studies have reported that genes encoding LEA proteins are induced to maintain the stability of the membranes and proteins and to provide detoxification and alleviation of cellular damage under conditions of dehydration (Roychoudhury, Paul & Basu, 2013;Shi et al., 2020;Shinozaki & Yamaguchi-Shinozaki, 2007;Tunnacliffe & Wise, 2007). In a recent transcriptome comparison of drought-resistant and drought-sensitive sorghum genotypes, six of the 25 top-induced genes in drought-resistant genotypes encoded LEA proteins (Abdel-Ghany et al., 2020). Our results were similar (Tables 2, 3 Roychoudhury, Paul & Basu, 2013;Shinozaki & Yamaguchi-Shinozaki, 2007). Blue and pink rounded boxes represent genes in the ABA-dependent and ABA-independent pathway, respectively. Rounded boxes with no border indicate genes were consistent with the reported studies, whereas rounded boxes with red dotted lines indicate genes reported in previous studies were not identified in our study.
Full-size DOI: 10.7717/peerj.11136/ fig-10 4). Twenty-nine up-regulated genes were identified, out of which five genes were annotated to LEA proteins. Among these, four genes were up-regulated at all time points after PEG treatment, including two genes homologues to LEA protein D-29, one gene homologous to LEA protein 2, and one gene homologous to LEA protein 46. Their expression levels were high, with an absolute fold change of 20 to 90 times, suggesting that LEA proteins were involved in the drought stress response in Tartary buckwheat as well.

Transcriptional regulatory network in ABA-dependent and ABA-independent pathways under drought stress of Tartary buckwheat
Plants have two regulatory pathways, the ABA-dependent pathway and the ABAindependent pathway, to manage drought responsive genes (Roychoudhury, Paul & Basu, 2013;Shinozaki & Yamaguchi-Shinozaki, 2007). To better understand the ABAdependent and ABA-independent pathways of Tartary buckwheat under drought stress, a transcriptional regulatory network was constructed based on the previously reported regulatory network and data obtained from our study (Fig. 10). The ABA-dependent pathway under drought stress begins with ABA biosynthesis. ABA biosynthesis begins from the cleavage of carotenoid precursors catalyzed by NCED, Clp protease regulatory subunit that is induced by water stress in Arabidopsis thaliana (Nakashima et al., 2010). We identified a homolog of ERD1 (FtPinG0001990100.01) which was clustered to C2 based on its expression pattern (Table S2 and Fig. 5). Interestingly, some NAC and bZIP TFs were also clustered to C2, suggesting they may be co-expressed and in the upstream of ERD1 (Table S2, Figs. 5, 8 and 10).
The ABA-dependent and ABA-independent pathways were found to participate in the regulation of drought stress of Tartary buckwheat.

CONCLUSION
We investigated the physiological changes and the gene expression changes in a time-course manner under drought stress simulated by 20% PEG treatment. A total of 1,190 DEGs were identified and the genes encoding LEA proteins were listed on the top up-regulated DEGs. All DEGs were grouped into six clusters, in which genes showed definite expression patterns and were involved in specific biological processes based on GO annotation. Further analyses of the ABA and TFs revealed they were also involved in the drought stress response in Tartary buckwheat. We proposed ABA-dependent and ABA-independent pathways in the regulation of drought stress of Tartary buckwheat. This is the first study using a large-scale sequencing method to unravel the transcriptomic changes under drought stress in Tartary buckwheat, which identified massive genes and a gene regulatory network in response to drought stress. Our study provides candidate genes for further functional studies.