Stress-induced changes in the expression of antioxidant system genes for rice (Oryza sativa L.) and bread wheat (Triticum aestivum L.)

Sampling and clustering of experiments

The GEO Datasets (Edgar, Domrachev & Lash, 2002; Barrett et al., 2012) were used as a source of data for the expression profiles of antioxidant genes for rice (O. sativa) and bread wheat (T. aestivum) obtained in control and stress (water deficiency and cold) conditions (Barrett et al., 2005). The initial sample consisted of 22 experiments for rice (GSE21651, GSE23211, GSE24048, GSE25176, GSE26280, GSE37940, GSE38023, GSE41647, GSE4438, GSE54466, GSE64576, GSE65024, GSE65025, GSE67373, GSE6901, GSE71680, GSE74465, GSE78972, GSE80811, GSE81462, GSE83378, GSE83912), and seven experiments for wheat (GSE23889, GSE30436, GSE31762, GSE42214, GSE45563, GSE47090, GSE79522), including data from microarrays and NGS platforms (see Table S1).

Each experiment contains genes expression profile corresponding to control and to one or several variants of stress conditions that differ in time and/or intensity of cold and/or water deficiency.

From these experiments, the expression profiles of antioxidant genes for seven classes of AOS enzymes were extracted. NGS data were normalized by the FPKM method. Microarray and NGS data were subsequently normalized to the expression of housekeeping genes: the UBQ10 gene (OS02G06640) was used for rice and the tubulin 4 alpha chain gene (TRIAE_CS42_5BL_TGACv1_405584_AA1330850) was used for bread wheat. Then, for each of the selected antioxidant genes in each experiment, the average expression was calculated for all replicates in control conditions, which was used to normalize the expression values of this gene in both control and stress conditions. The obtained normalized values reflect the amount of deviation of the expression of each gene from its average expression level under control conditions. Thus, for all selected experiments, we obtained a set of replicas consisting of the normalized values of the expression for each gene under study.

For rice, all experimental points were clustered by the Ward algorithm based on the obtained profiles using the Past program (Hammer, Harper & Ryan, 2001). Then, according to the clustering, we selected experiments that were characterized by an explicit stress response by the following rule: if more than 50% of stress replicates fell into a cluster with control replications, then such an experiment was not taken for further analysis, since it contains an implicit response to stress or noisy data. Thus, the final sample for rice included 13 experiments: GSE21651, GSE25176, GSE26280, GSE37940, GSE38023, GSE41647, GSE54466, GSE64576, GSE65024, GSE65025, GSE6901, GSE78972, GSE81462 (Table S1) containing 200 biological replicates, including 70 control and 130 water deficient and cold replicas. The resulting sample also was clustered by the Ward algorithm.

For bread wheat, we did not preselect experiments because of the small sample size. Therefore, we analyzed changing of expression levels of AOS genes in individual experiments and performed statistical analysis between control and stress replicas to investigate statistically significant response AOS genes.

Statistical analysis

Clustering and principal component analysis (PCA) was obtained using the Past program (Hammer, Harper & Ryan, 2001). Statistical data processing on clusters for rice and on types of effects inside individual experiments for wheat was carried out using a t-test with the Benjamini–Yekutieli correction. Samples of expression characteristics in stress clusters for rice and replicas for wheat were compared against control replicas and E-value were calculated. The genes of up- and down-regulation in stress were marked by these values.

Genetic material and growing conditions

For an experimental study, bread wheat T. aestivum L. cv. Saratovskaya 29 (S29) and cv. Yanetzkis Probat (YP) were used. The cv. S29 possesses an increased tolerance to stress conditions and the cv. YP, on the contrary, is less tolerant (Doroshkov, Pshenichnikova & Afonnikov, 2011; Doroshkov et al., 2016). Therefore, this pair of genotypes was chosen to study stress-induced changes in the expression of the AOS genes in plant leaves.

