Genome-wide transcriptome analysis of roots in two rice varieties in response to alternate wetting and drying irrigation

To Song,Debtosh Ds,Feng Yng,Moxin Chen,Yun Tin,Cholin Cheng,Cho Sun,Weifeng Xu, Jinhu Zhng,e,f,

aShenzhen Research Institute,The Chinese University of Hong Kong,Shenzhen 518057,Guangdong,China

bSouthern Regional Collaborative Innovation Center for Grain and Oil Crops in China, Hunan Agricultural University, Changsha 410128,Hunan,China

cCollege of Horticulture,Nanjing Agricultural University,Nanjing 210095,Jiangsu,China

dCenter for Plant Water-use and Nutrition Regulation and College of Life Sciences, Joint International Research Laboratory of Water and Nutrient in Crop,Fujian Agriculture and Forestry University, Fuzhou 350002,Fujian,China

eDepartment of Biology, Hong Kong Baptist University,Kowloon,Hong Kong,China

fState Key Laboratory of Agrobiotechnology,The Chinese University of Hong Kong,Shatin,Hong Kong,China

ABSTRACT Alternate wetting and drying (AWD) irrigation has been widely used as an efficient rice production method to obtain better yield without continuous flooding (CF) of the paddy field. However, how this practice affects gene expression to regulate rice physiology and morphology is largely unknown. In this study, we used two rice varieties, Nipponbare, a lowland rice cultivar, and Gaoshan 1, an upland cultivar, and found that root dry weight(RDW) and root oxidation activity (ROA) in both cultivars substantially increased in response to AWD.We then analyzed the differences in transcriptome profiles of their roots irrigated in AWD vs. CF conditions. AWD responsive genes are mainly involved in lignin biosynthetic pathway and phytohormone signal transduction pathway and belong mainly to bHLH,bZIP,NAC,WRKY,and HSF transcription factor families.We discussed how these differentially expressed genes may contribute to the morphological adaptations observed in roots exposed to AWD. This analysis also provides useful information to explain the similarities and differences in adaptation to AWD irrigation between the two rice ecotypes.

1. Introduction

