Biotic and abiotic stress-responsive genes are stimulated to resist drought stress in purple wheat

Journal of Integrative Agriculture 2020, 19(1): 33–50 Available online at www.sciencedirect.com

ScienceDirect

RESEARCH ARTICLE

Biotic and abiotic stress-responsive genes are stimulated to resist drought stress in purple wheat LI Xiao-lan1, 3, LÜ Xiang1, WANG Xiao-hong1, PENG Qin2, ZHANG Ming-sheng1, REN Ming-jian2 1

School of Life Sciences, Guizhou University/State Engineering Technology Institute for Karst Desertification Control/Key Laboratory of Plant Resources Conservation and Germplasm Innovation in Mountainous Region, Ministry of Education, Guiyang 550025, P.R.China 2 School of Agriculture, Guizhou University/Guizhou Sub-Center of National Wheat Improvement Center, Guiyang 550025, P.R.China 3 Special Key Laboratory of Microbial Resources and Drug Development from Higher Education Institution of Guizhou Province/ Research Center for Medicine and Biology, Zunyi Medical University, Zunyi 563006, P.R.China

Abstract Triticum aestivum L. cv. Guizi 1 (GZ1) is a drought-tolerant local purple wheat cultivar. It is not clear how purple wheat resists drought stress, but it could be related to anthocyanin biosynthesis. In this study, transcriptome data from droughttreated samples and controls were compared. Drought slightly reduced the anthocyanin, protein and starch contents of GZ1 grains and significantly reduced the grain weight. Under drought stress, 16 682 transcripts were reduced, 27 766 differentially expressed genes (DEGs) were identified, and 379 DEGs, including DREBs, were related to defense response. The defense-response genes included response to water deprivation, reactive oxygen, bacteria, fungi, etc. Most of the structural and regulatory genes in anthocyanin biosynthesis were downregulated, with only TaDFR, TaOMT, Ta5,3GT, and TaMYB-4B1 being upregulated. TaCHS, TaF3H, TaCHI, Ta4CL, and TaF3’H are involved in responses to UV, hormones, and stimulus. TaCHS-2D1, TaDFR-2D2, TaDFR-7D, TaOMT-5A, Ta5,3GT-1B1, Ta5,3GT-3A, and Ta5,3GT-7B1 connect anthocyanin biosynthesis with other pathways, and their interacting proteins are involved in primary metabolism, genetic regulation, growth and development, and defense responses. There is further speculation about the defense-responsive network in purple wheat. The results indicated that biotic and abiotic stress-responsive genes were stimulated to resist drought stress in purple wheat GZ1, and anthocyanin biosynthesis also participated in the drought defense response through several structural genes. Keywords: transcriptome, purple wheat, drought, anthocyanin, differentially expressed genes, defense response, stress

1. Introduction Received 29 September, 2018 Accepted 19 February, 2019 LI Xiao-lan, Tel: +86-851-28642444, E-mail: [email protected]; Correspondence REN Ming-jian, Tel: +86-851-83855894, E-mail: [email protected] © 2020 CAAS. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http:// creativecommons.org/licenses/by-nc-nd/4.0/). doi: 10.1016/S2095-3119(19)62659-6

Arid and semiarid soils account for approximately onethird of the global land area. Drought is a major stress for plants, and the drought-induced imbalance of such stresses as reactive oxygen species, osmotic stress, and water loss severely affect the growth and development of plants and reduce crop yields (Fang et al. 2017), which

34

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

seriously threatens human survival and population growth. However, plants have a series of defense mechanisms to resist drought stress injury and adapt to environmental changes for survival. Plants resist drought stress through a variety of defense mechanisms. To avoid or reduce drought stress injury, plants regulate stomata opening and closing, osmotic homeostasis, antioxidant systems, hormone signaling, miRNA family, light protection, and metabolic pathways (Mahajan and Tteja 2005; Shi et al. 2018). This regulatory network involves multiple genes and pathways. Therefore, drought stress alters the expression levels of a large number of genes and the biosynthesis of metabolites. Anthocyanin is a secondary metabolite that is beneficial to human health and contributes to stress resistance in plants. The anthocyanin has anticancer, anti-inflammatory, and antitumor activities in humans and attracts insect pollinators and provide antioxidant potential in plants (Li et al. 2017). Anthocyanin biosynthesis has been studied in many plants, and the basic molecular mechanism of its synthesis is clear. Anthocyanin biosynthesis is controlled by structural and regulatory genes. The structural genes include phenylalanine ammonia-lyase (PAL), 4-coumarateCoA ligase (4CL), chalcone synthase (CHS), chalcone isomerase (CHI), flavanone 3-hydroxylase (F3H), flavonoid 3’-hydroxylase (F3’H), flavonoid 3’5’-hydroxylase (F3’5’H), dihydroflavonol 4-reductase (DFR), anthocyanidin synthase (ANS), glucosyltransferase (GT), methyltransferase (MT), and acyltransferase (AT). Among these genes, CHS and ANS encode the key enzymes (Springob et al. 2003; Sasaki et al. 2013). The regulatory genes encode MYB, bHLH, and WD40 transcription factors (Jaakola 2013). Under drought stress, the expression levels of TaCHS, TaCHI, TaF3H, TaDFR, and TaANS in wheat leaves were upregulated, the anthocyanin content was increased, the plants had stronger drought resistance, but there were differences among different cultivars (Ma et al. 2014). As a part of the plant antioxidant system, anthocyanin biosynthesis was affected by drought stress. Purple wheat grains biosynthesize and accumulate anthocyanin (Jiang et al. 2018). Triticum aestivum L. cv. Guizi 1 (GZ1) is a local purple wheat cultivar that has stable yield and can resist drought stress in agricultural production. Therefore, anthocyanin biosynthesis in GZ1 may be involved in drought stress resistance. To test the above hypothesis, we used RNA sequencing to explore the molecular mechanisms of drought resistance, which involve multiple genes and regulatory networks. RNASeq technology is a reliable method to obtain gene coding sequences, their relative expression levels, and annotation information through database comparison (Marguerat and Bahler 2010; Jazayeri et al. 2015). In this study, 40 mmol L–1 PEG 6000 solution was used to simulate drought stress

on GZ1. Then, differentially expressed genes (DEGs) were identified by comparing the transcriptome data between controls and treated samples. These defense-responsive genes were filtered by P-value, DEGs annotation, and classification. The differential expression profiles of structural and regulatory genes involved in anthocyanin biosynthesis were analyzed, and the accuracy of transcriptome expression levels was verified by qRT-PCR. Moreover, protein-protein interaction analysis was used to determine the structural enzymes that interacted with other proteins under drought stress. This information will contribute to the selection of stress-resistant purple wheat and our understanding of the molecular mechanisms associated with stress tolerance.

2. Materials and methods 2.1. Plant growth and stress treatment Triticum aestivum L. cv. Guizi 1 (GZ1, 2n=6x=42, AABBDD) seeds were collected from Guizhou Sub-Center at the National Wheat Improvement Center, Guiyang, China. Seeds were soaked in distilled water until germination and transplanted to soil. The potting soil consisted of humus soil and yellow soil at a ratio of 1:1 (v:v), with 15 kg of soil and five plants per pot, and six pots for each treatment. A biological replicate consisted of two pots, and six pots corresponded to three biological replicates. Seedlings were grown in an environmental chamber with a 16-h/8-h light/dark photoperiod at 22°C/18°C, an irradiance of 150 μmol m−2 s−1 and 60% humidity. When the plants grown to 15 days post anthesis (dpa), the drought group was treated with 40 mmol L–1 PEG 6000 solution (1 L per pot), and the control group was treated with distilled water (1 L per pot). After 24 h, six or seven spikes were taken from each biological replicate, and three kernels were taken from the middle part of each spike. The grains of each biological replicate were mixed and ground under liquid nitrogen, then weighed to approximately 50 mg for extraction of total RNA.

2.2. RNA isolation, cDNA library, and sequencing Total RNA was extracted with the RNAiso Plus Kit (TaKaRa, Japan) according to the manufacturer’s instructions, and its concentration and integrity were determined with an Agilent 2100 Bioanalyzer System with Agilent RNA 6000 Nano Kit (Agilent Technologies Inc., USA). The concentration of total RNA was higher than 1 000 ng μL–1, and the RNA integrity number (RIN) was higher than 7, which indicated good total RNA integrity. A NanoDrop 2000 spectrophotometer (Thermo Scientific, USA) was used to detect the OD values; the OD260/OD280 was ≥1.8, and the OD260/OD230 was ≥1.8 (the specific measurement data are detailed in Appendix A and

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

Table 1). High-quality mRNA from three biological replicates were pooled and used for transcriptome sequencing. After fragmentation of oligo-dT magnetic bead-enriched mRNA, cDNA was synthesized for library construction (NEBNext® UltraTM II RNA Library Prep Kit, Gene Co., Ltd., China). The Illumina HiSeq 4000 platform was used for sequencing (BGI Technologies Co., Ltd., Wuhan, China). High-quality clean data with a Q30>93% were obtained by removing linker sequence and low-quality data. The average output from each sample was 6.67 GB. The dataset was submitted to the NCBI database (accession number: SRP148475).