Germinated seeds of plants for both varieties were planted in groups of at least six plants of each variety in plastic containers (PET, size 56 × 39 × 28 cm) with drainage and filled by expanded clay balls. Plants were grown up to the four-leaf stage in a greenhouse (Laboratory of Artificial Plant Growth) at the Institute of Cytology and Genetics SB RAS under a 18/6 h photoperiod (artificial supplementary lighting by Sylvania SHP-TS 600w, illumination of plants was carried out as 15–20 × 10³ lx); the temperature condition was maintained at 14–16 °C at night and 20–23 °C during the day. Watering with Knop solution (Hershey, 1994) was carried out daily, until complete saturation of the substrate. Knop solution with additional microelements was used for watering once a week.

Stress experiment design and sampling

For stress-modeling experiment for each experimental group we selected containers in which at least four to five plants of each variety grew at the same stage of development. Selected containers with wheat plants were simultaneously either subjected to cold stress (stress groups) or left without stress exposure in described previously growing conditions (control group). Cold stress was simulated by placing containers with plants of stress groups in a climatic chamber with constant temperature +4 °C and same as in greenhouse day/night lighting conditions for six (short-term) and 24 h (long-term cold stress conditions).

Each plant material sample was obtained by pooling of the third leaf middle part from three plants of same cv. grown in one container. The total weight of each sample was 100–200 mg. In this way three samples of plant material per each experimental group of plants was taken. After collection, the samples were immediately frozen in liquid nitrogen and stored at −70 °C until RNA isolation. Sampling was carried out simultaneously for the control and both stress groups of plants.

RNA isolation and RT-qPCR analysis

Total RNA was extracted from the collected samples using Plant RNA MiniPrep Kit (Zymo Research, Irvine, CA, USA) with additional treatment with RNase-Free DNase Set (QIAGEN, Hilden, Germany) according to the manufacturers instructions and protocols. Evaluation of the quantity and purity of RNA was performed on a Nanodrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA), integrity was assessed using agarose gel electrophoresis. For cDNA synthesis, an OT-M-MuLV-RH Reagent Kit (Biolabmix, Novosibirsk, Russia) was used, the reaction was performed according to the manufacturer’s standard protocol with oligo(dT)₁₆ primers in 20 μl per reaction volume.

For the experimental analysis, nine orthologs from key classes of APX, SOD and CAT enzymes were selected. The nomenclature of selected genes corresponds to our previous paper (Doroshkov & Bobrovskikh, 2018).

The primers were designed complementary to consensus sequence for all three homeologues of the CAT, APX and SOD genes. Design of the primers was made with the online service PrimerQuest Tool (Integrated DNA Technologies, Coralville, IA, USA) based on the Primer3 algorithm (Untergasser et al., 2012) with default conditions (Ta_opt = 60 °C, option “qPCR 2 primers with intercalating dyes”), and we also used the database of nucleotide sequences EnsemblPlants release 37 (Kersey et al., 2015). The amplicon location was within 2,000 bp from the 3′-end and the last intron-exon boundary was chosen in cases where this was possible. Specificity of the primers was checked in silico using the database of nucleotide sequences (NCBI Resource Coordinators, 2019) and algorithms NCBI BLAST, Primer-BLAST (Altschul et al., 1990; Ye et al., 2012). The gene encoding 18S rRNA was used as a reference gene. The primer sequence for the 18S rRNA gene was used from Beales et al. (2007). The specificity test was carried out in a similar way as for the other primers. The sequences of primers used in this work are shown in Table 1.

Oligonucleotides were synthesized in the Biosset Ltd. (Novosibirsk, Russia). Verification of the amplicon specificity was carried out by imaging on agarose gel electrophoresis and building a melting curve using a Lightcycler 96 (Roche Diagnostics, Indianapolis, IN, USA). Optimization of the amplification conditions was not required.