As one of the earliest domesticated food crops, rice (Oryza sativa L.) is an important staple food for half of the world population [1,2]. China shares about 19% of the global rice planting area and contributes to 32% of the global rice production (FAO, http://www.fao.org/faostat/zh/#data), and the amount of rice irrigation water accounts for 65% of annual agricultural water consumption in China [3]. Unfortunately,large parts of China face either physical or economic water scarcity [4]. Moreover, agricultural irrigation water will gradually and rapidly deplete due to fierce demand for water resources from urban and industrial sectors, and with increasing global commercialization it seems that industry will receive priority over irrigation [5]. In China, over 95% of rice is grown under traditional continuous flooding (CF)irrigation which expends a large amount of labor, time and energy due to the increased pumping of water in flooded fields [6]. Alternate wetting and drying (AWD) irrigation has been developed as a novel water-saving technique and has been adopted in many countries such as China, Bangladesh,India and Vietnam [7,8]. By reducing the required number of irrigation events,AWD irrigation can reduce water consumption by up to 30%-35% in comparison to CF irrigation [9,10].Among rice varieties, while lowland rice is grown in deep water irrigation and plants possess shallow-thin roots,upland rice cultivation involves dry fields and plants acquire deeper thicker roots [11,12]. Upland rice varieties have evolved to adapt to the drought-prone regions with increased drought tolerance due to the long-term natural and human selection[13]. Thus, global water shortage has forced farmers to go on with upland rice cultivation.

The underlying genetics behind this evolution of upland rice from lowland varieties is currently unexplored,except for some transcriptome and few quantitative trait loci (QTL)mapping studies [14-16]. In a cross between shallow-rooted lowland and deep-rooted upland variety, QTL were searched for root growth and were mapped to (DRO1, QRO1, QRO2,QBRT3.1-3.2, QBRT8.2, and QRT9.1-9.2) [12,17-19]. To analyze adaptation process of upland rice to aerobic conditions, a study based on 5779 single nucleotide polymorphisms (SNPs)was conducted and found that upland rice ecotypes have most robust roots (long and thick) and very high number of robust root alleles than other rice ecotypes[20].SNP genotyping analyses between upland and lowland landraces found that two potentially drought-resistance genes (ARAG1 and OsGL1-8) underlying root architecture-associated drought avoidance may have undergone directional selection in upland rice [21]. Comparative transcriptome analysis similar to our approach was used to state the molecular mechanism of stress adaptation in upland rice, and ecotype differentiation genes were enriched in ROS alleviating gene families such as peroxidases, glutathione related classes and in phytohormone metabolism and signaling factors and transcription factors[22-24].

In flooded soils, roots develop in the shallow soil layer,which favors nutrient uptake from the floodwaters.In aerobic soils, by contrast, root growth is more dispersed [25]. AWD irrigation can increase root biomass,root length density,total root absorption area, active root absorption area and zeatin(Z) + zeatin riboside (ZR) content in roots when compared to those under CF irrigation which may contribute to higher grain yield and water productivity[26].Rice shows greater root activity under AWD irrigation when compared to CF irrigation[27]. More developed root systems of upland rice helped the plants maintain a higher water status than that maintained by lowland rice when the plants were subjected to soil drying[28]. Phytohormones may play important regulatory role in crop growth by acting as signaling molecules to regulate crop physiological function and metabolism [29]. Accumulation of Abscisic acid (ABA) regulates auxin transport in the root tip,which enhances proton secretion (protons cause cell extension by cell wall loosening and then cause root growth) for maintaining root growth under moderate water stress [30].Many transcription factors (TFs) are also reported to regulate root growth [31,32]. However, regulation of these phytohormones and TFs during root growth in lowland and upland rice under AWD is yet to be understood fully.

In recent times, utilization of transcriptome approaches such as microarrays and next-generation sequencing have helped researchers to fish out genome-wide gene expression changes. RNA-Seq uses deep-sequencing technology to provide a far more precise measurement of the levels of transcripts and their isoforms than other methods such as microarrays [33]. RNA-seq allows gene discovery and global gene expression profiling, for example, to identify key signaling components of pathogen-resistance pathways [34].With the rapidly decreasing cost,RNA-Seq has been applied to numerous studies in response to abiotic stress treatments[35-38]. In our previous study, both upland and lowland rice showed commonality in root phenotype when irrigated under AWD [28]. In this study, we used RNA-Seq to explore molecular reasons for this common phenotypic adaptation to AWD. Based on the comparative analysis of the gene expression under different irrigation regimes,the information obtained hereby revealed useful clues to explain the similarities and differences in AWD irrigation transcriptome response between upland rice and lowland rice.

2. Materials and methods

2.1.Plant materials and growth conditions

Pot experiments were conducted in the greenhouse of Shenzhen Research Institute, the Chinese University of Hong Kong, Shenzhen, Guangdong, China (22°32′N, 113°56′E), during the early rice-growing season(Mar-June),as mentioned in Wang et al. [39]. Each pot (30 cm in height and 34 cm in diameter) was filled with 15 kg air-dried fine-grained soil collected from paddy fields.The experimental soil was sandy loam and the main properties were as follows: pH 4.91 (1:5,soil:water); total N, 1.73 g kg−1; total P, 0.72 g kg−1; total K,29.3 g kg−1; organic C, 22.5 g kg−1; alkali-hydrolysable N,216 mg kg−1; Olsen-P, 36.7 mg kg−1; and exchangeable K,115.0 mg kg−1.In the current study,two japonica rice varieties,Nipponbare(a lowland rice variety)and Gaoshan 1(an upland rice variety),were used in the experiment.Seeds were sown in the nursery beds on Mar 31, and 4-week-old seedlings were transplanted to pots on April 28, with each pot having six seedlings.Three conventional fertilizers were used:urea(46%N) as nitrogen fertilizer, superphosphate (12% P2O5) as phosphorus fertilizer, and potassium chloride (60% K2O) as potassium fertilizer. Before transplanting, basal fertilizers were applied as 1.0 g N, 1.2 g P2O5, and 0.9 g K2O per pot. N topdressing was applied at two stages, the panicle initiation stage(0.8 g N per pot)and thereafter,15 days later(0.8 g N per pot).

The two irrigation regime treatments were CF and AWD.The irrigation regimes were applied from 14 days after transplanting (DAT), when the seedlings had recovered from transplanting injury, up to 7 days before harvest. While CF treatment consisted of allowing a consistent level of water in the field to a depth of 2-5 cm,AWD involved water application only when the soil water potential in 15-20 cm soil layer reached −15 kPa. Water potential was measured using tensiometers (two tensiometers in each AWD irrigated pot)at midday each day.According to the measurement,when the soil water potential reached −15 kPa under the AWD irrigation, water was applied to a depth of 2-3 cm in the corresponding AWD pots. The quantity of water applied was quantified based on the number of cups used to add water.In total,15 and 9 irrigation events were applied before sampling date in CF and AWD irrigation,respectively.The heading date of Nipponbare under AWD and CF was between June 18-26.The heading date of Gaoshan 1 under AWD and CF was between June 20-28.

2.2. Root dry weight (RDW) and root oxidation activity (ROA)

Rice plants were sampled destructively from the 12 dedicated pots at heading stage (June 23), i.e. 3 replicates × 2 irrigation regimes × 2 varieties. Each pot contained three root samples which were collected to analyze RDW,ROA and transcriptome respectively. For RDW, all of the above-ground tissues were separated by cutting the rice plants at soil surface and roots were collected by gently removing the surrounding saturated soil and washing over a 2000 μm mesh size sieve. Samples were dried in an oven at 70°C until constant weight,and then RDW was measured.

To measure the ROA,rice roots were sampled destructively from the 12 pots as mentioned above.The roots were rinsed to remove soil and then cut from their nodal bases.The ROA was determined according to the method described in Ramasamy et al. [40]. The ROA was measured by determining the oxidation of alphanaphthylamine (α-NA) and was expressed as μg α-NA per gram dry weight (DW) per hour (μg α-NA g−1DW h−1).

2.3. RNA sequencing and data analysis

2.3.1. RNA extraction and quality check

CF and AWD RNA samples were collected at the heading stage 2 h after AWD pot rewatering (to ensure that the soil water potential reaches a value of 0). Rice roots were sampled destructively from the 12 pots as mentioned above.The roots were rinsed quickly, cut from their nodal bases and rapidly frozen in liquid nitrogen.Total RNA was extracted using EZNA plant RNA Kit (Omega Bio-Tek, GA, USA) in accordance with the manufacturer’s instructions. RNA purity was checked using the kaiaoK5500 Spectrophotometer (Kaiao, Beijing,China). RNA integrity and concentration were assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system(Agilent Technologies, CA, USA). RNA concentration for library preparation was measured using Qubit RNA Assay Kit in Qubit 3.0 to preliminary quantify and then dilute to 1 ng μL−1.

2.3.2. Library preparation for RNA sequencing

For library preparation, 2 μg of total RNA were input. NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) was used to generate sequencing libraries following the manufacturer’s recommendations. In short, mRNA was purified from the input total RNA using poly-T oligo-attached magnetic beads.This was followed by RNA fragmentation with addition of divalent cations under increased temperature. Then, first strand cDNA synthesis was performed using random hexamer primers and remaining RNA was degraded using RNase H. Subsequently, a second strand cDNA synthesis was performed. Resulting fragments were purified with QiaQuick PCR kit and followed by terminal repair, A-tailing and adapter addition. Finally, PCR was performed to finish library preparation.