2.3. Genome mapping, novel transcript prediction, and differentially spliced gene detection Clean reads were mapped to the wheat reference genome sequence (iwgsc_refseqv1.0_High Conf_CDS_2017Mar13. fa and iwgsc_refseqv1.0_High Conf_2017Mar13.gff3, the International Wheat Genome Sequencing Consortium, http://www.wheatgenome.org/) using hierarchical indexing for spliced alignment of transcripts (HISAT; v2.0.4, USA; Kim et al. 2015). The average mapping ratio of each sample was 84.00%, and the uniform alignment ratio between samples showed that the data were comparable between samples. The unmapped clean reads were used for transcript reconstruction using StringTie (v1.0.4, USA) (Pertea et al. 2015), and comparison of the reconstructed transcripts with reference annotation information was performed using Cuffcompare (Cufflink tool, v2.2.1, USA) (Trapnell et al. 2012). The coding potentiality of the novel transcripts was determined using the CPC (Coding Potential Calculator) Software (v0.9-r2, China; Kong et al. 2007). Next, novel coding transcripts were merged with reference transcripts to get a complete reference genome; downstream analysis was based on this reference. Table 1 Summary of Illumina HiSeq 4000 transcriptome sequencing and assembly for Guizi 1 (GZ1) Item1) Total bases (Mb) Raw reads Clean reads Clean reads Q20 (%) Clean reads Q30 (%) Total mapping ratio (%) Total gene number Known gene number Novel gene number Total transcript number Known transcript number Novel transcript number N50 length (bp) 1)

Control

Drought

58.3 45.1 98.39 95.25 76.46 73 294 64 442 8 852 88 995 61 771 27 224 1 945

53.44 44.19 98.52 95.47 87.75 62 771 55 737 7 034 72 313 49 499 22 814 1 844

N50 length, median length of all non-redundant sequences.

35

Differences in the relative abundance of subtypes of the same gene indicated that gene was regulated by splicing. Differentially spliced genes (differential isoform relative abundance between samples) were detected using rMATS (replicate multivariate analysis of transcript splicing, v3.0.9, USA; FDR<0.05, FDR is the false discovery rate; Shen et al. 2014; Wang J et al. 2017). The MATS statistical model was used to calculate the P-value and FDR. A gene showing a difference in isoform ratio between two conditions and an FDR≤0.05 was defined as a significantly differentially spliced gene.

2.4. Differentially expressed genes (DEGs) filtering and their functional annotation Gene expression levels were calculated using RNA-Seq by Expectation Maximization (RSEM) (Li and Dewey 2011). The significance of the difference in gene expression between the control and drought-treated samples was determined by the PossionDis method. The FDR was used to identify the P-value threshold in multiple tests. The DEGs were called when the fold change was greater than 2.00 and the FDR less than 0.001. If the log2fold change was greater than 1 and FDR was less than 0.001, the DEG was considered to be upregulated, and if the log2fold change was less than –1 and FDR was less than 0.001, the DEG was considered to be downregulated. Gene Ontology (GO, http://www.geneontology.org) functional classification and enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG, http:// www.genome.jp/kegg/) biological pathway classification, and enrichment analysis, and the NCBI nonredundant protein database (Nr, http: //www.ncbi.nlm.nih.gov/RefSeq/) comparison and annotation were performed on the DEGs. The detailed analytical methods were described previously (Wang et al. 2018).

2.5. Filtering defense-responsive genes, as well as structural and regulatory genes, in anthocyanin biosynthesis Defense-responsive genes were filtered by using keywords such as response, responsive, and defense response, then candidate genes were identified using GO, KEGG, and Nr annotation. These genes were counted and classified according to their P-value, specific stresses, possible encoded proteins, and homologs in other plants. Using the KEGG and Nr annotations and previous reports, anthocyanin biosynthesis-related genes were selected. The structural genes included TaPAL, Ta4CL, TaCHS, TaCHI, TaF3’5’H, TaF3H, TaDFR, TaANS, TaOMT, Ta3GT, TaUGT, Ta5,3GT, Ta5AT, TaCT (agmatine coumaroyl-transferase),

36

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

and TaHT. The regulatory genes included TaMYB and TaMYC. A list of these genes can be found in Appendix A and Table 2. Their expression levels were calculated using the fragments per kb per million reads (FPKM) value, and the log10FPKM values were used to draw a heatmap of gene expression.

2.6. Quantitative real-time PCR for TaCHS and TaANS Total RNA was extracted and evaluated as described above. First-strand cDNA synthesis was performed with 1 μg of total RNA using the PrimeScript RT Master Mix Kit (TaKaRa, Japan) according to the manufacturer’s recommended procedure (Li et al. 2014). Primers for each gene were designed with Beacon Designer version 8.20 (Premier Biosoft international, Palo Alto, CA, USA). The primers used for amplifying TaCHS, TaANS, and 18S ribosomal RNA (rRNA, GenBank accession no.: AY049040.1) are listed in Appendix B and Table 1. qRT-PCR was conducted on a Thermal Cycler Dice Real-Time System CFX96 Touch (Bio-Rad Laboratories Co., Ltd., USA) using SYBR Premix Ex Taq II (Tli RNAseH Plus; TaKaRa, Japan). The PCR conditions were as follows: (1) predenaturation at 95°C for 3 min; (2) 38 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 45 s; (3) a final extension at 72°C for 4 min. Baseline and threshold cycles (CT) were automatically determined using the CFX Manager Software (Bio-Rad Laboratories Co., Ltd., USA). Relative gene expression with respect to18S rRNA was determined using the 2–ΔΔCT method (Livak and Schmittgen 2001).

2.7. Protein-protein interaction networks in anthocyanin biosynthesis Double index alignment of next-generation sequencing data

(DIAMOND) was used to map the DEGs to the STRING database to obtain the interactions between DEG-encoded proteins using homology with known proteins (von Mering et al. 2005; Buchfink et al. 2015). The interaction networks related to anthocyanin biosynthesis structural gene-encoded proteins were filtered, and the network was created using the Cytoscape Software (v3.5.1, Germany).

2.8. Physiological and biochemical analysis The protein and starch contents in mature grains were determined using an InfraTec 1241 grain analyzer (FOSS, Denmark), and the thousand-grain weight was determined by a Wanshen SC-G automatic test and thousand-grain weight system (Hangzhou Wanshen Detection Technology Co., Ltd., China). The methods for total anthocyanin extraction and anthocyanin content determination were described previously (Abdel-Aal and Hucl 1999).

2.9. Statistical analysis The data were statistically analyzed in Excel 2013 and IBM SPSS Statistics 22.0 (SPSS Inc., USA) using the Duncan method (P<0.05) and a one-way ANOVA mathematical model. The standard error is presented for data performed in triplicate. The data were plotted using Origin 8.0 (OriginLab Corporation, USA) and Excel 2013.

3. Results 3.1. Phenotype observation and physiological and biochemical properties To explore the effects of drought stress on purple wheat grains in the early grain-filling stage, wheat plants were

Table 2 Defense-responsive structural genes in anthocyanin biosynthesis under drought stress Stress response GO:0009411//response to UV

GO:0009416//response to light stimulus GO:0006950//response to stress

GO:0009725//response to hormone

GO:0060918//auxin transport

Structural gene TaCHS-5D1; TaCHS-5D2 TaCHS-2A2; TaCHS-2B1; TaCHS-2B2; TaF3H-2A; TaF3H-2B; TaF3H-2D; TaCHI-5A; TaCHI5B-T1; TaCHI-5B-T2; TaCHI-5D TaCHS-2D1-T1; TaCHS-2D1-T2 TaF3’H-6D2 Ta4CL-6D2 TaCHS-5D1; TaCHS-5D2 TaCHS-2A2; TaCHS-2B1; TaCHS-2B2 TaCHS-2D1-T1; TaCHS-2D1-T2 TaCHS-5D1; TaCHS-5D2 TaCHS-2A2; TaCHS-2B1; TaCHS-2B2 TaCHS-2D1-T1; TaCHS-2D1-T2 TaCHS-5D1; TaCHS-5D2 TaCHS-2A2; TaCHS-2B1; TaCHS-2B2 TaCHS-2D1-T1; TaCHS-2D1-T2

37

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

watered with 40 mmol L–1 PEG6000 solution at 15 dpa. The treated plants did not show significant phenotypic changes and continued to grow. Mature grains were harvested, and the protein, starch, and anthocyanin contents and thousand-grain weight were determined. The results are shown in Fig. 1. The protein, starch, and anthocyanin contents were slightly lower, whereas the thousandgrain weight was obviously lower compared with that of the control (P<0.05). Drought stress during the early grain-filling stage had weak effects on grain nutrition but had a strong influence on grain weight. The grains may have undergone various regulatory mechanisms at the genetic level to resist drought. To explore these regulatory responses, transcriptome sequencing on the grains 24 h after drought treatment was performed.

3.2. GZ1 sequence data and transcriptome assembly Sequencing the control and treated grain transcriptomes using the Illumina HiSeq 4000 platform yielded 45 098 658 and 44 185 084 high-quality clean reads, respectively. Additional sequence data information is presented in Table 1. Through transcriptome assembly, the clean reads generated 88 995 and 72 313 transcripts for the control and treated grains, respectively. After genome mapping, there were 27 224 and 22 814 novel transcripts, and 8 852 and 7 034 novel genes for the control and treated samples, respectively. Comparison of the sequencing data and transcriptome assembly results showed that there were

C

Some effects of drought stress can be inferred based on the differential expression of genes between the control and treated samples. DEGs were filtered (fold change>2.00 and FDR<0.001), annotated, and classified. In total, 27 766 DEGs were obtained, with 53.37% of DEGs having P-values from 1E-4 to 1E-15, 19.51% having P-values from 1E-15 to 1E-30, and 1.93% having a P-value less than 1E-200 (Fig. 2-A). As the P-value decreased, the number of DEGs also gradually decreased, resulting in a pyramid-type expression pattern. When the P-value was less than 1E-200, the number of upregulated DEGs (325, 60.75%) was greater than the number of downregulated DEGs (210, 39.25%). The opposite was observed for the rest of the DEGs (Fig. 2-B and C). DEGs with P-values less than 1E-200 may be especially responsive to drought stress. Although there were few of these genes, most were upregulated, and they were also the most significant DEGs. The DEGs were used for functional annotation by comparison with the Nr, GO, and KEGG databases (Fig. 2-C). In total, 27 526 DEGs (99.14%) were annotated in the Nr database, among which 226 were transcription

D

25

Starch (%)

20

Protein (%) CK

3.3. Functional annotation and classification of the DEGs

15 10 5 0

Drought

CK

Drought

E

CK

Drought

Anthocyanin content (mg kg–1)

B

60 50 40 30 20 10 0

CK

Drought

F 100 90 80 70 60 50 40 30 20 10 0

CK

Drought

Thousand-grain weight (g)

A

more genes and transcripts in the control than in the droughttreated samples. This finding means that drought reduced the number of genes that were transcribed in treated grains, which may result in a decrease in metabolic activity.

40 35 30 25 20 15 10 5 0

*

CK

Drought

Fig. 1 Phenotypic observation and partial physical and chemical parameters of the control (CK) and treated samples. A, the observed phenotype of CK and treated plants. B, the observed phenotype of CK and treated mature grains. C, protein content. D, starch content. E, anthocyanin content. F, thousand-grain weight. The standard error was calculated from triplicate data. * represents significant differences (P<0.05).

38

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

A 1E-75 to 1E-100 1.80% 1E-60 to 1E-75 2.02%

B

1E-100 to 1E-200 2.27% 1E-200 to 0 1.93% 1E-3 to 1E-4 7.76%

1E-45 to 1E-60 3.73%

1E-15 to 1E-30 19.51%

P-value

1E-30 to 1E-45 7.61%

1E-4 to 1E-15 53.37%

1E-200 to 1E-0

210 535 325

1E-100 to 1E-200

416 630 214

1E-75 to 1E-100

382 500 118

1E-60 to 1E-75

446 562 116

1E-45 to 1E-60

852 1 035 183

1E-30 to 1E-45

300

1E-15 to 1E-30

Down

12 420 2 404

419 0

30 000

27 766

27 526 22 028

20 000 15 000

12 407

12 932

Species

Number of genes

25 000

11 887

10 000 4 003

5 000

226

0 Total

Nr

KEGG

Cellular Molecular Biological component Function Process

Annotations

GO and GO KEGG Common annotations

Nr-TF

14 824

1 737 2 156

2 000 4 000 6 000 8 000 10 000 12 000 14 000 16 000 18 000 20 000

Number of genes

Phyllostachys edulis

28

Volvox carteri f. nagariensis

32

Triticum turgidum subsp. durum

34

Cynara cardunculus var. scolymus

34

Emiliania huxleyi CCMP1516

35

Beta vulgaris subsp. vulgaris

40

Triticum monococcum

58

Hordeum vulgare

86

Oryza sativa Indica Group

96

Oryza brachyantha

129

Zea mays

224

Setaria italica

243

Sorghum bicolor

267 428

Oryza sativa Japonica Group

3 798

4 716 5 420

704

1E-3 to 1E-4

D

Up

1 813 2 113

1E-4 to 1E-15

C

Total

621

Other species

3 179

Brachypodium distachyon

3 347

Triticum aestivum

4 717

Triticum urartu

6 656

Hordeum vulgare subsp. vulgare

7 161

Aegilops tauschii 0

1 000 2 000 3 000 4 000 5 000 6 000 7 000 8 000

Number of genes

Fig. 2 Functional annotation and classification of the differentially expressed genes (DEGs). A, P-value distribution of the DEGs. B, upregulation and downregulation statistics for the DEGs with different P-value distributions. C, DEG annotation statistics. D, species distribution of sequences. Nr, the NCBI nonredundantprotein database; KEGG, Kyoto encyclopedia of genes and genomes; GO, Gene Ontology; TF, transcription factor.

factors. A total of 22 028 DEGs (79.33%) were annotated in the KEGG database. The annotation from the GO database was as follows: (1) 12 407 DEGs (44.68%) were cellular components, (2) 12 932 (46.57%) were molecular functions, (3) 11 887 (42.81%) were biological processes, and (4) 4  003 (14.41%) were associated with all three aspects. In total, 3 799 DEGs (13.68%) were annotated in both the GO and KEGG databases. The DEGs were comprehensively annotated, which provided a good basis for further study of drought defense mechanisms. According to the functional annotation from the Nr database, the top five plant sources of the DEGs were as follows: (1) 7 161 (25.79%) from Aegilops tauschii, (2) 6 656 (23.97%) from Hordeum vulgare subsp. Vulgare, (3) 4 717(16.99%) from Triticum urartu, (4) 3 347 (12.05%) from

Triticum aestivum, and (5) 3 179 (11.45%) from Brachypodium distachyon, reaching a total of 25 060 (90.25%) (Fig. 2-D). This information might aid the prediction of the evolutionary relationships among GZ1 and other wheat varieties.

3.4. GO and KEGG enrichment of the DEGs To further classify the DEGs to identify defense-responsive genes, GO enrichment analysis was performed for the biological process, cellular component, and molecular function categories. It was found that the number of downregulated genes was greater than the number of upregulated genes in all three categories (Fig. 3). However, many upregulated genes were annotated as being involved in metabolic processes, cellular processes, and responses

39

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

Number of genes

12 000

Up

10 000

Total

Down

8 000 6 000 4 000 2 000

B

N

V

C

M

M

P

N

M

R

T

S

M

S

N

T E

S

E M

M

A

C

O M

C

S

C O M P C

N

S

C B

B

I

R M

D

B

R

M

R

C

D

R

L

L

C

S

G

V

C

0

Biological process

Cellular component

Molecular function

Fig. 3 Gene Ontology (GO) enrichment of the differentially expressed genes (DEGs).

to stimuli in biological processes; cells, cell parts, and organelles in cellular components; and binding and catalytic activity in molecular functions. This information may contribute to further studies on gene function and determining how to narrow down the search for defenseresponsive genes. KEGG enrichment analysis was performed to better understand the possible outcomes of changes in expression of the DEGs. KEGG enrichment was divided into the following five categories: (1) metabolism (12 659, 59.96%), (2) genetic information processing (5 272, 24.97%), (3) environmental information processing (1 296, 6.14%), (4) cellular processes (1 008, 4.77%), and (5) organismal systems (876, 4.15%; Fig. 4-A). The KEGG term “metabolism” was associated with the largest number of DEGs, indicating that the effects of drought stress on metabolism was the greatest. A bubble chart was made according to the enrichment number and adjusted P-value (Q-values). As shown in Fig. 4-B, there was a large enrichment for DEGs in RNA transport, mRNA surveillance, biosynthesis of amino acids, carbon metabolism, and so on. The number of downregulated genes was greater than the number of upregulated genes in these pathways (Fig. 4-C). The protein and starch contents are the most important part of the wheat grain which accounted for approximately 70% of grain weight (Fig. 1). However, genes involved in amino acid biosynthesis, carbon metabolism, and the citrate cycle (TCA cycle) pathways were downregulated under drought stress. These pathways are necessary for protein and starch biosynthesis. The downregulated genes in these pathways may result in a small decrease in protein and starch content, but a significant decrease in thousand-grain weight.

3.5. Defense-responsive DEGs under drought Based on the annotation for biological process, a filter for defense-responsive DEGs induced by drought stress, specifically “defense response” was used. A total of 323 defense-responsive DEGs were obtained and classified according to their P-values (Fig. 5). Among them, 45.99% of the defense-responsive DEGs had P-values from 1E-4 to 1E-15, 20.99% from 1E-15 to 1E-30, and 9.26% less than 1E-200. Additionally, 71 defense-responsive DEGs were upregulated, and 252 were downregulated. The upregulated defense-responsive DEGs were mainly distributed at the two ends of the P-value range; in particular, 27 of them (38.89%) had P-values less than 1E-200, and 12 (16.67%) had P-values from 1E-100 to 1E-200. In contrast, the downregulated defense-responsive DEGs were mainly distributed in ranges with larger P-values; for example, 192 of them (76.19%) had P-values from 1E-4 to 1E-30. The defense-responsive DEGs displayed relatively concentrated P-values. The smaller the P-value is, the greater the difference in the DEG expression levels between the control and treated samples, indicating that the sensitivity of these genes’ response to drought was higher. Genes with P-values less than 1E-200 may be very important for stress defense, especially the upregulated genes. To explore in which defense reactions these related DEGs participated, their biological process annotations were classified, and the results are shown in Appendix C and Figs. 1–7. These defense reactions can be divided into nine categories, including defense, immune, cellular, sugars, light, biotic stress, hormone, abiotic stress, and compound and stimulus responses. These defense reactions involved a wide range of biological process, indicating that drought

40

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

A

Transport and catabolism

1 008

Signal transduction

1 039

Membrane transport

Non-homologous end-joining Nicotinate and nicotinamide metabolism

257

Inositol phosphate metabolism 2 024

Translation

4 951 1 927

Carbohydrate metabolism

Nucleotide metabolism

1 012 905

Energy metabolism

511

Metabolism of other amino acids

507

Glycan biosynthesis and metabolism

486

Environmental adaptation

750

Pantothenate and CoA biosynthesis

Q-value

N-Glycan biosynthesis

Environmental information processing Genetic information processing Metabolism

2.0E-06

Nucleotide excision repair

1.5E-06

ABC transporters

1.0E-06

Biosynthesis of amino acids

5.0E-07

Citrate cycle (TCA cycle) Fructose and mannose metabolism

646

Metabolism of cofactors and vitamins

Metabolism of terpenoids and polyketides

Cellular processes

Organismal systems

840

536

500

Mismatch repair

Global and overview maps

Lipid metabolism

250

Butanoate metabolism

521

Amino acid metabolism

Gene number

C5−Branched dibasic acid metabolism

993

Transcription

Biosynthesis of other secondary metabolites

Endocytosis

1 734

Folding, sorting and degradation

Replication and repair

B

Valine, leucine and isoleucine biosynthesis Carbon metabolism Alanine, aspartate and glutamate metabolism mRNA surveillance pathway Ribosome biogenesis in eukaryotes

338

RNA transport 876 2 000

0.3

0.4

0.5

0.6

Rich factor

4 000 6 000 8 000

Number of genes

C

DEGs number of the most enriched pathway

DEGs number

800 Down

600

Up

400 200

mRNA

surve

RNA

transp ort Endo cytos is illance pathw ay Carbo n meta Biosy bolism nthes is of a mino acids Purin e meta Glyco bolism lysis/G Ribos lu coneo ome b genes iogen is esis in eukary otes Phosp A B C tran hatidy Carbo s linosit porters n fixa ol sig tion in naling photo syste synth m etic o rganis Citrate ms Fructo cycle (TCA se an cycle) d man nose metab Nucle olism otide excisio n repa Pyruv ir ate m Inosit etabo ol pho lism sphate metab N−Gly olism can b iosyn Alanin thesis e, asp DNA artate replic and g ation lutam ate m etabo lism Mism atch re pair Prote asom Fatty Nicoti e acid b nate a iosyn nd nic thesis otinam ide m etabo Butan lism Panto oate m thena etabo Valine te and lism , leuc CoA b ine an iosyn d isole thesis ucine biosy nthes Propa C5−B is noate ranch metab ed dib olism asic a c id me Non− tabolis homo m logou s end −joinin Sulfur g relay syste m

0

Fig. 4 Kyoto encyclopedia of genes and genomes (KEGG) enrichment of differentially expressed genes (DEGs). A, statistics for the KEGG enrichment of DEGs. B, the number and Q-values of the DEGs in the pathway enrichment. C, the upregulated and downregulated DEGs in the pathways.

stress induced a complex network of defense response reactions, which included responses to abiotic and biotic stresses, signaling pathways, metabolism, etc. In addition, these defense-responsive genes encoded more than 123 proteins, for example, interleukin-1 receptor-associated kinase 4, chaperonin GroEL, myb proto-oncogene protein, and peroxidase (more details in Appendix D), which indicated that the proteins involved in these defenseresponsive reactions were also abundant. However, the defense-responsive DEGs that may be the most sensitive to drought stress should be considered. Candidates that had a P-value of less than 1E-200 (29 DEGs) were selected. These defense-responsive genes were involved in responses to reactive oxygen species,

osmotic stress, water deprivation, plant-type hypersensitive reactions, light intensity, hormones, metal ions, bacteria, fungi, etc. (Fig. 6, more detailed information in Appendix E and Table 1). Among them, 27 DEGs (93.10%) were upregulated. Of them 12 DEGs (41.38%) were associated with defense response to fungi, of which 11 were upregulated and only one was downregulated. Drought defense response was also observed, including osmotic stress response (6, 20.69%) and water deprivation response (3, 10.34%). These results suggested that drought stress may induce GZ1 grains to increase their defense responses not only to abiotic stress but also to biotic stress, specifically the defense response to fungi. In addition, the data of significant DEGs (fold change>2.00, FDR<0.001, total

41

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

149 131

160

Up Total Down

120 100

68

80

11 Down Total

1E

-4

-3

to

to

1E

-4

-1 5

Up

1E

-3 0 1E to

-1 5

11

18

7

-4 5 1E to

61

15

17

2

-6 0 to

1E

1E

-4 5

to -6 0

to -7 5 1E

1E

-1 0 1E

1E to 0

-1 0

-7 5

0

0 -2 0

0 to 0 -2 0 1E

13

16 3

2

0

1E

9

5

12

20

7

1E

5

7

1E

19

2

27

-3 0

29

40

1E

60

1E

Number of genes

140

E-value distribution of defense related genes

Fig. 5 Defense-responsive differentially expressed gene (DEG) statistics.

14

13

10

12 11

8 6

Response to reactive oxygen species

Defense response

4

Response to water deprivation

Plant-type hypersensitive response

Response to light intensity

Response to hormone

Defense response, incompatible interaction

3

Response to metal ion

Defense response to bacterium

Response related biological process 1E-200 to 0

Up

Down

Up

Pigment accumulation in tissues in response to UV light

1

Total

1

Down

Up

Down

1

Total

Up

1

Down

Total

Up

Down

Total

Up

Down

Up

1

Down

Up

Down

1

Total

Up

Total

Down

Up

Total

Down

4

3

1

Response to osmotic stress

5

4

3

1

Up

Down

1

Total

1

Up

1

Total

0

3

Total

4 2

5 4

Total

5

Total

6

Down

Number of genes

Total Up Down

12

12

Defense response to fungus

Fig. 6 Defense-responsive differentially expressed genes (DEGs) with a P-value of less than 1E-200.

FPKM≥5) and highly DEGs (fold change>2.00, FDR<0.001,

encapsulating structure that may perform shuttling functions

Max FPKM≥100) were in Appendix E.

in the plastid and nuclear lumen and is cold-inducible, and

Analysis of the function of these genes based on

one pathogenesis-related protein with nuclease activity

biological process annotation allowed us to determine their

located on an intracellular membrane-bounded organelle.

expression levels in defense response. For example, for

Additionally, similar information is shown in Appendix E.

a water deprivation response, two late embryogenesis

These results will help to elucidate how purple wheat

abundant proteins, Lea14-A and Lea14-1, located on the

responds to defense stress. For example, pathogen

membrane, were upregulated. For osmotic stress response,

recognition receptors of plants are predominantly receptor

there were three interleukin-1 receptor-associated kinase 4

serine/threonine (Ser/Thr) kinases that are evolutionarily

located on an intracellular membrane-bounded organelle

related to animal interleukin-1 receptor-associated kinase

(two were upregulated, and one was downregulated), one

(IRAK)/pelle-soluble kinases (Liu et al. 2018). Plants and

glycine-rich RNA-binding protein located on the external

animals mediate early steps of the innate immune response

42

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

Protein name

Gene ID

Dehydration-responsive element-binding protein 2A

TraesCS3A01G099100

Dehydration-responsive element-binding protein 2A

TraesCS3B01G115400

Dehydration-responsive element-binding protein 2C

TraesCS2B01G529200

Dehydration-responsive element-binding protein 2G

TraesCS1B01G256000

Dehydration-responsive element-binding protein 2G

TraesCS1A01G244800

Dehydration-responsive element-binding protein 2G

TraesCS1D01G244500

Dehydration responsive element binding protein

TraesCS2D01G515900

Dehydration responsive element binding protein

TraesCS2A01G514300

CK

Drought

3.5

0

Drought responsive element binding transcription factor TraesCS7B01G028700 Drought responsive element binding transcription factor TraesCS7D01G127600 Drought responsive element binding transcription factor TraesCS7A01G128800 Low temperature-responsive RNA-binding protein

TraesCS4A01G293000

Low temperature-responsive RNA-binding protein

TraesCS4D01G018500

Low temperature-responsive RNA-binding protein

TraesCS4B01G020300

Heat-responsive transcription factor

TraesCS7A01G270100

Heat-responsive transcription factor

TraesCS7B01G168300

Heat-responsive transcription factor HSF85

TraesCS4B01G383300

Heat-responsive transcription factor HSF85

TraesCS4D01G358900

Heat-responsive transcription factor HSF85

TraesCS5A01G548700

Heat-responsive transcription factor HSF85

TraesCS4B01G385800

–3.5

Fig. 7 Heatmap for stress-responsive differentially expressed genes (DEGs) using a filter condition “responsive” from the Nr (the NCBI Non-Redundant Protein Database) annotation.

through pathogen recognition receptors, which commonly associate with or contain members of a monophyletic group of IRAK family that include Drosophila Pelle, human IRAKs, rice XA21, and Arabidopsis FLS2 (Dardick and Ronald 2006). Hence, interleukin-1 receptor-associated kinase 4 might be involved in the innate immune response to pathogens in plants, which indicated that the genes related to biotic stress response in GZ1 might also respond to drought stress. According to the Nr annotation, the keyword “responsive” was used as a filter, and 56 additional DEGs were obtained, of which 20 DEGs encoded 11 dehydration-responsive element-bind proteins/transcription factors (DREBs), three low temperature-responsive RNA-bing proteins, and six heat-responsive transcription factors (Fig. 7). Of them 18 DEGs were upregulated, and only two DREBs were downregulated. These two DREBs were annotated by KEGG as chromodomain-helicase-DNA-binding protein 4, but the other nine DREBs were EREBP (ethyleneresponsive element binding proteins)-like factors. The

other 26 DEGs encoded auxin/ethylene-responsive proteins, chitin-inducible gibberellin-responsive proteins, and calcium-dependent protein kinases (more information is shown in Appendix F). The results revealed that DREBs were differentially expressed under drought stress, and there were more upregulated DREBs than downregulated ones. Moreover, low temperature-responsive RNA-binding proteins and heat-responsive transcription factors were upregulated under drought stress. These results indicated that drought, low temperature, and heat may have a common response mechanism, and auxin- and ethylene-responsive proteins were also induced by drought stress.

3.6. DEGs in secondary biosynthesis Purple wheat can biosynthesize anthocyanin, which distinguishes purple wheat from white wheat. Anthocyanin biosynthesis is a secondary metabolism pathway, and anthocyanin is a compound involved in stress resistance in plants. Therefore, DEGs from secondary metabolic

43

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

pathways were analyzed to determine whether drought stress affected anthocyanin biosynthesis. Among 13 secondary metabolic pathways, the number of upregulated genes was greater than the number of downregulated genes in seven pathways. DEG enrichment was the highest for phenylpropanoid biosynthesis and the second highest for flavonoid biosynthesis (Fig. 8). Anthocyanin biosynthesis involves the phenylpropanoid biosynthesis pathway, the flavonoid biosynthesis pathway, and the anthocyanin biosynthesis pathway. The number of downregulated genes was greater than that of upregulated genes in the first two pathways, though the opposite was observed for the anthocyanin biosynthesis pathway. This may be the reason that the anthocyanin contents in GZ1 grains under drought stress were lower than that in the control. Therefore, we analyzed related genes involved in anthocyanin biosynthesis to obtain more details on the effects of drought stress.

3.7. Expression levels and defense-responsive characteristics of regulatory and structural genes in anthocyanin biosynthesis Anthocyanin biosynthesis is controlled by regulatory and structural genes. In purple wheat, only MYB and bHLH control anthocyanin biosynthesis; their two loci, Pp3 and Pp-D1, are on chromosomes 7D and 2A, respectively. Candidate genes encoding these proteins have been 600

reported in a recent study (Jiang et al. 2018). To explore the mechanisms of anthocyanin biosynthesis under drought stress, we identified regulatory and structural genes based on the KEGG and Nr annotations and results reported in other studies, and created heatmaps according to their log10FPKM values. As shown in the heatmaps, the structural genes TaPAL, Ta4CL, TaCHS, TaF3H, TaF3’H, and TaANS were downregulated, whereas TaDFR, TaOMT, and TaGTs were upregulated. The regulatory genes TaMYB and TaMYC and specifically TaMYB-7D1 and TaMYC-2A1, which have been reported to be purple grain regulatory genes, were downregulated (Fig. 9). From these results, most of the structural and regulatory genes were downregulated, which may have led to the inhibition of anthocyanin biosynthesis. To verify the gene expression levels from the transcriptome, the relative expression levels of the key structural genes TaCHS and TaANS were analyzed by qRTPCR. Comparison of the FPKM values between the control and drought-treated samples showed that TaCHS and TaANS were both down-regulated. The qRT-PCR results also showed the same effect (Appendix B and Fig. 1). The relative expression level trends were consistent with the FPKM value for TaCHS and TaANS, so the FPKM value was used to represent gene expression level. Moreover, the defense response annotations of the structural and regulatory genes were also analyzed, and

548

Up

Number of genes

500 400

Total

Down

348

300 200

200 130

100

71

46

68 44 41 40 52 35 52 22 17 24 20 10 31 12

36 32 22 27 5 4

3

26 23 1623 7

4 4

3 3

Ph

en bi ylpr Fl o o av on syn pa St th no oi ilb es id d en bi oi os is d, yn gi dia th ng ry es er lh is In e o do l b pta le io n al sy oi ka nt d a lo he n id si d bi Is s os of la yn vo th no es id is bi os yn Fl av th es on is py e an rid T bi d in ro o e pa Is sy fla al n oq nt vo ka e, he no ui lo pi no si l id pe s lin b rid e i o al i s n yn e ka th an lo es d id Be is b io nz s ox yn ai th no es id is M bi on os ob yn ac th es ta m is An bi os th yn oc th ya es ni is n bi os yn C th af es fe in is e m e Be ta bo ta lis la in m bi os yn th es is

0

59

Secondary metabolic pathway

Fig. 8 Kyoto encyclopedia of genes and genomes (KEGG) enrichment of differentially expressed genes (DEGs) in secondary biosynthesis.

44

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

CK

Drought

CK

Ta3GT-5A

TaPAL-1B1-T2

TaCHS-2B-T2

Ta3GT-1D

TaPAL-1B2

TaCHS-2D1-T2

TaUGT-7B

TaPAL-1B4-T1

TaCHS-5D1

Ta5,3GT-2A

TaPAL-1B4-T2

TaCHS-5D2

Ta5,3GT-3A

TaPAL-1B5-T1

TaCHS-2B1

Ta5,3GT-5B1

CHS

TaCHS-2D2

Ta5,3GT-6D

TaPAL-1D3

TaCHS-2A2

Ta5,3GT-7A-T1

TaPAL-2B1-T1

TaCHS-2D1-T1

Ta5,3GT-7A-T2

TaPAL-2B2-T1

TaCHS-2B2

Ta5,3GT-3B3

TaPAL-2B2-T2

TaCHS-2B-T1

Ta5,3GT-1A1

TaPAL-6A1

TaCHS-2D3

TaPAL-6B1-T1

TaCHS-2A1-T2

TaPAL-5B2

TaCHI-5B-T2

TaPAL-1A1-T2

Ta5,3GT-3B1 Ta5,3GT-5A

TaPAL-1B5-T2

TaCHI-5B-T1

Ta5,3GT-7B1

TaPAL-1B3

TaF3H-2B

Ta5,3GT-new2

F3H

TaPAL-1A2 TaPAL-1A1-T1 F3’5’H

TaPAL-4A1

TaF3H-2A

Ta5,3GT-1B-1

TaF3H-2D

Ta5,3GT-new1

TaF3'5'H-4D-T1

Ta5,3GT-3B2

TaF3'5'H-4D-T2

Ta5,3GT-1A2

TaPAL-1D4-T2

TaF3'5'H-new

Ta5,3GT-5D2

TaPAL-2A1

TaF3'H-6D1-T2

Ta5,3GT-2B

TaPAL-1D2-T1

TaF3'H-6D2

Ta5,3GT-5B2

TaPAL-2D2

TaF3'H-5A2

Ta5,3GT-1B2

TaPAL-2A2

TaF3'H-2B

Ta5AT-2D-T2

TaPAL-1D1-T1

TaF3'H-2A

Ta5AT-2D-T3

TaF3'H-5B

TaCT-6B

TaPAL-1B6

TaF3'H-5A3

TaHT-5B

TaPAL-5B1

TaF3'H-5A1

TaCT-1A

TaPAL-6B1-T3

TaF3'H-6D1-T1

TaPAL-new

TaF3'H-2D-T1

TaPAL-2B1-T2

TaF3'H-2D-T2

TaCT-4B Acyltransferases TaCT-5B

TaPAL-1D4-T1

TaDFR-2B-T1

TaCT-4D

TaPAL-2D1

TaDFR-2B-T3

TaCT-4A1

TaPAL-6B1-T4

TaDFR-2D1-T1

TaCT-4A2

TaPAL-6B1-T2

TaDFR-3B-T2

Ta5AT-2D-T1

TaPAL-6D1-T2

TaDFR-2D1-T2

TaCT-2D

TaPAL-6D1-T1

TaDFR-2B-T2

Ta5AT-4A

F3'H

TaPAL-2D3

DFR

Ta4CL-2A1

–2

Ta5,3GT-7B2

TaCHI-5D

TaPAL-1D2-T2

TaCT-2A

TaDFR-7B

TaMYB-5A1

Ta4CL-6D2

TaDFR-7A

Ta4CL-6B1-T1

TaDFR-3D

Ta4CL-6B1-T2

TaDFR-2D2

TaMYB-7D1-T1

Ta4CL-6D1

TaDFR-2A

TaMYC-U1

Ta4CL-6A1

TaDFR-7D

TaMYC-U2

TaDFR-3B-T1

TaMYC-4B1

TaANS-6A2

TaMYC-5A2

TaANS-6B

TaMYC-4B2

0

2

ANS

TaANS-6A1 TaANS-6D

TaOMT-5A-T2 MethylTaOMT-5A-T1 transferases TaOMT-5D

Drought

Glycosyl- Ta5,3GT-5D1 transferases Ta5,3GT-3D

TaCHI-5A

CHI

TaPAL-1A3

4CL

CK

TaCHS-2A1-T1

TaPAL-1D1-T2

PAL

Drought

TaPAL-1B1-T1

MYB

MYC

TaMYB-7D1-T2 TaMYB-4B1

TaMYC-U3-T1 TaMYC-U3-T2 TaMYC-U3-T3 TaMYC-2D1 TaMYC-4B3 TaMYC-2B1 TaMYC-2A1

Fig. 9 Heatmap for the regulatory and structural genes in anthocyanin biosynthesis. Ta, Triticum aestivum L. PAL, phenylalanine ammonia-lyase; 4CL, 4-coumarate-CoA ligase; CHS, chalcone synthase; CHI, chalcone isomerase; F3H, flavanone 3-hydroxylase; F3’5’H, flavonoid 3’5’-hydroxylase; F3’H, flavonoid-3’-hydroxylase; DFR, dihydroflavonol 4-reductase; ANS, anthocyanidin synthase; OMT, caffeoyl-CoA O-methyltransferase; 3GT, anthocyanidin 3-O-glucosyltransferase; UGT, cyanidin 3-O-rutinoside 5-O-glucosyltransferase-like; 5,3GT, anthocyanidin 5,3-O-glucosyltransferase; 5AT, anthocyanin 5-aromatic acyltransferase; CT, agmatine coumaroyl-transferase; HT, shikimate O-hydroxycinnamoyltransferase-like; MYB, myeloblastosis transcription factor; and MYC, basic helix-loop-helix transcription factor. 1D, 5A, 7B, etc., represent the chromosomal location. The same chromosome with different numbers represent the corresponding copies of one gene located on the same chromosome. T1, T2, and T3 represent the alternative splicing transcripts of a corresponding gene or copy. new represents a new gene that could not be mapped to the wheat reference genome.

certain structural genes were found to be involved in stress response (Table 2). TaCHS responded to UV and hormones (auxin transport); TaF3H and TaCHI responded to UV; Ta4CL and TaF3’H responded to light stimulus; TaANS responded to stimuli. These data suggest that anthocyanin

biosynthesis may also be affected by stress, especially light stimulus such as UV. However, most of these structural genes were downregulated, with only TaCHI being slightly upregulated, which indicated that drought stress may reduce their defense response.

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

3.8. Protein-protein interaction sites in anthocyanin biosynthesis and their defense-responsive characteristics under drought stress The above analysis explored the drought-regulated expression responses of structural and regulatory genes involved in anthocyanin biosynthesis. To explore how proteins encoded by these genes interact with other proteins under drought stress, protein-protein interaction analysis was performed. There were seven proteinprotein interaction sites encoded by structural genes under drought conditions, including TaCHS-2D1, TaDFR2D2, TaDFR-7D, TaOMT-5A, Ta5,3GT-1B1, Ta5,3GT-3A, and Ta5,3GT-7B1 (Fig. 10). Among them, TaDFR-2D2 interacted with Ta5,3GT-3A, Ta5,3GT-7B1, and TaOMT-5A through other proteins; TaDFR-7D interacted with Ta5,3GT1B1 through other proteins. TaDFR-2D2 had the largest number of interaction proteins (327), followed by TaDFR-7D (166). The information from these interaction sites and networks showed that there was not only an internal interaction between TaDFR and Ta5,3GT and TaOMT but also an external interaction relationship with other proteins. Moreover, the effects of drought on protein-protein interaction relationships were primarily shown for proteins encoded by late structural genes, as only TaCHS-2D1 belonged to proteins encoded by the early structural genes, which had two interaction proteins among the seven sites. To investigate whether these interacting proteins respond to stress, their annotations were analyzed. Of them, 90 proteins were annotated as being involved in defense

45

response reactions (Appendix G). These proteins respond to hormones, other organisms, metal ions, oxidative stress, toxic substances, and abiotic stimuli (such as cold and radiation). The hormone responses included auxin transport, auxin metabolic processes, ethylene metabolic processes, and brassinosteroid metabolic processes. Moreover, proteins that were predicted to interact were involved in many metabolic processes and some developmental processes. The metabolic processes included the metabolism of phytosteroids, phosphatecontaining compounds, disaccharides, DNA, asparagine, fatty acids, and so on. The developmental processes included reproduction, multicellular organism, cell growth, and gametophyte development. The data showed that the interacting proteins were involved in primary metabolism, genetic regulation, and growth and development, which indicated that TaCHS-2D1, TaDFR-2D2, TaDFR-7D, TaOMT-5A, Ta5,3GT-1B1, Ta5,3GT-3A, and Ta5,3GT7B1 connected anthocyanin biosynthesis with other pathways under drought stress. The effects of drought stress on the other pathways may be transmitted to anthocyanin biosynthesis through these structural genes, which subsequently affects anthocyanin biosynthesis and accumulation.

3.9. Defense-responsive network against drought stress in purple wheat According to the annotation of defense-responsive genes and protein-protein interaction networks, a speculative

Ta5,3GT-1B-1

TaOMT-5A

Ta5,3GT-7B1

Ta5,3GT-3A

TaDFR-7D

TaDFR-2D2

TaCHS-2D1

Fig. 10 Protein-protein interaction sites in anthocyanin biosynthesis. The dots represent proteins that are predicted to have protein-protein interaction relationships. The curved lines represent the direction of the protein-protein interaction. The arrows indicate the structural enzymes in anthocyanin biosynthesis.

46

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

defense response mechanism of purple wheat against drought stress is shown in Fig. 11. The defensive responserelated signal transduction transmits drought stress signals to intracellular actors, stimulating biotic and abiotic stressresponse networks, affecting primary and secondary metabolism, and producing protective metabolites against drought stress. The defense mechanism of purple wheat against drought is a very complicated regulatory network, and this speculative model requires many more subsequent verification experiments to verify it.

obtained, and 379 DEGs including DREBs were related to defense responses (Figs. 2, 5, and 7). These defense responses included responses to biotic stress, abiotic stress, hormones, etc. Twenty-nine of these DEGs with P-values less than 1E-200 were involved in responses to reactive oxygen, osmotic stress, water deprivation, plant-type hypersensitive reactions, light intensity, hormones, metal ions, bacteria, and fungi (Fig. 6). Most of the structural and regulatory genes in anthocyanin biosynthesis were downregulated, with only TaDFR, TaOMT, Ta5,3GT, and TaMYB-4B1 being upregulated under drought stress (Fig. 9). TaCHS, TaF3H, TaCHI, Ta4CL, and TaF3’H were involved in responses to UV, hormones, and stimulus (Table 2). TaCHS-2D1, TaDFR-2D2, TaDFR-7D, TaOMT-5A, Ta5,3GT1B1, Ta5,3GT-3A, and Ta5,3GT-7B1 connected anthocyanin biosynthesis with other pathways through protein-protein interaction networks, and their interacting proteins were involved in primary metabolism, genetic regulation, growth and development, and defense responses (Fig. 10). A speculative model of the defense response network against drought stress in purple wheat was created (Fig. 11).

4. Discussion Drought is a common environmental stress factor for wheat (Mahajan and Tteja 2005). However, purple wheat GZ1 can resist drought stress. In our study, drought slightly reduced the anthocyanin, protein, and starch contents of GZ1 grains, but significantly reduced grain weight (Fig. 1). Drought also reduced the number of gene transcripts (Table 1) and suppressed the gene expression of many metabolic pathways (Figs. 3 and 4). A total of 27 766 DEGs were

Anthocyanin biosynthesis TaCHS, TaDFR, TaDMT, TA5, 3GT

Defense Drought

Secondary metabolism

Defense

Primary metabolism

Stress-responsive transcription factors

Drought

Defense

Biotic and abiotic stressresponsive proteins

Drought

Nucleus Defense responserelated signal transduction

Cell wall

Cytoplasm

Membrane

Fig. 11 Defense response pattern of drought stress in purple wheat. Dotted arrows indicate the action path; grey arrows indicate drought stress; white arrows indicate the defense response; the spiral in the nucleus represents the genetic material.

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

These results indicated that purple wheat GZ1 initiated biotic and abiotic stress-responsive genes to resist drought stress and that anthocyanin biosynthesis also participated in drought defense response through several structural genes.

4.1. Defense-responsive genes will contribute to candidate genes for crop breeding of stress resistance The defense-responsive genes from the GZ1 may act as candidate genes for crop breeding of stress resistance. Based on the defense response classification, gene IDs, sequences, and functional annotation were obtained, which provides information for further study of their functions and for the design of molecular markers for stress-resistant breeding. For example, there were seven genes related to water deprivation responses (Appendix E). Among them, there were four upregulated genes; one encoded a predicted protein from Hordeum vulgare subsp. Vulgare and three encoded LEA proteins (late embryogenesis abundant protein Lea14-A). The gene expression levels and P-values of three LEAs were used to evaluate their role in defense response to water deprivation. The activity of two LEA genes from Triticum urartu may be higher than that from Aegilops tauschii. Additionally, a previous study reported that the upregulated LEA genes improve drought and salt tolerance (Wang H Q et al. 2017). Of the three downregulated genes, which all encoded cysteine proteinase, two were from Ae. tauschii and one was from Agropyron mongolicum. Another previous study reported that upregulation of cysteine proteinase inhibitor genes enhances salt, drought, oxidative stress, and cold tolerance (Zhang et al. 2008). Therefore, the downregulation of cysteine proteinase genes may enhance GZ1 drought tolerance. The various sources of defenseresponsive genes may permit differential defense-response abilities. Based on this information, the selection of parents and molecular marker design can be accurately considered for the purpose of stress resistance in hybrid breeding (Leng et al. 2017). In addition, transgenic approaches of defense-responsive genes may be used for improving crop performance (Teng et al. 2017). At the same time, the results also showed that the purple wheat GZ1 may have a dual defense response mechanism to drought stress, wherein upregulated genes enhance tolerance and downregulated stress-sensitive genes reduce negative effects. However, the details of this mechanism and the functions of defenseresponsive genes require to further study.

4.2. Drought-induced defense response to multiple stresses Plants may share some defense mechanisms to resist

47

various stresses. Different stresses can have similar injury effects on plants, and plants resist these injuries with similar defense response mechanisms (Cadiz and Kwon 2008). For example, biotic and abiotic stresses can cause an imbalance in ROS, and plants resist excessive ROS by increasing the biosynthesis of antioxidants such as polyphenols and anthocyanin and the expression of antioxidant enzyme genes such as superoxide dismutase and peroxidase (Hanaka et al. 2018). Therefore, if a plant has one type of defense mechanism, it may use that mechanism to resist several stresses. In our study, drought induced not only the differential expression of genes involved in defense responses to oxidative stress, osmotic stress, and water deprivation but also genes involved in responses to such stimuli as cold, fungi, plant-type hypersensitive reactions, light intensity, hormones, metal ions, bacteria, and heat (Fig. 6). Moreover, drought stress induced the upregulation of DREBs, low temperature-responsive RNAbinding proteins, and heat-responsive transcription factors, which indicated that these transcription factors/proteins may be upregulated simultaneously to resist drought stress (Fig. 7). A previous study reported that salt stress induces the differential expression of genes related to cold, heat, drought, salt, wounds, and light (Goyal et al. 2016). Interestingly, drought stress not only caused a defense response to abiotic stress but also promoted a defense response to biotic stress in GZ1 grains, as 41.38% of the DEGs that had P-values less than 1E-200 were related defense responses to fungi, among which 91.67% were upregulated. These genes encoded basic endochitinase B, peroxidase, DNA-directed RNA polymerase II subunit RPB1, interleukin-1 receptor-associated kinase 1 and 4, and the nuclear pore complex protein Nup214 (Appendix E). In previous studies, the basic endochitinase was shown to confer protection against chitin-containing fungal pathogens (Kim et al. 2010), and an increase in the activity of peroxidase contributes to tomato defense against the nonhost pathogen Magnaporthe grisea (Uma and Podole 2015). These results indicate that these proteins may enhance the defense response to fungi in purple wheat. This information is useful for filtering and the functional identification of disease-resistant genes in wheat. Moreover, this may be one of the reasons why GZ1 has strong disease resistance in agricultural production. However, further study is required to determine which diseases these genes resist.

4.3. Drought-regulated anthocyanin biosynthesis through multiple pathways and several structural genes Anthocyanin is a secondary plant metabolite, and its biosynthesis is regulated by primary metabolism, hormonal

48

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

signals, and external environmental factors (Liu et al. 2018). However, the regulatory network of anthocyanin biosynthesis is highly complex and not fully understood. In a previous study, abscisic acid (ABA) was shown to be involved in the upregulation of anthocyanin biosynthesis under drought stress. ABA increases the expression of anthocyanin synthesis-related genes through the regulation of microRNA156 under drought stress (Gonzalez-Villagra et al. 2017). However, exogenous ABA has no effect on foliar anthocyanin accumulation in purple rice, which is different from plant reproductive organs (Ploenlap and Pattanagul 2015). Biotic stress is involved in the upregulation of anthocyanin biosynthesis under drought stress. Anthocyanin biosynthesis in vitro in grapevine (Vitis vinifera) leaves that were infected with grapevine leafrollassociated virus 3 is enhanced by drought stress (Cui et al. 2017). Drought stress also enhances the upregulation of biosynthesis- or signaling-related genes for JA and anthocyanin in whitefly-infested maize plants (Park et al. 2015). Moreover, UGT79B2 and UGT79B3, two Arabidopsis rhamnosyltransferase genes, are strongly induced by cold, salt and drought stresses and their ectopic expression increases anthocyanin accumulation and antioxidant activity. In contrast, overexpression of these genes in a mutant that cannot synthesize anthocyanin did not increase stress tolerance (Li et al. 2017). These data show that the ABA signaling pathway and biotic stress responses are involved in the upregulation of anthocyanin biosynthesis-related genes under drought stress, and UGT responds to drought stress but is not necessary to increase stress tolerance. In this study, we found that upregulation and downregulation of anthocyanin biosynthesis-related genes occurred at the same time. Some genes were upregulated while others were downregulated. TaDFR-2D2, TaDFR-7D, TaOMT-5A, Ta5,3GT-1B1, Ta5,3GT-3A, and Ta5,3GT-7B1 were protein-protein interaction sites that were identified as having a role under drought stress. The interacting proteins were involved in stress responses, primary metabolism, genetic regulation, hormone metabolism and transport, growth and development. Hormone signalingrelated genes were involved in auxin transport, auxin metabolic processes, ethylene metabolic processes, and brassinosteroid metabolic processes. The associated stress responses included responses to metal ions, oxidative stress, toxic substances, abiotic stimuli (such as cold and radiation), etc., but did not include drought stress responses. This indicated that the role of drought stress on anthocyanin biosynthesis in GZ1 grains was indirect. Moreover, ABA signaling pathway-related proteins were not found among the interacting proteins. However, the auxin, ethylene, and brassinosteroid signaling pathways were involved in the regulation of anthocyanin biosynthesis.

There may be several reasons for these results. First, the anthocyanin content of GZ1 grains was partially reduced, which was inconsistent with other reports that showed drought enhances anthocyanin synthesis. Second, different plants may have different defense-responsive mechanisms to drought stress, and polyploid wheat may have more complicated drought defense mechanisms. Third, drought may increase or decrease anthocyanin biosynthesis through different mechanisms. Fourth, protein-protein interaction analysis of anthocyanin biosynthesis under drought stress can identify more interaction networks and candidate genes. Nevertheless, more studies are needed to elucidate the role of drought in anthocyanin biosynthesis. We also found that drought stress increased anthocyanin biosynthesis in purple wheat GZ1 in other grain developmental stages under different treatments (unpublished data). The defense-responsive mechanisms in these stages may be different. A comparison between purple wheat and white wheat will also be considered to explore whether the defense response of purple wheat is stronger than that of white wheat. These problems will be solved gradually in future studies.

5. Conclusion Overall, our results revealed that purple wheat GZ1 initiated biotic and abiotic stress-responsive genes to resist drought stress and provided 379 defense-responsive candidate genes encoding more than 123 proteins. Anthocyanin biosynthesis-related structural genes participated in multiple stress responses, and seven protein-protein interaction sites of structural enzymes were involved in drought stress-related regulation, which indicated that anthocyanin biosynthesis was affected by drought stress. Moreover, a speculative defense response network against drought stress in purple wheat was created. These results will help to elucidate the stress-resistant mechanisms in purple wheat and are beneficial for the breeding of new excellent crop varieties.

Acknowledgements This study was supported by the grants from the National Key R&D Program of China (2017YFD0100901-4 and 2016YFC0502604), the National Natural Science Foundation of China (31660390), the Major Special Project of Science and Technology Program in Guizhou, China (2017-5411-06 and 2017-5788), the Construction Project of State Engineering Technology Institute for Karst Desertification Control, China (2012FU125X13), the Innovation Talents Team Construction of Science and Technology in Guizhou, China (2016-5624), the

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

Major Research Project of Innovation Group in Guizhou, China (2016-023), and the Graduate Innovation Fund of Guizhou University, China (2017025), and the Science and Technology Project in Guizhou, China (2019-4246). The authors thank Qin Jiangkun, Deng Lei, Wang Liaohong, Qian Xiaokang, Wang Taiyu, and Wei Jinshan (Guizhou University) for their assistance in collecting samples. Appendices associated with this paper can be available on http://www.ChinaAgriSci.com/V2/En/appendix.htm

References Abdel-Aal E S M, Hucl P. 1999. A rapid method for quantifying total anthocyanins in blue aleurone and purple pericarp wheat. Cereal Chemistry, 76, 350–354. Buchfink B, Xie C, Huson D H. 2015. Fast and sensitive protein alignment using DIAMOND. Nature Methods, 12, 59–60. Cadiz N M, Kwon T R. 2008. Drought and salinity stress interactions: An overview. Asia Life Sciences, 17, 337–347. Cui Z H, Bi W L, Hao X Y, Li P M, Duan Y, Walker M A, Xu Y, Wang Q C. 2017. Drought stress enhances up-regulation of anthocyanin biosynthesis in Grapevine leafroll-associated virus 3-infected in vitro Grapevine (Vitis vinifera) leaves. Plant Disease, 101, 1606–1615. Dardick C, Ronald P. 2006. Plant and animal pathogen recognition receptors signal through non-RD kinases. PLoS Pathogens, 2, 0014–0028. Fang Y, Du Y, Wang J, Wu A, Qiao S, Xu B, Zhang S, Siddique K H M, Chen Y. 2017. Moderate drought stress affected root growth and grain yield in old, modern and newly released cultivars of winter wheat. Frontiers in Plant Science, 8, 672. Gonzalez-Villagra J, Kurepin L V, Reyes-Diaz M M. 2017. Evaluating the involvement and interaction of abscisic acid and miRNA156 in the induction of anthocyanin biosynthesis in drought-stressed plants. Planta, 246, 299–312. Goyal E, Amit S K, Singh R S, Mahato A K, Chand S, Kanika K. 2016. Transcriptome profiling of the salt-stress response in Triticum aestivum cv. Kharchia Local. Scientific Reports, 6, 27752. Hanaka A, Lechowski L, Mroczek-Zdyrska M, Strubinska J. 2018. Oxidative enzymes activity during abiotic and biotic stresses in Zea mays leaves and roots exposed to Cu, methyl jasmonate, and Trigonotylus caelestialium. Physiology and Molecular Biology of Plants, 24, 1–5. Jaakola L. 2013. New insights into the regulation of anthocyanin biosynthesis in fruits. Trends in Plant Science, 18, 477–483. Jazayeri S M, Munoz L M M, Romero H M. 2015. RNA-Seq: A glance at technologies and methodologies. Acta Biologica Colombiana, 20, 23–35. Jiang W, Liu T, Nan W, Chamila Jeewani D, Niu Y, Li C, Wang Y, Shi X, Wang C, Wang J, Li Y, Gao X, Wang Z. 2018. Two transcription factors TaPpm1 and TaPpb1 co-regulate the anthocyanin biosynthesis in purple pericarp of wheat. Journal of Experimental Botany, 69, 2555–2567

49

Kim D, Langmead B, Salzberg S L. 2015. HISAT: A fast spliced aligner with low memory requirements. Nature Methods, 12, 357–360. Kim K H, Kamal A H M, Shin K H, Choi J S, Park C S, Heo H Y, Woo S H. 2010. Wild relatives of the wheat grain proteome. Journal of Plant Biology, 53, 344–357. Kong L, Zhang Y, Ye Z Q, Liu X Q, Zhao S Q, Wei L, Gao G. 2007. CPC: Assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Research, 35, W345–W349. Leng P F, Thomas L, Xu M L. 2017. Genomics-assisted breeding - A revolutionary strategy for crop improvement. Journal of Integrative Agriculture, 16, 2674–2685. Li B, Dewey C N. 2011. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12, 323. Li P, Li Y J, Zhang F J, Zhang G Z, Jiang X Y, Yu H M, Hou B K. 2017. The Arabidopsis UDP-glycosyltransferases UGT79B2 and UGT79B3, contribute to cold, salt and drought stress tolerance via modulating anthocyanin accumulation. The Plant Journal, 89, 85–103. Li X, Yang X, Hu Y, Yu X, Li Q. 2014. A novel NAC transcription factor from Suaeda liaotungensis K. enhanced transgenic Arabidopsis drought, salt, and cold stress tolerance. Plant Cell Reports, 33, 767–778. Li X, Zhang M, Ren M, Lü X, Ji N, Wang X. 2017. Molecular regulation mechanism of anthocyanin synthesis in purple wheat. Plant Physiology Journal, 53, 521–530. Livak K J, Schmittgen T D. 2001. Analysis of relative gene expression data using real-time quantitative PCR and the 2(–Delta Delta CT) method. Methods, 25, 402–408 Liu Y, Tikunov Y, Schouten R E, Marcelis L F M, Visser R G F, Bovy A. 2018. Anthocyanin biosynthesis and degradation mechanisms in solanaceous vegetables: A review. Frontiers in Chemistry, 6, 52. Ma D, Sun D, Wang C, Li Y, Guo T. 2014. Expression of flavonoid biosynthesis genes and accumulation of flavonoid in wheat leaves in response to drought stress. Plant Physiology and Biochemistry, 80, 60–66. Mahajan S, Tteja N. 2005. Cold, salinity and drought stresses: An overview. Archives of Biochemistry and Biophysics, 444, 139–158. Marguerat S, Bahler J. 2010. RNA-seq: From technology to biology. Cellular and Molecular Life Sciences, 67, 569–579. von Mering C, Jensen L J, Snel B, Hooper S D, Krupp M, Foglierini M, Jouffre N, Huynen A M, Bork P. 2005. STRING: Known and predicted protein-protein associations, integrated and transferred across organisms. Nucleic Acids Research, 33, D433–D437. Park Y S, Bae D W, Ryu C M. 2015. Aboveground whitefly infestation modulates transcriptional levels of anthocyanin biosynthesis and jasmonic acid signaling-related genes and augments the cope with drought stress of maize. PLoS ONE, 10, e0143879. Pertea M, Pertea G M, Antonescu C M, Chang T C, Mendell J T, Salzberg S L. 2015. StringTie enables improved

50

LI Xiao-lan et al. Journal of Integrative Agriculture 2020, 19(1): 33–50

reconstruction of a transcriptome from RNA-Seq reads. Nature Biotechnology, 33, 290–295. Ploenlap P, Pattanagul W. 2015. Effects of exogenous abscisic acid on foliar anthocyanin accumulation and drought tolerance in purple rice. Biologia, 70, 915–921. Sasaki N, Matsuba Y, Abe Y, Okamura M, Momose M, Umemoto N, Nakayama M, Itoh Y, Ozeki Y. 2013. Recent advances in understanding the anthocyanin modification steps in carnation flowers. Scientia Horticulturae, 163, 37–45. Shen S, Park J W, Lu Z X, Lin L, Henry M D, Wu Y N, Zhou Q, Xing Y. 2014. rMATS: Robust and flexible detection of differential alternative splicing from replicate RNA-Seq data. Proceedings of the National Academy of Sciences of the United States of America, 111, E5593–E5601. Shi G Q , Fu J Y, Rong L J, Zhang P Y, Guo C J, Xiao K. 2018. TaMIR1119, a miRNA family member of wheat (Triticum aestivum), is essential in the regulation of plant drought tolerance. Journal of Integrative Agriculture, 17, 5–14. Springob K, Nakajima J, Yamazaki M, Saito K. 2003. Recent advances in the biosynthesis and accumulation of anthocyanins. Natural Product Reports, 20, 288–303. Teng W, He X, Tong Y P. 2017. Transgenic approaches for improving use efficiency of nitrogen, phosphorus and potassium in crops. Journal of Integrative Agriculture, 16, 2657–2673.

Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley R D, Pimentel H, Salzberg S, Rinn L J, Pachter L. 2012. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nature Protocols, 7, 562–578. Uma B, Podile A R. 2015. Apoplastic oxidative defenses during non-host interactions of tomato (Lycopersicon esculentum L.) with Magnaporthe grisea. Acta Physiologiae Plantarum, 37, 26. Wang F, Dong X, Tang X, Tu L, Zhao B, Sui N, Fu D, Zhang X. 2018. Comparative transcriptome analysis revealing the effect of light on anthocyanin biosynthesis in purple grains of wheat. Journal of Agricultural and Food Chemistry, 66, 3465–3476. Wang H Q, Wu Y C, Yang X B, Guo X R, Cao X Y. 2017. SmLEA2, a gene for late embryogenesis abundant protein isolated from Salvia miltiorrhiza, confers tolerance to drought and salt stress in Escherichia coli and S. miltiorrhizo. Protoplasma, 254, 685–696. Wang J, Pan Y, Shen S, Lin L, Xing Y. 2017. rMATSDVR: rMATS discovery of differential variants in RNA. Bioinformatics, 33, 2216–2217. Zhang X, Liu S, Takano T. 2008. Two cysteine proteinase inhibitors from Arabidopsis thaliana, AtCYSa and AtCYSb, increasing the salt, drought, oxidation and cold tolerance. Plant Molecular Biology, 68, 131–143.

Executive Editor-in-Chief ZHANG Xue-yong Managing editor WANG Ning