Real-time PCR was performed using a standard three-step protocol on a Lightcycler 96 Instrument (Roche Diagnostics, Indianapolis, IN, USA) with low-profile 96-well plates (#PCR-96-LP-FLT-C; Corning Axygen Scientific, Union City, CA, USA), Ultra-Clear optically transparent films (#UC-500; Corning Axygen Scientific, Union City, CA, USA), and BioMaster HS-qPCR SYBR Blue reagent kit (2×) (Biolabmix, Novosibirsk, Russia). Pre-denaturation was 7 min at 95 °C; denaturation was 18 s at 95 °C; primer annealing was 20 s at 60 °C; elongation was 15 s at 72 °C with detection at the end of the cycle; 45 cycles. The total volume of one qPCR reaction was 20 μl. One qPCR reaction contains two μl of diluted with water (1:5) cDNA and forward/reverse primers to final concentration 0.25 μM each. The reaction with the addition of ddH₂O instead of cDNA was used as a negative control (NTC). Three biological and three technical replicates were used in each qPCR experiment.

Relative expression calculation

The genes of interest (GOI) relative expression was determined by ΔΔC_q method. qPCR data represent the differential expression of AOS genes under cold stress conditions. For this the raw data from LightCycler 96 Instrument was analysed by LinRegPCR software (Ruijter et al., 2009) to determine C_q and efficiency of reactions and then relative ratio of GOI expression in stress conditions was calculated according proposed by Pfaffl (2004) multiple samples efficiency corrected mathematical model (see Eq. (1)) in LibreOffice Calc (The Document Foundation, Berlin, Germany).

(1) $$\mathrm{Ratio}=\frac{{(E_{\mathrm{target}})}^{\mathrm{\Delta }{C}_{q\mathrm{target}}\left({\mathrm{MEAN}}_{\mathrm{control}}-{\mathrm{MEAN}}_{\mathrm{treat}}\right)}}{{(E_{\mathrm{reference}})}^{\mathrm{\Delta }{C}_{q\mathrm{reference}}\left({\mathrm{MEAN}}_{\mathrm{control}}-{\mathrm{MEAN}}_{\mathrm{treat}}\right)}}$$

Further data analysis was performed with R (R Development Core Team, 2018), R-studio IDE (RStudio Team, 2015) and xlsx R package (Dragulescu & Arendt, 2018). The plots representing expression ration were generated with RColorBrewer (Neuwirth, 2014), ggplot2 (Wickham, 2016) and ggpubr (Kassambara, 2018) R packages.

Cold and water deficiency response experiments for rice

For rice, as a result of selection and preprocessing of experiments outputs (see detailed description in “Materials and Methods”) and their further cluster analysis, the replicates were separated into a control cluster, five cold response clusters (CR), and six water deficiency response (WDR) clusters (Fig. 1A; Table S2). The PCA showed that WDR more often combines experiments into clusters, but at the same time, CR clusters are more compact and even all together occupy an area similar in size to the control cluster (Figs. 1B and 1C). The smaller scatter of points inside the clusters and the positions of the clusters relative to each other indicates a more systemic response of the AOS to cold stress. Also, most often, replications from different experiments were combined into one cluster, which shows the similarity of their response to stress and indicates the characteristic behavior pattern of antioxidant genes.

The contribution of individual AOS genes in CR and WDR

For rice, we investigated 35 genes encoding AOS components. In CR clusters, we observed 31 AOS genes, which characterized by explicit up- or down-regulation. Also, we obtained eight genes (marked below by asterisk), which are significantly downregulated in all cold-response clusters.

For CR clusters, the following genes are upregulated:

• APX: OS04G14680

• GPX: OS04G46960

• GR: OS10G28000

• MDHAR: OS08G44340

For CR clusters, the following genes are downregulated:

• APX: OS04G35520, OS02G34810, OS12G07820*, OS12G07830, OS07G49400*, OS08G43560

• CAT (all genes): OS06G51150, OS03G03910, OS02G02400

• DHAR (all genes): OS05G02530*, OS06G12630*

• GPX: OS06G08670, OS03G24380, OS11G18170*

• GR: OS02G56850*, OS03G06740*

• MDHAR: OS09G39380*, OS02G47800, OS02G47790, OS08G05570

• SOD (all genes): OS04G48410, OS06G05110, OS06G02500, OS03G22810, OS07G46990, OS08G44770, OS05G25850.

In WDR clusters, we observed ten AOS genes, which characterized by clear up- or down-regulation. Also, we obtained three genes, two of them were downregulated (marked below by asterisk), and gene OS06G12630 (DHAR) has a mixed regulation.

For WDR clusters, the following genes are upregulated:

• GR: OS10G28000 (also observed for CR)

• MDHAR: OS08G44340 (for five from six WDR clusters, which is also observed for CR)

For WDR clusters, the following genes are downregulated:

• APX: OS08G41090, OS07G49400 (also observed for CR), OS08G43560 (also observed for CR), OS12G07830 (also observed for CR),

• CAT: OS03G03910, OS02G02400* (for both also observed for CR),

• GPX: OS06G08670, OS11G18170* (for both also observed for CR),

• MDHAR: OS02G47800, OS02G47790,

• SOD: OS03G11960, OS03G22810 (for all also observed for CR).

We compared upregulated genes during stress response with their characteristics studied in Doroshkov & Bobrovskikh (2018). Gene OS04G14680; APX H (highly conservative with low expression level) significantly upregulated in four clusters of stress-response (CR1, CR5, WDR1, WDR6). Gene OS06G51150; CAT A (highly conservative, expressed highly in roots) significantly upregulated in three WDR clusters (WDR1, WDR2, WDR5). Gene OS04G46960; GPX E (highly conservative with high expression level) significantly upregulated in CR clusters (CR1, CR2, CR4, CR5). SOD genes OS04G48410; SOD G (low conservative); OS07G46990; SOD E (highly conservative with high expression level); OS05G25850, SOD C (low conservative with low expression level) are significantly upregulated in WDR clusters (WDR1, WDR5, WDR6; WDR1, WDR2, WDR5, WDR6; WDR1, WDR2, WDR5, WDR6, accordingly). Also, some genes that were not characterized in the previous study were upregulated: OS10G28000; GR A (eight clusters); OS05G02530, DHAR B (WDR5, WDR6); OS08G44340, MDHAR A (seven clusters). So, four of six genes which founded upregulated are highly conservative, which can indicate the importance of this characteristic for pre-selection of important AOS genes

Cold and water deficiency response experiments for bread wheat

For bread wheat (T. aestivum) due to the small sample size (seven experiments), we did not preselect experiments. Therefore, we analyzed separated experiments with specific stress conditions separately (Table S3).

By analyzing the microarray experiments, we obtained data only for some AOS genes because platform GPL3802 had only some samples of AOS genes. Also, microarray probe measure total expression of all homeological genes, which imposes additional restrictions on measurement accuracy. Platform included probes for APX B, APX G, GPX D, SOD C, SOD F, CAT A, CAT B, GR A, GR B, DHAR B. For the other hand, analysis of NGS experiment GSE79522 provide data about each copy of AOS genes in bread wheat.

Stable and most interesting data about regulation of AOS genes during stress was obtained from experiments GSE79522 and GSE23889. In water-deficient experiments GSE79522 we revealed a significant upregulation of CAT A, APX G, SOD D, CAT C, MDHAR A significantly and downregulation of APX F and DHAR A. Analysis of cold stress treatment of experiment GSE23889 showed significantly downregulation of following genes: CAT A, CAT B, GR A, GR B, DHAR B between control and stress replicas. By analyzing the expression of antioxidant genes depending on the time of acclimatization on GSE23389 (see link) (2, 14, 21, 35, 42, 56, 70 days) was significantly revealed upregulation for GPX-E (in most cases for spring Manitou), SOD-F (some time points for every genotype), GR-B for winter Manitou (2 days, 21 days) and spring Manitou (2 days), DHAR-B (in most cases of winter Manitou and spring Manitou and downregulation of GR A for spring Norstar (2, 35, 56, 70 days).

So, the following analysis on bread wheat data gives an incomplete picture about regulation of AOS genes in stress conditions. Therefore, we decided to study AOS genes in our experiment, which include different groups of genes of AOS of bread wheat.