2.3.3. Library examination, clustering and sequencing

Library insert size was quantified using StepOnePlus Real-Time PCR System (Library valid concentration ＞ 10 nmol L−1).Sample clustering was performed on a cBot cluster generation system using HiSeq PE Cluster Kit v4-cBot-HS (Illumina).Subsequently, libraries were sequenced on an Illumina platform to obtain 150 bp paired-end reads.

2.3.4. Transcriptome analysis

Reads were processed through quality check using FastQC protocol and high-quality reads (referred to as “clean reads”)were obtained for downstream mapping to reference genome of rice (ftp://ftp.ensemblgenomes.org/ pub/plants/release-24/fasta/oryza_sativa/) in HISAT2 v2.0.5 [41]. Following this, raw read count for each gene in each sample was obtained with HTSeq v0.6.0, and FPKM (Fragments Per Kilobase Million mapped reads) was calculated to estimate the expression level of genes in each sample. Finally, DESeq2 v1.6.3 was used for differential gene expression analysis between samples.Briefly, it estimates the gene expression level by linear regression calculating the fold changes for sample comparisons and then calculates the p-value with Wald test. Finally,the p-value was corrected by the BH method to give q-value.Genes with q ≤ 0.05 and |log2fold change| ≥ 1 were identified as differentially expressed genes (DEGs). Functional enrichment analyses were performed using hierarchical clustering, functional classification and pathway analyses [42].

2.4. Quantitative real-time PCR validation

Validation of RNA-seq results was performed for 18 randomly chosen genes using quantitative RT-PCR. The sequences of the gene primers are provided in Table S1.Total RNA was extracted as mentioned in the RNA-Seq experiment and first-strand cDNA was generated using a TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix kit (TransGen Biotech). Quantitative RT-PCR was carried out in the ABI Step One Plus Real-Time PCR System(Applied Biosystems, USA) utilizing the SYBR Premix Ex Taq RT-PCR kit (Takara Bio, Japan). The rice gene Actin1(Os03g0718100) was used as a housekeeping gene for expression normalization.

2.5. Statistical analysis

Analyses were conducted using IBM SPSS Statistics ver. 18(SPSS Inc.,IL,USA)with Duncan’s Multiple Range test(DMRT)post-hoc test, and the significance level was set to P ＜0.05.Data are means ± SE from three independent biological replicates.

3. Results

3.1. Both rice cultivars display higher RDW and ROA under AWD

For the two rice cultivars, RDW and ROA were higher under AWD irrigation when compared to CF irrigation (Fig. 1).Compared to the lowland rice cultivar Nipponbare,the upland rice cultivar Gaoshan 1 had higher RDW under both AWD and CF irrigation regimes (Fig. 1-A), and higher ROA under AWD irrigation(Fig.1-B).

3.2. Cultivar type was the major factor contributing to transcriptome variation followed by irrigation scheme

A principal component analysis(PCA)was conducted with the clean data of 12 samples to assess transcriptional variation among the cultivar samples under different irrigation regimes and assess within-sample replicate variation (Fig. 2-A). The replicates for each treatment-cultivar combination clustered closely with each other,suggesting excellent replication in the RNASeq experiment. Along PC1, major variation was seen at the level of the cultivars.PC2 showed a large variation among irrigation regimes but only in the lowland rice, suggesting that irrigation regimes strongly affected the transcriptome of lowland rice when compared to the upland rice.

3.3. Nipponbare displayed strong transcriptome modulation between different irrigation schemes

Differences in gene expression were examined to decipher the genes that may participate in root growth regulation under the two contrasting irrigation regimes. When comparing AWD samples with CF ones, there were 6597 and 526 DEGs in the roots of Nipponbare and Gaoshan 1, respectively (Fig. 2-B). For Nipponbare, 2406 genes were up-regulated and 4191 genes were down-regulated under AWD irrigation compared to CF irrigation. For Gaoshan 1, comparison of AWD irrigation with CF irrigation revealed 354 genes that were up-regulated and 172 genes that were down-regulated. When comparing between Nipponbare and Gaoshan 1, we found DEGs of Nipponbare were 12.5-fold to those of Gaoshan 1. More DEGs were downregulated under AWD in Nipponbare. Contrastingly, in Gaoshan 1 more DEGs were up-regulated under AWD.Venn diagrams were created to analyze the gene expression similarities and differences between Nipponbare and Gaoshan 1. There were 129 common up-regulated DEGs and 25 common downregulated DEGs between the two cultivars(Fig. 2-C; Table S2, Table S3).

Fig.1- Root dry weight (A)and root oxidation activity(B) of a lowland(Nipponbare) and an upland(Gaoshan 1)variety under two contrasting irrigation treatments.CF,continuous flooding;AWD,alternate wetting and drying irrigation.Letters above the error bars indicate significant difference at P ＜0.05.

3.4. Functional classification of the DEGs by GO and KEGG pathway analysis

To further characterize genes affected under AWD, the DEGs were subjected to Gene Ontology (GO) enrichment and KEGG pathway analyses. DEGs of Nipponbare (Nip-CF vs. Nip-AWD)and Gaoshan 1 (GS-CF vs. GS-AWD) were classified in three main GO categories: “biological processes”, “cellular component” and “molecular function” (Figs. 3, 4). In the biological process category, the GO terms most enriched were in cellular process, metabolic process, response to stimulus and biological regulation for both cultivars. In the cellular component category, cell part, organelle, membrane part, membrane,organelle part, extracellular region, macromolecular complex and cell junction were significantly enriched for Nipponbare and Gaoshan 1. In the molecular function category, the top four function terms involving most DEGs were binding,catalytic activity, transporter activity and transcription regulator activity for both Nipponbare and Gaoshan 1.

To further explore the biological functions of the DEGs,an enrichment analysis based on KEGG pathway database was performed. The top 20 KEGG pathways involving most DEGs were listed (Fig. 5). For Nipponbare, the DEGs were mostly enriched in the terms such as“Phenylpropanoid biosynthesis,Biosynthesis of amino acids, Carbon metabolism, Ribosome and Plant hormone signal transduction”(Fig.5-A);whereas for Gaoshan 1, the DEGs were mostly enriched in “Protein processing in endoplasmic reticulum, Plant-pathogen interaction, Ubiquitin mediated proteolysis, Spliceosome and Phenylpropanoid biosynthesis” (Fig. 5-B). For the two cultivars,the common enriched pathways were:Phenylpropanoid biosynthesis, Protein processing in endoplasmic reticulum,Carbon metabolism, Plant hormone signal transduction,Plant-pathogen interaction, Starch and sucrose metabolism,MAPK signaling pathway-plant,and Endocytosis.

Fig.2- Changes in gene expression of Nipponbare and Gaoshan 1 under AWD and CF irrigation.(A)Principal component analysis(PCA)for the transcriptome samples.(B)Number of up-and down-regulated DEGs in the two cultivars for AWD vs.CF transcriptome comparison.(C)Venn diagram showing the intersection of up-and down-regulated DEGs sets as shown in B.CF,continuous flooding;AWD,alternate wetting and drying irrigation.

3.5. Lignin biosynthetic pathway

We obtained phenylpropanoid biosynthesis as one of the enriched pathways in both cultivars in response to AWD vs.CF comparison. Lignin is one of the important products of phenylpropanoid metabolic pathway and root morphological changes as observed in the cultivars in response to AWD may be attributed to the genes displaying expression changes to it. In our transcriptome,65 and 4 DEGs involved in the process of lignin biosynthesis and encoding enzymes such as cinnamoyl-CoA reductase (CCR), cinnamyl alcohol dehydrogenase (CAD), peroxidase (PRX), caffeoyl shikimate esterase (CSE), Hydroxycinnamoyl-CoA quinate/shikimate hydroxycinnamoyltransferase (HCT) and caffeoyl-CoA Omethyltransferase(CCoAOMT)were detected in the roots of Nipponbare and Gaoshan 1, respectively (Fig. 6). Two DEGs and one DEG encoding CCR were up-regulated in the roots of Nipponbare and Gaoshan 1, respectively. One DEG and two DEGs encoding CAD were up- and down-regulated in the roots of Nipponbare, respectively. One DEG encoding CAD was up-regulated in the roots of Gaoshan 1. Five DEGs and 52 DEGs encoding PRX were up- and down-regulated in the roots of Nipponbare, respectively. Two DEGs encoding PRX were up-regulated in the roots of Gaoshan 1. One DEG encoding HCT,one DEG encoding CSE and one DEG encoding CCoAOMT were down-, up- and down-regulated respectively in the roots of Nipponbare. On the other hand, no significant DEGs encoding HCT, CSE and CCoAOMT were detected in the roots of Gaoshan 1.

3.6.Protein processing in endoplasmic reticulum(ER)pathway

The quality of synthesized proteins is monitored and controlled by the ER pathway and hence variations in this may hint towards extreme changes in plant environment. Therefore, ER pathway may act as a sensor to these changing conditions.The ER quality control system appears to have been usurped to serve as an environmental sensor and stress responder in plants[43].In the study, 31 and 8 DEGs involved in protein processing in ER pathway were detected in the roots of Nipponbare and Gaoshan 1(Fig. 7), respectively. There were 24 and 6 DEGs related to ERassociated degradation (ERAD) for Nipponbare and Gaoshan 1,respectively (Fig. S1), and four Heat Shock Protein 20 (HSP20)family genes (Os03g0245800, Os03g0266300, Os03g0266900,Os03g0267000) which participated in ERAD were up-regulated for both varieties. In addition, genes encoding eukaryotic translation initiation factor 2-alpha kinases (EIF2AKs) were down-regulated for both varieties,and 4 genes encoding protein disulfide-isomerases(PDIAs)were down-regulated in the roots of Nipponbare (Fig. 7-A), while one gene encoding UDP-glucose:glycoprotein glucosyltransferase (HUGT)was down-regulated in the roots of Gaoshan 1(Fig.7-B).

3.7.Plant hormone signal transduction pathway

Abscisic acid receptor pyrabactin resistance (PYR)/pyrabactin resistance like (PYL)/regulatory component of ABA receptor(RCAR), protein phosphatase 2C (PP2C) and serine/threonineprotein kinase SRK2(SNRK2)play essential roles in ABA signal transduction pathway [44]. In this study, one DEG encoding PYR/PYL/RCAR, one DEG encoding PP2C and one DEG encoding SNRK2 were detected in the roots of Nipponbare(Fig.8-A),and one DEG encoding PYR/PYL/RCAR were detected in the roots of Gaoshan 1 (Fig. 8-B). Genes involved in cytokinin (CK) signal transduction pathway were also differently regulated in both varieties (Fig. 8). Genes encoding histidine-containing phosphotransfer protein (AHP) and twocomponent response regulator ARR-A family (ARR-A) were down-regulated and two-component response regulator ARRB family(ARR-B)were up-regulated in the roots of Nipponbare(Fig.8-A),while gene encoding ARR-A was up-regulated in the roots of Gaoshan 1 (Fig. 8-B). Interestingly, DEGs involved in brassinosteroid (BR) and auxin signal transduction pathway were specifically detected for Nipponbare roots(Fig.8-A),and most of these genes were down-regulated.

Fig.6-Expression changes for genes involved in protein processing in ER pathway in the roots of Nipponbare(A)and Gaoshan 1(B)when subjected to AWD irrigation compared to CF.For each rectangular shape,red shapes indicate up-regulated DEGs,while green shapes indicate down-regulated DEGs.ER,endoplasmic reticulum;ERAD,ER-associated degradation;PDIAs,protein disulfide-isomerases;EIF2AKs;eukaryotic translation initiation factor 2-alpha kinases;HUGT,UDP-glucose:glycoprotein glucosyltransferase;CF, continuous flooding;AWD,alternate wetting and drying irrigation.

3.8. Changes in the expression profile of transcription factors(TFs)

TFs play important regulatory roles in plant signaling responses [45]. In this study, the AWD treatment led to a number of TFs being differentially expressed in roots (Fig. 9-A).In total,328 differentially expressed TFs were identified for Nipponbare (164 TFs were up-regulated and 164 TFs were down-regulated), and 26 differentially expressed TFs (20 TFs were up-regulated and 6 TFs were down-regulated) were identified for Gaoshan 1. These TFs belonged to bHLHs (basic helix-loop-helix), bZIPs (basic region-leucine zipper), C2H2s and C3Hs (C2H2 and C3H zinc-finger proteins), MYB (v-myb avian myeloblastosis viral oncogene homolog), NACs (NAM/ATAF/CUC), WRKYs (WRKY proteins) and HSF (heat shock transcription factor) family. In the Venn intersections of TF DEGs in Nipponbare and Gaoshan 1 (Fig. 9-B), we obtained 6 common up-regulated TFs (Os03g0745000, HSF family;Os06g0565200, HSF family; Os08g0546800, HSF family;Os06g0258500, MYB_related; Os12g0634500, bZIP family and Os02g0606200, DBB family) and 1 common downregulated TF(Os03g0388600, MYB family).

3.9.Expression changes from RNA-Seq and qRT-PCR strongly correlated for above mentioned candidate genes

From the above-discussed genes belonging to lignin biosynthetic pathway, hormone signal transduction pathway and transcriptional regulation,some genes were selected for gene expression verification. To validate the expression fold changes obtained in our RNA-Seq dataset, 18 selected genes from lowland rice and upland rice were analyzed by qRT-PCR(Fig. S2). The expression fold changes from qRT-PCR strongly correlated with those obtained from RNA-Seq,with a correlation coefficient of R=0.91,indicating that the RNA-Seq data is reliable enough to plan for future functional studies utilizing these candidate genes.

Fig.7-Expression changes for genes involved in Lignin biosynthetic pathways in the roots of Nipponbare(A)and Gaoshan 1(B)when subjected to AWD irrigation compared to CF.For each rectangular shape,red shapes indicate up-regulated DEGs,while green shapes indicate down-regulated DEGs.Abbreviations:PTAL,Phe and Tyr ammonia lyase;PAL,Phe ammonia lyase;C4H,cinnamate 4-hydroxylase;4CL,4-coumarate CoA ligase;HCT,Hydroxycinnamoyl-CoA quinate/shikimate hydroxycinnamoyltransferase;CSE,caffeoyl shikimate esterase;CCR,cinnamoyl-CoA reductase;CCoAOMT,caffeoyl-CoA Omethyltransferase;COMT,caffeic acid O-methyltransferase; CAD,cinnamyl alcohol dehydrogenase;PRX,peroxidase. CF,continuous flooding;AWD,alternate wetting and drying irrigation.

4. Discussion

4.1. Response of root system to AWD irrigation

Roots play critical roles in plant growth and development,such as anchor plants in the soil, absorb and transport water and nutrients [46,47]. The root system, being the plant organ directly in contact with the soil,is the first line of defense for maintaining plant productivity under soil abiotic stresses[48].The roots are primarily responsible for the adaptation and responses to various stress situations through complex interactions between the genes [49]. A long and thick root system,the ratio of root to shoot weight and root penetration ability of upland rice have been reported to contribute to yield under water deficit conditions[50].Zhang et al.[10]concluded that moderate AWD (rewatered when soil water potential reached −15 kPa) can enhance rice root growth and improve grain yield. In this study, RDW and ROA were higher under AWD irrigation for both lowland rice and upland rice cultivar when compared with CF irrigation(Fig.1),which is consistent with the results of previous studies [26,51]. We hypothesized that a stronger root system was induced to acquire enough water and nutritional elements from soil when rice was under AWD irrigation.When compared with Nipponbare,Gaoshan 1 had higher RDW under both irrigation regimes, and higher ROA under AWD irrigation, which indicated Gaoshan 1 may have a better adaptability to AWD irrigation than Nipponbare.In a downstream transcriptome analysis, 6597 DEGs in roots of Nipponbare and 526 DEGs of Gaoshan 1 were detected(Fig.2-A). Dehydration avoidance refers to the case that plants constitutively avoid excess loss of water through more developed root systems [52]. We hypothesized that more developed root systems in Gaoshan 1 helped to get water in deep layers of soil when plants were under AWD treatment,which may explain why less DEGs were enough to cope with the changes of soil water environment under AWD irrigation.

4.2. Response of protein processing in ER pathway to AWD irrigation

Proteins are modified and folded in the ER, and one-third of newly synthesized proteins are misfolded [53]. Misfolded proteins can be folded with the assistance of ER chaperones,or destroyed by ERAD [54]. Abiotic stress can induce the accumulation of unfolded proteins in the ER and cause ER stress[55].ER response pathway was reported to regulate root growth when Arabidopsis were under water stress [56]. In our work, 24 and 6 genes involved in ERAD were differentially expressed in the roots of Nipponbare and Gaoshan 1 (Fig. 6),respectively, and four common HSP20 family genes involved in ERAD were up-regulated in the roots of both rice varieties.Genes involved in ER stress recovery (EIF2AKs) were downregulated in both rice varieties, and 4 genes involved in protein targeting (PDIAs) were specifically down-regulated in Nipponbare, while one gene involved in protein correctly folded (HUGT) was specifically down-regulated in Gaoshan 1.These showed that protein processing in ER pathway may be involved in the AWD irrigation adaption of rice roots.

Fig.8-Expression changes for genes involved in the hormone signal transduction pathways in the roots of Nipponbare(A)and Gaoshan 1(B)when subjected to AWD irrigation compared to CF.For each rectangular shape,red shapes indicate up-regulated DEGs,while green shapes indicate down-regulated DEGs.PYR/PYL,abscisic acid receptor PYR/PYL family;PP2C,protein phosphatase 2C;SNRK2,Sucrose non-fermenting-1-related protein kinase 2; ABF,ABA responsive element binding factor;AUX1,auxin influx carrier;TIR1,transport inhibitor response 1;AUX/IAA,auxin-responsive protein IAA;ARF,auxin response factor;BAK1,brassinosteroid insensitive 1-associated receptor kinase 1; BRI1,protein brassinosteroid insensitive 1; BSK,BRsignaling kinase;BZR1/2,brassinosteroid resistant 1/2;TCH4,xyloglucan:xyloglucosyl transferase TCH4;CYCD3,cyclin D3;CRE1,arabidopsis histidine kinase 2/3/4(cytokinin receptor);AHP,histidine-containing phosphotransfer protein;ARR-A,twocomponent response regulator ARR-A family;ARR-B, two-component response regulator ARR-B family;CF,continuous flooding;AWD,alternate wetting and drying irrigation.

Fig.9- Proportion of differentially expressed TF genes.(A)Number of TF genes in the up-and down-regulated DEGs of Nipponbare and Gaoshan 1 in the AWD vs.CF comparison.(B) Venn diagram showing the intersection of up-and downregulated TF gene sets as shown in A. CF,continuous flooding;AWD,alternate wetting and drying irrigation.

4.3. Response of lignin biosynthetic pathway to AWD irrigation

Phenylpropanoids are a series of secondary metabolites derived from phenylalanine, and work as structural and signaling molecules [57,58]. Lignin is one of the most important products of phenylpropanoid metabolic pathway.Lignin is a complex aromatic polymer primarily composed of p-hydroxyphenyl(H)-,guajacyl(G)-and syringyl(S)-units,and it is the main component of secondary cell walls and provides structural integrity and stiffness to plant body [59,60]. Lignin biosynthesis can be induced by various biotic and abiotic stress conditions, such as pathogen infection, low temperatures, water deficit, light, UV-B radiation and mineral deficiency [61,62]. Drought can induce tight and thick cell walls in Norway spruce roots,and more lignin increases the mechanical strength of cell walls, and cell wall lignification reduces water loss and cell dehydration[63,64].Lignin biosynthesis in the root of apple improves drought adaptation[65].In this study,genes encoding the key enzymes(CCR,CAD,PRX and so on)involved in the biosynthesis of lignin were differentially expressed in the roots of both Nipponbare and Gaoshan 1(Fig.7),indicating that lignin biosynthetic pathway participated in the regulation of AWD irrigation adaption in both rice varieties.

4.4. Response of plant hormones signal transduction pathways to AWD irrigation

Natural plant hormones including ABA, auxin, BR, CK,ethylene, gibberellins and more, act as signaling molecules to play crucial roles in mediating plant defense response against biotic and abiotic stresses [66,67]. Plant hormones were also reported to regulate root growth and development[30,56,68-70]. In this study, DEGs encoding PYR/PYL/RCAR,PP2C, SNRK2 (ABA signal transduction pathway) and AHP,ARR-A, ARR-B (CK signal transduction pathway) were detected in the roots of Nipponbare (Fig. 8-A), and DEG only encoding PYR/PYL and ARR-A were detected in the roots of Gaoshan 1 (Fig. 8-B). DEGs involved in BR and auxin signal transduction pathway were also detected in the roots of Nipponbare, while not in the roots of Gaoshan 1. The results showed that plant hormone signaling pathways may contribute to AWD related drought resistance.Nipponbare regulated more phytohormone signaling genes than Gaoshan 1 to adapt AWD condition. ABA can regulate root growth and this regulation interfaces with auxin and brassinosteroids [65].Phosphorylation promotes ARR-A and ARR-B function, respectively repressing and activating the CK response. ARR-B can also intensify Aux/IAA transcription,thus reinforcing the Aux/IAA feedback loop[71].This in turn may partially explain the increased root growth in response to AWD in the lowland cultivar Nipponbare(30%)as compared to the upland Gaoshan 1 (20%).

4.5. Response of TFs to AWD irrigation

TFs regulating drought-responsive gene transcription have been identified by previous studies, such as those from bHLH,MYB, bZIP, NAC and WRKY families [72,73]. bHLH TFs were also involved in root cell elongation and regulation of epidermal cell fate in roots [74,75]. Some MYB genes are involved in drought response through their regulation of lateral root growth, such as AtMYB77 protein is involved in promoting lateral root growth through interaction with the PYL8 ABA receptor [76]. The ABA-dependent bZIP TFs play central roles in regulating stress tolerance, and bZIP11 TF was reported to function as a regulator of primary root growth[77,78]. Transgenic rice plants over-expressing OsNAC2 and OsNAC3 TFs displayed high tolerance to dehydration, high salinity and/or cold stress, and OsNAC6 TF was involved in root structural adaptions and nicotinamide biosynthesis for drought tolerance [79-81]. Overexpression of OsWRKY30 in rice plants has been found to enhance drought tolerance and root growth regulation [82-84]. In our study, 328 and 26 TFs were differentially expressed in the roots of Nipponbare and Gaoshan 1 when they were under AWD vs. CF treatment(Fig. 9-A), respectively. Six TFs (Os03g0745000,

Os06g0565200, Os08g0546800, Os06g0258500, Os12g0634500,and Os02g0606200) and one TF (Os03g0388600) were upregulated and downregulated in both rice varieties, respectively. These TFs (especially the common differentially expressed TFs) may participate in the regulation of root adaption to AWD irrigation.

4.6. Common up- and down-regulated DEGs in Nipponbare and Gaoshan 1

A GO enrichment analysis was carried out for common upregulated (129; Table S2) and common down-regulated (25;Table S3) DEGs shared between Nipponbare and Gaoshan 1.Genes involved in the carbohydrate and cell wall polysaccharide metabolism were highly enriched in the common upregulated DEGs for Nipponbare and Gaoshan 1 (Table S4). No significant terms were obtained for the common downregulated DEGs for Nipponbare and Gaoshan 1. This contained many chitinase, glycoside hydrolase and peptidase enzyme encoding genes (Table S5).

Chitinases have been implicated in the root growth in Arabidopsis and rice previously. An Arabidopsis gene encodes chitinase-like protein (AtCTL1) when mutated reduced root growth significantly [85]. Rice homolog of this gene called OsCTL1 provides mechanical strength to tissues by inducing cellulose content. It was shown that OsCTL1 mutation causes easy breakage of tissues [86]. It was found that chitinase-like protein encoding gene (OsCLP) in rice, when mutated results in severe root growth reduction while overexpression improved the growth [87]. It is highly probable that the cell wall polysaccharide genes mentioned in Table S5 get up-regulated in both upland and lowland rice to promote root growth in AWD irrigation perhaps to ensure better water acquisition.

Individual look into the two gene lists suggested candidate genes which might promote AWD led root growth.Os10g0177200, OsDSR-1 mutation results in reduced drought tolerance while overexpression leads to increased sensitivity to ABA and improved drought tolerance possibly by scavenging ROS species[88].It was also reported in promoting lateral root formation in the context of potassium level by acting as a calcium sensor[89].Os05g0560900,OsGA2ox8 may be involved in regulating GA concentrations during crown root initiation and is also expressed in lateral roots [90]. Os03g0245800,OsHSP26 gene confers tolerance to oxidative and heat stress[91]. Os05g0202800, OsMT3b may regulate root growth as well by scavenging ROS as another gene OsMT2b when overexpressed displayed large adventitious and big lateral roots while silencing lead to serious reduction in root growth[92,93]. These group of metallothionein genes act as ROS scavengers and are downregulated in waterlogged conditions such as those in CF irrigation[94].Os04g0653600,OsPLT1 gene is a gene belonging to the PLETHORA gene family which has been investigated deeply for their role in root formation in Arabidopsis [95]. This such, a general root growth regulator maybe controlled at gene expression level to modulate root biomass. Detailed tissue-specific investigations maybe conducted for AWD vs. CF irrigated roots to validate its potential role. Finally, Os02g0606200, OsBBX4 may be regulating light dependent root growth as has been shown for rice and for Arabidopsis BBX4 gene[96,97].Altogether these genes make a preliminary list of candidates which can be pursued in functional studies in the future to enhance root biomass and hence water uptake while using less water for irrigation.

5. Conclusions

We found that RDW and ROA increased when Nipponbare(lowland rice cultivar) and Gaoshan 1 (upland rice cultivar)were irrigated under AWD regime. Gaoshan 1 had a higher RDW under CF and AWD regimes, and higher ROA under AWD irrigation than Nipponbare. The results of transcriptome analysis (AWD vs. CF) showed that 6597 and 526 genes were differentially expressed in the root of Nipponbare and Gaoshan 1, respectively. Theses DEGs play roles in lignin biosynthesis, phytohormone signal transduction (ABA, CK and so on) and transcriptional regulation (such as bHLHs,bZIPs, NACs, WRKYs, and HSF families), which may contribute to root phenotypic adaptation to AWD irrigation. The results of this work may contribute to explain the similarities and differences in AWD irrigation adaption between the two rice cultivars.

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2020.01.007.

Acknowledgments

This work was supported by the National Key Research and Development Program of China (2017YFE0118100 and 2018YFD02003025), the National Natural Science Foundation of China (31761130073 and 31872169), the China Postdoctoral Science Foundation (2019M663122), the Shenzhen Virtual University Park Support Scheme to Shenzhen Research Institute of the Chinese University of Hong Kong, the Natural Science Foundation of Hunan Province (2019JJ50263) and the Hong Kong Research Grant Council(AoE/M-05/12,AoE/M-403/16,CUHK14122415,14160516, and 14177617).

Author contributions

Jianhua Zhang and Weifeng Xu designed experiments. Tao Song, Debatosh Das, Feng Yang, Yuan Tian, Chaolin Cheng and Chao Sun performed experiments. Tao Song, Debatosh Das, Feng Yang and Moxian Chen analyzed data. Tao Song,Debatosh Das and Feng Yang wrote the manuscript. Jianhua Zhang and Weifeng Xu critically commented and revised it.