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

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₂O₂ solution, rinsed three times with double distilled water (ddH₂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₂O₂ 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.

Reads mapping and analysis of differentially expressed genes (DEGs)

Raw reads were filtered using Trimmomatic-0.38 (Bolger, Lohse & Usadel, 2014) with the following settings: trimmomatic-0.38.jar PE -threads 30 input_forward.fq.gz input_reverse.fq.gz output_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:140. Quality control of all of the acquired clean reads was conducted using FastQC (https://github.com/s-andrews/FastQC) with the following parameters: fastqc -t 30 *.paired. Then the clean reads were mapped to the Tartary buckwheat genome data using Hisat2 (Kim, Langmead & Salzberg, 2015) with the following parameters: hisat2 -x genome -1 output_forward_paired -2 output_reverse_paired -S *.sam. The FPKM value for each gene was calculated using cufflinks (Trapnell et al., 2010). To evaluate the replicate reproducibility, PCA and Pearson correlation analyses were performed in R language (https://www.R-project.org/).

DEGs were identified with the DESeq2 package in R language, with the false discovery rate (padj) <0.05 and absolute value of Log₂Ratio ≥ 1 as the thresholds. Gene annotation was performed by a local BLASP against the NR database, with the threshold of 1e−5. GO items were enriched by agriGO (http://bioinfo.cau.edu.cn/agriGO/) and the KEGG pathway was enriched by KOBAS 3.0 (http://kobas.cbi.pku.edu.cn/kobas3/?t=1). The significantly enriched GO items and KEGG pathways were visualized in R language.

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/cgi-bin/dev/primer3plus.cgi) was used to select gene-specific primers (Table S1). qRT-PCR was performed using the SYBR® Premix Ex Taq™ 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₂O₂ 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₂O₂ 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₂O₂ 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 analyzed DEGs using the pair-wise comparison of samples after different treatment times with 0 h PEG as the control to illustrate the transcriptome changes of Tartary buckwheat seedlings under drought stress. One thousand one hundred-ninety genes were up or down regulated by drought stress, among which 177, 558, and 315 DEGs were up-regulated in 1 h vs. 0 h, 3 h vs. 0 h, and 6 h vs. 0 h, whereas 36, 176, and 382 DEGs were down-regulated in 1 h vs. 0 h, 3 h vs. 0 h, and 6 h vs. 0 h (Fig. 4 and Table S2). Of these DEGs, 86 genes were co-up-regulated at 1 h vs. 0 h, 3 h vs. 0 h, and 6 h vs. 0 h, whereas 6 genes were co-down-regulated at 1 h vs. 0 h, 3 h vs. 0 h, and 6 h vs. 0 h. There were 53, 115, and 6 genes up-regulated in the comparisons of 1 h vs. 0 h, 3 h vs. 0 h, 3 h vs. 0 h, and 6 h vs. 0 h, and 1 h vs. 0 h and 6 h vs. 0 h, respectively, whereas 14, 77, and 0 genes were down-regulated in the comparisons of 1 h vs. 0 h and 3 h vs. 0 h, 3 h vs. 0 h and 6 h vs. 0 h, and 1 h vs. 0 h and 6 h vs. 0 h, respectively. 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 oxidation–reduction 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 the drought-responsive-genes. Six clusters (C1-C6) were identified (Fig. 5 and 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₂O₂ 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 treatment; however, there was no significantly enriched GO biological processes in this cluster.

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 (9′-cis-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), including 4 DEGs homologous to NCED (FtPinG0009412200.01, FtPinG0006853200.01, FtPinG0003131500.01, and FtPinG0000246400.01), 2 DEGs homologous to AREB (FtPinG0002143600.01 and FtPinG0003196200.01), 5 DEGs homologous to PP2C (FtPinG0007629100.01, FtPinG0009574600.01, FtPinG0004850700.01, FtPinG0002889200.01, and FtPinG0006346100.01), 2 DEGs homologous to PYR (FtPinG0001214600.01 and FtPinG0007802700.01), and 1 DEG homologous to OST1 (FtPinG0003981600.01). Interestingly, 15 of the 17 DEGs showed highly consistent expression patterns and were classified to C5, with the exception of FtPSY2 (C4) and FtPYR1 (C1) (Table S2). Among these, two AREB homologs, FtPinG0002143600.01 and FtPinG0003196200.01, have been identified and named as FtbZIP83 and FtbZIP5 respectively in Tartary buckwheat previously (Li et al., 2019; Li et al., 2020).

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₂O₂ 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, Vallabhaneni & Wurtzel, 2008; Li et al., 2019; Li et al., 2020; Trivedi, Gill & Tuteja, 2016).

Confirmation of the transcriptome data by qRT-PCR

The expression patterns of 31 DEGs were verified by qRT-PCR, including five LEA protein encoding genes, seven genes in ABA biosynthesis (FtNCED1, FtNCED2, FtNCED3, FtNCED4, FtPSY1, FtPSY2, and FtBCH), 16 genes in ABA signal pathway (FtABRE1, FtABRE2, FtPP2C1, FtPP2C2, FtPP2C3, FtPP2C4, FtPP2C5, FtPYR1, FtPYR2, FtOST1, FtRD29B, FtRD26, FtDREB1D, FtHB1, FtHB2, and FtRD22), and three TFs in the ABA-independent pathway (FtDREB2, FtRD19, and FtERD1). The results showed that the expression patterns obtained by qRT-PCR were highly correlated with those obtained by transcriptome data, with a Pearson R correlation ranging from 0.76 to 0.98 (Fig. 9), suggesting that our transcriptome results were reliable.

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₂O₂ 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₂O₂ 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₂O₂ 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₂O₂ 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 and 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 ABA-independent pathway, to manage drought responsive genes (Roychoudhury, Paul & Basu, 2013; Shinozaki & Yamaguchi-Shinozaki, 2007). To better understand the ABA-dependent 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, which is the key regulatory step in the ABA biosynthetic pathway (Chernys & Zeevaart, 2000; Seo & Koshiba, 2002). We found that four genes (FtNCED1, FtNCED2, FtNCED3, and FtNCED4) encoding NCED in ABA biosynthesis and three genes (FtPSY1, FtPSY2, and FtBCH) encoding key enzymes (PSY and BCH) in carotenoid biosynthesis were up-regulated (Fig. 7). All of these genes were found to be induced by drought stress or other abiotic stresses in previous studies (Du et al., 2010; Frey et al., 2012; Li, Vallabhaneni & Wurtzel, 2008). There are at least five signal pathways that regulate drought stress in the ABA-dependent pathway in plants. The primary pathway is the ABRE2 mediated signal pathway (Roychoudhury, Paul & Basu, 2013; Shinozaki & Yamaguchi-Shinozaki, 2007). We identified two genes homologous to ABRE2 (FtbZIP83 and FtbZIP5), and eight genes (FtPYR1, FtPYR2, FtPP2C1, FtPP2C2, FtPP2C3, FtPP2C4, FtPP2C5, and FtOST1) in the upstream of ABRE2 as DEGs (Fig. 7). Among these, FtbZIP83 and FtbZIP5 were reported to improve drought/salt tolerance via an ABA-mediated pathway in Transgenic Arabidopsis recently (Li et al., 2019; Li et al., 2020). ABRE2 binds to the DRE/CRT element in the promoter of RD29B and RD20A to activate their expression in response to stress (Yamaguchi-Shinozaki & Shinozaki, 2005). We identified a gene (FtPinG0004870600.01) homologous to RD29B that was induced under drought stress in Tartary buckwheat, however, the homolog of RD20A was not among the identified DEGs (Fig. 10). The ABRE2 mediated signal pathway in addition to the MYB2 and MBC2, RD26, DREB1D, and HB mediated signal pathways are also important constituents in plant ABA-dependent pathway (Roychoudhury, Paul & Basu, 2013; Shinozaki & Yamaguchi-Shinozaki, 2007). We found four homologs of RD26 (FtPinG0002173200.01), DREB1D (FtPinG0008274300.01), and HB (FtPinG0001748800.01 and FtPinG0008157500.01), all of which were up-regulated under drought stress (Table S2). Though we found no homologs of MYB2 or MBC2, one gene (FtPinG0002802300.01) homologous to its downstream gene, RD22, was identified in our data. It was down-regulated after drought stress and exhibited a different expression pattern compared to other genes in the ABA-dependent pathway (Table S2 and Fig. 5). Sixteen of the 20 DEGs mentioned above showed similar expression patterns and were up-regulated at 1 h and 3 h PEG, but down-regulated at 6 h PEG, and clustered to C5 (Table S2, Fig. 5 and Fig. 7), with the exception of FtPSY2, FtPYR1, FtPinG0008274300.01 and FtPinG0002802300.01. These results confirmed that genes in ABA-dependent regulatory pathway are co-regulators in response to drought stress in Tartary buckwheat. In addition, changes of ABA content corresponded with the molecular changes, and confirmed the role of ABA in response to drought stress and tolerance in buckwheat (Fig. 7).

There are at least two signal pathways, namely the DREB2 mediated signal pathway and the NAC and bZIP TFs mediated signal pathway, to regulate drought stress in the plant ABA-independent pathway (Fu et al., 2016; Roychoudhury, Paul & Basu, 2013; Shinozaki & Yamaguchi-Shinozaki, 2007). DREB2 belongs to the AP2-EREBP TF family and may transcriptionally activate the expression of RD19 and RD29A (Lee et al., 2010). We identified the homologues of DREB2 (FtPinG0007618600.01) and RD19 (FtPinG0002253700.01) but were unable to identify the homolog of RD29A (Table S2 and Fig. 10). Massive NAC and bZIP family TFs have been reported to function in the ABA-independent pathway by activating ERD1, which encodes a chloroplast-targeted 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.