Transcriptomics profiling in response to cold stress in cultivated rice and weedy rice

Accepted Manuscript Transcriptomics profiling in response to cold stress in cultivated rice and weedy rice

Shixin Guan, Quan Xu, Dianrong Ma, Wenzhong Zhang, Zhengjin Xu, Minghui Zhao, Zhifu Guo PII: DOI: Reference:

S0378-1119(18)31112-0 https://doi.org/10.1016/j.gene.2018.10.066 GENE 43322

To appear in:

Gene

Received date: Revised date: Accepted date:

27 June 2018 24 August 2018 24 October 2018

Please cite this article as: Shixin Guan, Quan Xu, Dianrong Ma, Wenzhong Zhang, Zhengjin Xu, Minghui Zhao, Zhifu Guo , Transcriptomics profiling in response to cold stress in cultivated rice and weedy rice. Gene (2018), https://doi.org/10.1016/ j.gene.2018.10.066

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Transcriptomics profiling in response to cold stress in cultivated rice and weedy rice Shixin Guan1,2#, Quan Xu1#, Dianrong Ma1, Wenzhong Zhang1, Zhengjin Xu1, Minghui Zhao1*, Zhifu Guo1,2*

PT

1. Rice research Institute, College of Agronomy, Shenyang Agricultural University, Shenyang

RI

110866, China

SC

2. Key Laboratory of Agricultural Biotechnology of Liaoning Province, College of Biosciences

NU

and Biotechnology，Shenyang Agricultural University, Shenyang, Liaoning 110866, China

*Corresponding authors:

MA

Zhifu Guo, E-mail: [email protected], Phone: +86-24-88487164, Fax: +86-24-88492799

D

Minghui Zhao, E-mail: [email protected], Phone: +86-24-88487183, Fax: +86-24-88487183

AC

CE

PT E

#The first two authors contributed equally to this paper.

1

ACCEPTED MANUSCRIPT

Abstract Weedy rice is an important germplasm resource for rice improvement because it has useful genes for many abiotic stresses including cold tolerance. We identified the cold tolerance and cold sensitivity of two weedy rice lines (WR 03-35 and WR 03-26) and two cultivated rice lines

PT

(Kongyu 131 and 9311). During the seedling stage of these lines, we used RNA-seq to measure changes in weedy rice and cultivated rice whole-genome transcriptome before and after cold

RI

treatment. We identified 14,213 and 14,730 differentially expressed genes (DEGs) in cold-tolerant

SC

genotypes (WR 03-35, Kongyu 131), and 9219 and 720 DEGs were observed in two cold-sensitive genotypes (WR 03-26, 9311). Many common and special DEGs were analyzed in

NU

cold-tolerant and cold-sensitive genotypes, respectively. Some typical genes related to cold stress such as the basic helix-loop-helix (bHLH) gene and leucine-rich repeat (LRR) domain gene etc.

MA

The number of these DEGs in cold-tolerant genotypes is more than those found in cold-sensitive genotypes. The gene ontology (GO) enrichment analyses showed significantly enriched terms for biological processes, cellular components and molecular functions. In addition, some genes related

D

to several plant hormones such as abscisic acid (ABA), gibberellic acid (GA), auxin and ethylene

PT E

were identified. To confirm the RNA-seq data, semi-quantitative RT-PCR and qRT-PCR were performed on 12 randomly selected DEGs. The expression patterns of RNA-seq on these genes

CE

corresponded with the semi-quantitative RT-PCR and qRT-PCR method. This study suggests the gene resources related to cold stress from weedy rice could be valuable for understanding the

AC

mechanisms involved in cold stress and rice breeding for improving cold tolerance.

Key Words: Weedy rice; Cold stress; RNA-seq; Transcriptome analysis; DEGs

2

ACCEPTED MANUSCRIPT

1.Introduction Rice (Oryza sativa L.) is the most important staple food in the world and nearly half of the world’s population depends on rice as a staple food. As a thermophilic cereal, rice is highly sensitive to low temperature stress. Thus, cold stress directly affects the growth, production, and

PT

distribution of rice. Rice increases its freezing tolerance in response to low temperatures, the phenomenon known as cold acclimation (Thomashow 1999). Comprehensive studies of the

RI

molecular mechanisms for the response to low temperature are of fundamental importance to rice

SC

biology.

Domestication served as an important instance of severe phenotypic divergence in response

NU

to selection for a long time. Plant domestication is the genetic modification of a wild species to create a new form of a plant altered to meet human needs (Doebley et al. 2006). The genetic basis

MA

of domestication-associated characteristics has been identified in many organisms, involving rice, maize (Zea mays L.), wheat (Triticum aestivum L.) and tomato (Solanum lycopersicum L.) and their wild or relative species (Izawa et al. 2009; Wang et al. 1999; Gegas et al. 2010; Frary et al.

D

2000; Parker et al. 2010). These studies indicate that it is important to research the genetic

PT E

differences between wild lines and cultivated lines. Weedy rice (Oryza sativa L.) is a botanical variety of rice that occurs in rice-planting areas

CE

worldwide. It is also called red rice and it competes aggressively with cultivated species and produces low spikelet fertility. A few loci related to the similar phenotype between weedy rice and

AC

domesticated or wild rice which explain whether weedy rice derives from domesticated rice or wild rice (Thurber et al. 2010; Subudhi et al. 2014). Weedy rice has many different types with significantly different morphologies and possibly different sources (Vaughan et al. 2001). Some weedy rice genotypes may have a better cold adaptation mechanism because their seedlings grow faster than cultivated rice even during early planting in temperate conditions when cold stress often occurs. Cold-tolerant weedy rice can be a good genetic source for developing cold-tolerant, weed-competitive rice cultivars. (Oh et al. 2004, Bevilacqua et al. 2015; Wang et al. 2013) . However, little is known about the molecular basis of the cold-tolerance mechanisms in weedy rice. Studying the mechanisms of cold tolerance or cold adaptation in weedy rice is significant for

3

ACCEPTED MANUSCRIPT rice breeding for cold tolerance. The development of DNA sequencing technology has supported research into the effect of domestication at the whole-genome level (Hufford et al. 2012; Qi et al. 2013; Li et al. 2013). RNA sequencing technology (RNA-Seq) and digital gene expression have provided some ideal methods of identifying cold tolerance genes, their expression profiles and functional genomic data in some

PT

cultivated lines and wild lines (Maia et al. 2017; Zhang et al. 2017). Transcriptome analysis in cold-tolerant rice has identified 121 genes that were induced by

RI

treatment with cold stress, and found that a reactive oxygen species (ROS)-bZIP1 regulon plays an important role in early responses to cold stress (Wang et al. 2017). Through comparing the

SC

transcriptomes of the cold-tolerant rice LTH and the cold-sensitive rice IR29, the expression of many stress-inducible genes was observed to be upregulated in the cold-tolerant genotype in

NU

response to cold treatment. Further analyses showed that the C-repeat binding factor (CBF) and MYBS3 regulons are associated with tolerance to cold stress in rice (Zhang et al. 2012). Yang et

MA

al. (2015) performed a comparative transcriptome analysis of the cold-tolerant rice TNG67 and the cold-sensitive rice TCN1 in response to low temperatures. The cold-tolerant genotype had more

D

rapid alterations in gene expression to tolerate cold stress than the cold-sensitive genotype, which

PT E

related to protein metabolism, modification, folding and defense responses. The transcription factors of TNG67, such as OsMyb, OsbHLH, OsIAA23, SNAC2 and OsWRKY1v2 may be good candidates for cold stress tolerance-related genes in rice (Yang et al. 2015). A comparative

CE

microarray data analysis in a specific cold-tolerant rice variety has revealed the complicated cross-talk between the CBF and ABA-responsive element (ABRE) regulons and the role of

AC

abscisic acid (ABA) signaling in rice cold tolerance (Hossain et al. 2010). Shen et al. (2014) compared two cold-tolerant wild rice cultivars (Dongxiang wild rice and Chaling wild rice) and other cultivated genotypes via RNA-seq technology under cold stress. The results showed that 318 differentially expressed genes (DEGs) were found to be common DEGs related to cold tolerance in the wild, cold-tolerant rice. Further, it was suggested that calcium signal transduction may play a dominating role in the cold stress response of rice and RNA helicases may be directly involved in temperature sensing in rice under cold stress (Shen et al. 2014). Previous transcriptome studies of cold responses in rice were performed with cultivated rice and a small quantity of wild rice, ignoring weedy rice, which is known to be highly related to cold tolerance. It is imperative to 4

ACCEPTED MANUSCRIPT study the genetic mechanisms of low temperature tolerance using weedy rice. In the present study, we used the Illumina sequencing platform to sequence the transcriptomes of one cold-tolerant and one cold-sensitive weedy rice genotypes, and one cold-tolerant and one cold-sensitive cultivated rice genotypes, and then compared the transcriptome differences between the weedy rice and cultivated rice strains The DEGs were

PT

screened and identified for the weedy rice and cultivated rice. Then pathway enrichment analysis of the DEGs and gene ontology (GO) enrichment analysis were used to indicate the biological

RI

processes involved in the response of rice to cold stress. The goals of this study were to identify and characterize potential candidate genes involved in cold stress tolerance in rice. This should

SC

contribute to our understanding of the cold response mechanisms in rice and weedy rice, and

NU

therefore assist efforts to improving cold tolerance in rice in future breeding programs.

MA

2.Materials and methods

2.1 Plant materials and temperature treatments

D

Four rice genotypes with contrasting sensitivities to cold were used in this study. These

PT E

genotypes included the cold-tolerant weedy rice 03-35 (named CTWR1 in the results), the cold-sensitive weedy rice 03-26 (named CTWR2 in the results), the cold-tolerant cultivated rice Kong-Yu 131 from Japan (named CTR in the results) and the cold-sensitive cultivated rice 9311

CE

(named SCR in the results). Mature non-dormant seeds of the four rice lines were soaked at 28°C for 48 h and then planted in the same plastic container filled with soil at the same spacing.

AC

Seedlings were grown at 30°C until the three-leaf stage. The three-leaf seedlings were transferred to a growth chamber and maintained at a temperature of 4°C (treatment) or 30°C (control) with a 12 h light:12 h dark photoperiod as the temperature treatment. After 48 h of treatment, leaf samples containing both treated and control samples were harvested with three technological replicates and two biological replicates. Each repetition of 0.1 g fresh leaves were sampled and cut into 2-cm pieces. The leaf samples were then rinsed in pure water to remove cellular proteins from the cut ends and placed in tubes with 10ml distilled water. The collected samples were flash-frozen in liquid nitrogen and stored at −80°C until further use. 5

ACCEPTED MANUSCRIPT 2.2 Electrical conductivity measurement To measure the relative electrical conductivity (REC) under cold stress, 0.5 g of fresh leaves from each sample material after being treated at 4°C (treatment) and 30°C (control) for 48 h were harvested and put into 20 mL of distilled water in a tube and immersed for 3 h at 25°C. The initial conductivity of the solution (R1) and water (R2) were both measured. The solution and water were

PT

then heated to 100°C for 30 min, and the final conductivity of the solution (R3) and water (R4) were both measured at room temperature. The REC was calculated via the equation: REC (%) =

RI

[(R1－R2)/(R3－R4)]×100.

SC

2.3 RNA isolation

Total RNA from all samples was isolated with an RNAiso Plus kit (TaKaRa, Japan)

NU

according to the manufacturer’s instructions. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked with the NanoPhotometer spectrophotometer

MA

(IMPLEN, CA, USA). RNA concentration was measured with a Qubit RNA Assay Kit in a Qubit 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed with the RNA Nano

PT E

D

6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

2.4 cDNA library construction and transcriptome sequencing In total, of 3 μg of RNA per sample was used as the input material for the RNA sample

CE

preparations. The sequencing libraries were prepared with the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s recommendations; index codes were

AC

added to attribute sequences to each sample. Clustering of the index-coded samples was performed on a cBot Cluster Generation System with the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the libraries were sequenced with an Illumina sequencing platform (Hiseq 2500) and 125-bp or 150-bp paired-end reads were generated (Beijing Novogene Bioinformatics Technology Co., Ltd). The raw reads in fastq format have been deposited to NCBI (accession:SRP131818).

2.5 Reads filtration and mapping Raw reads in fastq format were processed with in-house Perl scripts and clean reads were 6

ACCEPTED MANUSCRIPT obtained by removing reads containing adapters, reads containing poly-N and low-quality reads. The Q20, Q30 and GC content of the clean reads were calculated. All the downstream analyses were based on the high quality clean reads. The reference genome and gene model annotation files were downloaded from the rice genome database (http://rice.plantbiology.msu.edu/). An index of the reference genome was built with Bowtie v2.2.3 (Langmead et al. 2009) and paired-end clean

PT

reads were aligned to the reference genome with TopHat v2.0.12 (Trapnell et al. 2012).

RI

2.6 Quantification and differential expression analysis of genes

HTSeq v0.6.1 was used to count the read numbers mapped to each gene (Anders and Huber

SC

2010). Next, the FPKM of each gene was calculated from the read count mapped to this gene, the length of the gene and the sequencing depth to quantify gene expression. The DESeq R package

NU

(v.1.18.0) was used to perform differential expression analysis and hierarchical clustering of the two conditions (two biological replicates per condition)(Anders and Huber 2010). The P-values

MA

were adjusted via Benjamini and Hochberg’s method to control the false discovery rate. An

D

adjusted P-value of <0.05 and |log2 (fold change)| > 1 were set as the criteria to filter the DEGs.

PT E

2.7 GO and Kyoto Encyclopedia of Genes and Genomes Pathway enrichment analysis

GO enrichment analysis of differentially expressed genes was performed with the GOseq R

CE

package (Young et al. 2010). GO terms with corrected P-value of <0.05 were considered to be significantly enriched in DEGs. The Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG;

AC

http://www.genome.jp/kegg, accessed 16 February 2018) annotations were assigned according to the KEGG database. KOBAS software (Peking University, Beijing, China) was used to test the enrichment of DEGs in KEGG pathways (Mao et al. 2005).

2.8 Semi-quantitative RT-PCR and Quantitative real-time PCR validation Semi-quantitative RT-PCR and quantitative real-time PCR (qRT-PCR) were performed to

examine the expression patterns shown by the RNA-seq analysis. Twelve genes were selected on the basis of their potential functions in RNA-seq. The actin gene of rice was used as the internal control in this study. Complementary DNA was diluted fivefold, and 2 µl of the diluted cDNA was 7

ACCEPTED MANUSCRIPT added in each 25µl PCR mixture. Each reaction mixture contained 12.5 µl of 2x Taq PCR Mix (Aidlab, China) for semi-quantitative RT-PCR and 12.5µl of 2x Sybr Green qPCR Mix (Aidlab, China) for qPCR. The reaction conditions of semi-quantitative RT-PCR were as follows: pre-heating at 95 °C for 5 min, 28 cycles at 94 °C for 30 s, 56 °C for 30 s and 72 °C for 30 s, and final extension at 72 °C for 10 min. PCR products were analyzed with electrophoresis in a 1.5 %

PT

agarose gel. qRT-PCR was performed with 3 min pre-heating at 94 °C, followed by 40 cycles of 94 °C for 20 s, 56 °C for 20 s, and 72 °C for 30 s. The relative expression levels were calculated

RI

via the ∆∆Ct method. Each sample was analyzed in three biological replicates. The primer

SC

sequences used for PCR are listed in Table S3.

NU

3.Results

MA

3.1 Physiological performance of the four rice strains during cold stress The REC is widely used as an indicator to evaluate membrane damage. In order to verify the cold sensitivity of four rice genotypes, the REC of seedling leaves after treatment with cold stress

D

was measured. The REC values of the weedy rice 03-35（CTWR1） and the cultivated rice

PT E

Kong-Yu 131 (CTR) were smaller than those of the weedy rice 03-26 (CTWR2) and the cultivated rice 9311 (SCR) ( Figure 1). It was verified that the cold tolerance of cold-tolerant lines (CTWR1

CE

and CTR) was better than that of cold-sensitive lines (CTWR2 and SCR), as expected. The result

AC

is consistent with previous studies (Bevilacqua et al. 2015; Feng et al. 2017; Wang et al. 2013).

3.2 Transcriptome sequencing and read assembly Eight RNA samples from eight experimental conditions (control and cold-treated materials) were prepared to construct cDNA libraries with fragments 150~200 bp in length. The cDNA libraries were then subjected to transcriptome sequencing analysis with the Illumina HiSeq 2000 platform. A total of 941 million raw reads were generated from the weedy rice and cultivated rice. The raw reads in fastq format have been deposited to NCBI (accession:SRP131818). A total of 885 million clean reads 125 bp or 150 bp in length were generated. The average Q20 and Q30 values were 94.72% and 89.61%, respectively (Table 1). We mapped the clean reads

8

ACCEPTED MANUSCRIPT to the Oryza indica reference genome with Tophat v2.0.12. The clean reads that successfully mapped to the reference genome ranged from 82.04% to 91.46% (Table S1). The reproducibility of the RNA-sequencing data was assessed for the cold-treated and the control samples, and the results showed that the two replicates were in highly consistent for 0.841 < R2 < 0.994 (Figure S1),

PT

which indicated that the reproducibility of RNA-sequencing data was reliable.

3.3 Differentially expressed genes in response to cold stress

RI

To identify DEGs in the four rice genotypes, the fragment number for each gene was normalized to fragments per kb of transcript sequence per million base pairs sequenced

SC

(FPKM)(Trapnell et al. 2010). The fold change of DEGs was calculated using the FPKM value of the treatment vs. the control. In order to identify the gene expression patterns, we performed

NU

hierarchical clustering of the DEGs for the four rice lines under the cold treatment and control conditions according to the log2 (fold change) values > 1 (Figure 2).

MA

A larger number of DEGs were found in the cold-tolerant rice genotypes (CTWR1, CTR) than in the cold-sensitive rice (CTWR2, SCR) after cold treatment. The number of DEGs for

D

CTWR1, CTR, CTWR2 and SCR was 14213, 14730, 720 and 9219, respectively (Figure 3, Table

PT E

S2). For the cold-tolerant genotypes, 7315 and 7918 DEGs were upregulated in CTWR1 and CTR, respectively, whereas the number of upregulated DEGs was 554 and 5164 in the cold-sensitive genotypes CTWR2 and SCR, respectively. The number of upregulated DEGs is nearly equivalent

CE

to the number of upregulated DEGs in all genotypes except for CTWR2 (Figure 3). It is showed that the number of DEGs in cold-tolerant genotypes is more than those in cold-sensitive genotypes.

AC

Four hundred and seventy-eight common DEGs related to cold tolerance were detected in all four genotypes. These included 432 common upregulated DEGs and 46 common downregulated DEGs. A total of 4692 common DEGs were specifically shared by the two cold-tolerant rice genotypes. These included 1993 common upregulated DEGs and 2699 common downregulated DEGs (Fig. 3). In particular, 3405 and 3304 specific DEGs were detected in the cold-tolerant genotypes CTWR1 and CTR, respectively. In addition 1448 DEGs were upregulated and 1957 DEGs were downregulated in CTWR1, whereas 1883 DEGs were upregulated and 1421 DEGs were downregulated in CTR. Some typical genes related to abiotic stresses including cold stress, such as Myc-type TF with a basic helix-loop-helix (bHLH) domain (OS01G0952800, OS06G0653200, 9

ACCEPTED MANUSCRIPT OS03G0135700, OS02G0747900), leucine-rich repeat (LRR)-containing proteins with a LRR domain (OS01G0162300, OS02G0138000, OS05G0406800 and OS03G0637600), WRKY TF (OS09G0481700, OS05G0137500 and OS01G0730700), etc. In contrast, only 13 common DEGs were specifically shared by the two cold-sensitive rice genotypes.

PT

3.4 Gene annotation and functional classification The GO enrichment analysis was performed to obtain the functional information for the

RI

DEGs mentioned above. The DEGs were then assigned to three main categories in the GO database: biological processes, cellular components and molecular functions. A significantly

SC

enriched GO term was identified if the corrected P-value of <0.05. The top 30 significant GO terms for biological processes, cellular components and molecular functions were finally

NU

identified in all genotypes. In the CTWR1 comparison, 15, 13 and 2 functional terms were identified for biological processes, cellular components and molecular functions, respectively,

MA

whereas 21, 8 and 1 functional terms were observed in these three functional terms in the CTR comparison. For example, “cellular process” (GO: 0009987), “cellular metabolic process” (GO:

D

0044237) and “biosynthetic process” (GO: 0009058) in the biological process category were the

PT E

most highly enriched GO terms in the CTR and CTWR1 comparison, respectively. In contrast, there were 28, 1, and 1 functional terms for the three main GO categories in SCR. Specifically, there were 28 and 2 functional terms for biological processes and molecular functions in the

CE

CTWR2 comparison, but significant cellular components were not identified. In addition, the GO term “biosynthetic process” (GO: 0009058) in the biological process category was most highly

AC

enriched in the CTWR2 comparison (Figure 4). It is obvious that there were fewer functional terms for biological processes but more functional terms for cellular components in the cold-tolerant rice genotypes (CTR and CTWR1) than in the cold-sensitive rice genotypes (SCR and CTWR2). The KEGG enrichment analyses of the DEGs highlighted hormone signal transduction as being particularly influenced by cold temperatures. It was observed that the expression of genes related to signal transduction of auxin showed altered expression in response to cold treatment. As an example, the expression of a gene (Os02g0579000) encoding a NAC family transcription factor (NAC1) was upregulated after cold treatment in CTWR1 and CTR, respectively. In addition, the 10

ACCEPTED MANUSCRIPT expression of Os06g0191300, which encodes a MAP kinase kinase (BUD1), was downregulated under cold stress in CTR. The DEGs involved in the modulation of ABA and GA were also identified. Moreover, the transcript levels of genes involved in ROS scavenging, such as the copper chaperone (Os04g0573200) and OXIDATIVE STRESS 3 (OXS3) (Os02g0684400), were

PT

decreased in response to cold stress.

3.5 Identification of cold-stress-responsive transcription factors

RI

In order to analyze the regulatory network, we observed that some transcription factors (TFs) exhibited differential expression under cold stress, out of which, those encoding zinc finger

SC

proteins were the largest group, followed by MYB, AP2/ERF, bHLH, NAC, bZIP and WRKY proteins. There were 622, 640, 451 and 23 zinc finger TFs in CTWR1, CTR,SCR and CTWR2,

NU

respectively. In addition, several heat shock transcription factors (HSF) (12 in CTWR1, 11 in CTR, 10 in SCR and 1 in CTWR2) were found to be differentially expressed in response to cold stress.

MA

The TFs that were common to all genotypes include 20 zinc finger proteins, 11 MYB, 7 AP2/ERF, 6 NAC, 5 WRKY, 3 bHLH and 2 ZIP genes, as well as one HSF gene. In these common TFs, 54

D

were upregulated in all lines, and one was downregulated in the SCR and CTWR2 lines but

PT E

upregulated in the CTWR1 and CTR lines. There were 941, 999, 707 and 12 TFs found only in CTWR1, CTR, SCR, and CTWR2, respectively. In CTWR1, 617 TFs were upregulated and 324 TFs were downregulated; 676 upregulated and 323 downregulated TFs were detected in CTR. In

CE

SCR, 404 upregulated and 303 downregulated TFs were identified. For CTWR2, 11 TFs were upregulated and one TF was downregulated.

AC

On the whole, there were more TFs in the cold-tolerant rice genotypes (CTR and CTWR1) than in the cold-sensitive rice genotypes (SCR and CTWR2). The number of upregulated genes was more than downregulated genes in all genotypes.

3.6 Confirmation of RNA-seq data by qRT-PCR To confirm the RNA-seq data, semi-quantitative RT-PCR and qRT-PCR were performed on 12 randomly selected DEGs, all of which are predicted to be associated with cold tolerance. The functional annotations of these DEGs are listed in Table S3. As shown in Figure 5, the expression patterns of RNA-seq on these 12 genes corresponded with the semi-quantitative RT-PCR and 11

ACCEPTED MANUSCRIPT qRT-PCR method. This result indicates the validity of the RNA-seq study.

4. Discussion Cold stress is a major environmental factor limiting the productivity of rice. Therefore, it is vitally important to study the mechanisms of the response to cold stress in rice. Although breeding

PT

programs for cold-stress-tolerant rice have made some progress, the molecular and genetic mechanisms of cold response remain poorly characterized in rice, especially in weedy rice

RI

(Vinocur and Altman 2005; Sanghera et al. 2011). The weedy genotypes exhibit wide diversity in

SC

their ability to withstand cold stress. Some weedy rice genotypes have better cold tolerance than cultivated rice genotypes during early planting in low temperature condition (Oh et al.

NU

2004; .Bevilacqua et al. 2015; Wang et al. 2013, Li et al. 2017). Cold-tolerant weedy rice can be tapped for rice improvement, for example via marker-assisted breeding methods. Kongyu 131, a

MA

japonica rice variety from the cold zone, has the characteristics of cold tolerance, early maturity and wide adaptability (Feng et al. 2017).

In this study, we chose two weedy rice genotypes, CTWR1 (03-35, cold-tolerant) and

D

CTWR2 (03-26, cold-sensitive), and two cultivated rice genotypes CTR (Kong-Yu 131,

PT E

cold-tolerant) and SCR (9311, cold-sensitive), as materials. These were suitable for evaluating the transcriptomes of rice species and studying DEGs in response to cold stress treatment. According

CE

to the current guidelines for RNA-Seq research, the experiments using RNA-Seq to describe changes in transcriptomes must contain at least 3 biological replicates. However, we only

AC

performed 2 biological replicates because of the problem about amount of weed rice genotype in this study. According to the assessment of the RNA-sequencing data, the two replicates were in highly consistent for 0.841 < R2 < 0.994 (Figure S1), which indicated that the reproducibility of RNA-sequencing data was reliable. In order to further verification, we used semi-quantitative RT-PCR and qRT-PCR method to confirm. The expression patterns of RNA-seq on 12 randomly selected genes corresponded with the semi-quantitative RT-PCR and qRT-PCR method both.

4.1 Distribution and analysis of DEGs in cold-tolerant and cold-sensitive genotypes

12

ACCEPTED MANUSCRIPT In this study, approximately 941 million raw reads with lengths of 125-bp and 150-bp were obtained. We identified 14,213 and 14,730 DEGs in the two cold-tolerant genotypes CTWR1 and CTR, and 9219 and 720 DEGs were observed in the two cold-sensitive genotypes SCR and CTWR2.

These results indicate that more cold-related genes may participate in regulatory

network in cold-tolerant genotypes than in cold-sensitive genotypes. For common DEGs, 478

PT

common DEGs were detected in all four genotypes. These genes may play a basic role for resisting cold stress. We found that 4692 common DEGs were specifically shared by the two cold

RI

tolerant genotypes CTWR1 and CTR, but 3405 and 3304 specific DEGs were detected for CTWR1 and CTR. In contrast, 539 common DEGs were identified for the two cold-sensitive

SC

genotypes CTWR2 and STR, but 1408 and 29 specific DEGs were detected for STR and CTWR2. Only 13 common DEGs were specifically shared by the two cold-sensitive genotypes. It is showed

NU

that there are more common DEGs and specific DEGs in the cold-tolerant genotypes than in cold-sensitive genotypes. In these common DEGs and specific DEGs, there are some typical genes

MA

related to abiotic stresses including cold stress such as Myc-type TFs with a bHLH domain, LRR-containing proteins with a LRR domain, WRKY TFs, etc. Interestingly, these common

D

DEGs were specifically shared by the two cold-sensitive genotypes did not contain typical genes

PT E

related to abiotic stresses like those mentioned above. The GO enrichment analyses of the DEGs for all comparisons showed significantly enriched terms for biological processes, cellular components and molecular functions. It was obvious that

CE

there were fewer functional terms for biological processes but more functional terms for cellular components in the cold-tolerant rice genotypes (CTR and CTWR1) than in the cold-sensitive rice

AC

genotypes (SCR and CTWR2). Interestingly, most of the enriched GO terms for cellular components were shared between the cold-tolerant genotypes CTWR1 and CTR, but the enriched GO terms for biological processes were different. For cellular components, GO terms “intracellular part” (GO:0044424), “intracellular” (GO:0005622), “cytoplasm” (GO:0005737), “cell” (GO:0005623), “cell part” (GO:0044464), “organelle” (GO:0043226), “intracellular organelle” (GO:0043229) and “cytoplasmic part” (GO:0044444) were shared by the cold-tolerant genotypes CTWR1 and CTR. These results suggest that these terms may be essential for the cold-tolerant phenotype of CTWR1 and CTR genotypes. For biological processes, GO terms “cellular

process”

(GO:0009987),

“regulation

of

cellular

process”

(GO:0050794), 13

ACCEPTED MANUSCRIPT “nucleobase-containing compound metabolic process” (GO:0006139), “aromatic compound biosynthetic process” (GO:0019438), “heterocycle biosynthetic process” (GO:0018130), “organonitrogen

compound

metabolic

process”

(GO:1901564),

“biological

regulation”

(GO:0065007) and “nucleobase-containing compound biosynthetic process” (GO:0034654) were found only in CTR genotype. The results show that these terms may be responsible for the better

PT

cold-tolerant phenotype of CTR genotype. Thus these terms may play an important role in the various cold responses in these two comparison groups. As for the shared GO terms, the different

RI

genes and the number of genes from each term may also illuminate the lines’ different levels of cold tolerance (Li et al. 2016). Therefore, the differences in the number of genes for shared GO

4.2 Hormone signals under cold stress

NU

different cold responses of the two tolerant lines.

SC

terms and enriched GO terms in biological processes may generate the basis for the slightly

MA

Many hormones such as auxin, ABA, GA and ethylene, are known to play an important role in a number of adaptive responses to abiotic stresses in plants (Jain and Khurana 2009; Peleg and

D

Blumwald 2011). In the present study, the transcriptome data indicated that the expression of

PT E

genes related to several plant hormones involved in signal transduction pathways was altered after cold treatment.

Although the role of auxin following cold stress remains unclear, recent studies have shown

CE

that cold-stress-induced alterations in plants are closely associated with auxin responses. In Arabidopsis thaliana (L.) Heynh., some key regulatory components of cold signaling pathways

AC

are known to modulate growth and development via the control of auxin signaling (Shibasaki et al. 2009; Miura et al. 2011). Several DEGs involved in auxin signaling, such as BUD1 (Os06g0191300) and NAC1 (Os02g0579000), were detected in this study. BUD1 is known to be a negative regulator of polar auxin transport and its overexpression leads to activation of plant basal and systemic acquired tolerance (Zhang et al. 2008). AT1G56010 (a homolog of Os02g0579000), also known as NAC1, is involved in auxin signaling (Xie et al. 2000). AT1G56010’s homolog in Brassica napus L. participates in response to various biotic and abiotic stresses, such as pathogen infection, physical wounding, low temperature and drought (Hegedus et al. 2003). ABA is possibly the most frequently identified hormone in terms of integrating diverse stress 14

ACCEPTED MANUSCRIPT signals and regulating downstream stress responses (Tuteja 2007). It is known that the cold stress response may be ABA-independent in plants, but there is growing evidence to the contrary (Heino et al. 1990; Lang et al. 1994). The SnRK2 protein kinase (SAPK) family has been reported to participate in various regulatory pathways, such as ABA responsiveness (Kobayashi et al. 2004). In this study, three DEGs encoding these proteins were detected. SAPK1 (Os03g0390200) showed

PT

upregulation in SCR, CTWR1 and CTR; SAPK3 (Os10g0564500) showed downregulation in SCR and CTR. Interestingly, SAPK5 (Os04g0691100) was downregulated in CTWR1 but upregulated

RI

in CTR. In addition, we identified the enhanced expression of PP2C, a crucial component in the ABA signaling pathway. The result is identical to the study on A. thaliana (Tähtiharju and Palva

SC

2001), but differs from a similar analysis of Lilium lancifolium Thunb. (Wang et al. 2014). Therefore, whether and how ABA-mediated signaling participates in rice cold stress responses

NU

remains to be further determined.

The role of the growth-related phytohormone GA in the abiotic stress response is becoming

MA

increasingly clearly (Colebrook et al. 2014). For instance, DELLA proteins are known to repress GA signaling and can facilitate the degradation of GA via the ubiquitin proteasome (Zentella et al.

D

2007). Thus exposure of A. thaliana to cold stress facilitated the accumulation of DELLA and led

PT E

to DELLA-mediated growth restriction (Achard et al. 2008, Wang et al. 2017). The expression of DELLA protein-encoding genes were found to be upregulated in this study, which implying GA degradation may be promoted by cold stress. Therefore, it is possible that the plants slowed their

CE

growth to conserve energy as an adaptation to cold conditions. In conclusion, these results indicate that various hormonal signals participate in response to cold and regulating the expression of

AC

numerous downstream genes.

4.3 ROS signals in cold response ROS such as superoxide, hydroxyl radicals and hydrogen peroxide are produced during cold conditions (Suzuki and Mittler 2006). ROS accumulated following abiotic stresses have been found to be toxic to cellular functions at high levels (Jaspers and Kangasjärvi 2010). Repressed photosynthesis is a crucial factor that leads to ROS accumulation (Noctor et al. 2002). Plants have evolved a various set of antioxidant systems in order to prevent lethal damage resulting from ROS. Several enzymes, including catalase, glutathione transferase and superoxide 15

ACCEPTED MANUSCRIPT dismutase, have been characterized for their function in ROS scavenging. In this study, we also detected genes involved in ROS scavenging. For instance, AT1G12520, the homolog of Os04g0573200, encodes a copper chaperone for SOD1, which is upregulated in response to copper and senescence (Huang et al. 2012; Kuo et al. 2013). AT5G56550 (the Arabidopsis homolog of Os02g0684400), which encodes OXS3, is involved in heavy metal and oxidative

PT

stress (Blanvillain et al. 2009).

RI

4.4 Transcription factors involved in rice responses to cold stress

The TFs and their related regulatory genes are vital for cold tolerance and have a complex

SC

regulatory network in plants, which regulates enzymes, regulatory proteins and metabolites (Cushman and Bohnert 2000). In the current study, TF genes from eight TF families were

NU

identified, including the zinc finger, MYB, AP2/ERF, bHLH, NAC, bZIP, WRKY, and HSF families. The largest group of cold-inducible TFs belonged to the zinc finger family and was

MA

composed of 1736 members.

The expression of OsiSAP8 (Os06g0612800) has been found to be induced under cold stress

D

(Kanneganti and Gupta 2008). Expression of OsDOS (Os01g0192000) has been shown to be

PT E

induced under drought stress, although not under cold stress (Kong et al. 2006). The altered expression of zinc finger genes in our data suggests that they may be involved in responses to cold stress. OsMYB15 (Os02g0624300) was cold-induced in the SCR and CTR varieties, and its

CE

homologous gene in Arabidopsis (MYB15, At3g23250) participates in negative regulation of the cold tolerance pathway (Agarwal et al. 2006), implying their different functional regulatory roles

AC

in rice and Arabidopsis. A MYB gene, OsMYBS3 (Os10g0561400), was induced under cold stress in CTWR1, and CTR. OsMYBS3 has been shown to repress the DREB1/CBF-dependent cold signaling pathway in rice. It responds more slowly to cold stress than DREB1, revealing an additional new pathway that regulates cold adaptation in rice (Suh et al. 2010). Os01g0841500 (OsMYB3R2), identified as an upregulated gene CTWR1 and CTR, is reported to improve cold tolerance by mediating the CBF pathway and the cell cycle under cold stress (Ma et al. 2009). These results indicate that MYB family members may play a crucial role in the regulatory pathway of cold tolerance in rice. The AP2/EREBP family contains a large number of plant-specific TFs and participates in abiotic stress responses (Sharoni et al. 2011). Following cold 16

ACCEPTED MANUSCRIPT stress, 52, 67, 53 and 10 AP2/ERF genes were induced in CTWR1, CTR, SCR and CTWR2, respectively, suggesting these AP2/ERF genes play a dominant role in cold stress responses in rice. OsCBF3 (Os02g0677300) was highly upregulated in all genotypes, which is identical to previous reports that DREB/CBF contributes to the cold stress response through regulating many transcriptomic and metabolic changes (Yun et al. 2010; Dubouzet et al. 2003; Mizoi et al. 2012,

PT

Wang et al. 2017). ICE1, a bHLH protein in Arabidopsis, has been reported to specifically bind to the CBF3

RI

promoter region and to promote the expression of CBF3, which leads to improved tolerance to cold stress (Chinnusamy et al. 2003). A cold-induced bHLH gene (Os03g0741100) also identified

SC

to be induced by cold stress in all genotypes were in rice (Yun et al. 2010). The remaining bHLH TF genes were reported to be induced by cold stress for the first time in this study, and thus further

NU

analysis for functional characterization of these bHLH genes is required. NAC TF genes encode a set of plant-specific TFs involved in responses to biotic and/or abiotic

MA

stress (Zheng et al. 2009; Lin et al. 2007; Nakashima et al. 2007). Seventeen and 15 NAC genes have been found to play key roles in cold stress regulation in rice. The expression of NAC genes is

D

known to be upregulated by cold stress (Zhang et al. 2012). Our results are identical to this trend

PT E

for NAC TFs.

In plants, many bZIP TFs have been discovered to function in stress signaling, but few of these have been identified in rice. In this study, 41, 50, 35 and 3, DEGs of the bZIP genes were

CE

identified for each of CTWR1, and CTR, SCR and CTWR2 respectively. Among the differentially expressed bZIP TFs, OsABF1 (LOC_Os01g64730) and OsbZIP52 (LOC_Os06g45140) were

AC

previously identified and reported to be involved in responses to cold stress (Hossain et al. 2010; Liu et al. 2012).

Many WRKY genes have been found to be transcriptionally regulated in response to abiotic stress treatments in rice (Ramamoorthy et al. 2008; Berri et al. 2009). Three of the DEGs of the WRKY genes (Os03g0758900, Os11g0490900 and Os01g0826400) have been shown to be regulated by cold or H2O2 stress in a rice variety (Yun et al. 2010). Five WRKY genes were found in all genotypes in this study. In addition, 12, 11, 10 and 1 HSF genes showed differential expression in CTWR1, CTR, SCR and CTWR2, respectively. HSF plays an important role in heat shock responsiveness in plants (Kotak et al. 2007), but the role of the heat shock-related proteins 17

ACCEPTED MANUSCRIPT in the cold stress response of plants still requires to be explored. Weedy rice is an important resource for the genetic improvement of rice, like other wild species (Lu and Ellstrand 2014). A quantitative trait locus was identified for cold tolerance at the reproductive stage in a recombinant inbred line population involving a japonica weedy rice. It was suggested that weedy rice can be a vital resource for rice breeding programs to improve cold

PT

tolerance (Oh et al. 2004). In this study, we demonstrated that cold-tolerant weedy rice has more DEGs related to abiotic stresses than cold-sensitive rice. Many genes are closely related to cold

RI

tolerance, which could be considered as candidate genes or potential markers for cold tolerance. In

SC

addition, it has been shown that weedy rice has many desirable attributes such as higher root mass, rapid seedling growth, deeper roots etc. (Suh et al. 1997). Weedy rice is conspecific to cultivated

NU

rice and frequently closely related to its companion cultivated rice, facilitating development of populations for investigating the genetics of useful agronomic traits (Hoagland and Paul 1978;

MA

Akasaka et al. 2009; Ishikawa et al. 2005; Li et al. 2017). Thus the gene resources related to cold stress from weedy rice could be valuable for improving cold tolerance in rice breeding programs. Functional analysis of these genes related to cold stress in this study will be performed by both

5. Conclusion

PT E

D

gain- and loss-of-function analyses in the future.

CE

In this study, we provided an overview of the transcriptome of two weedy rice and two cultivated rice with different cold tolerance using an Illumina Hiseq 2500 platform under cold stress. 14213

AC

and 14730 DEGs were identified in cold-tolerant genotypes CTWR1 and CTR, while 720 and 9219 DEGs were observed in cold-sensitive genotypes CTWR2 and SCR, respectively. Meanwhile, numerous common and special DEGs with up-regulation and down-regulation were observed in cold-tolerant and cold-sensitive genotypes. These cold-responsive DEGs relate to TFs, ROS scavenging and plant hormones such as abscisic acid (ABA), gibberellic acid (GA), auxin and ethylene, which are mainly involved in biological processes, cellular components and molecular functions. Further study should be required to identify whether these DEGs are the genes responsible for exploring cold-tolerance mechanisms in the rice. Our findings provide the

18

ACCEPTED MANUSCRIPT gene resources about cold stress in weedy rice for rice breeding for improving cold tolerance.

Acknowledgements

RI

PT

We express our gratitude to the anonymous reviewers for helpful comments to improve the manuscript. This work was supported by the National Key R&D Program of China (Grant No. 2018YFD0200805), the Youth Science and Technology Innovation Personnel Training Project in Agricultural Field in Liaoning Province (Grant No. 2015038), the Scientific and Technological Project in Shenyang (Grant No. 17-231-1-34), the Natural Science Foundation of Liaoning Province (Grant No. 201602670) and the Open Projects Program of the Key Laboratory of Ministry of Agriculture and Ministry of Education of the People's Republic of China.

SC

References

Akasaka M, Ushiki J, Iwata H, Ishikawa R, Ishii T., 2009. Genetic relationships and diversity of weedy

NU

rice (Oryza sativa L.) and cultivated rice varieties in Okayama Prefecture, Japan. Breeding Science 59, 401-409.

Achard P., Gong, F., Cheminant, S., Alioua, M., Hedden, P., and Genschik, P., 2008. The cold-inducible

MA

CBF1 factor-dependent signaling pathway modulates the accumulation of the growth-repressing DELLA proteins via its effect on gibberellin metabolism. Plant Cell 20, 2117-2129 Agarwal, M., Hao, Y., Kapoor, A., Dong, C. H., Fujii, H., Zheng, X., et al., 2006. A R2R3 type MYB transcription factor is involved in the cold regulation of CBF genes and in acquired freezing

D

tolerance. J Biol Chem. 281, 37636- 37645 Biol. 11:R106

PT E

Anders, S., and Huber, W., 2010. Differential expression analysis for sequence count data. Genome Berri, S., Abbruscato, P., Faivre-Rampant, O., Brasileiro, A. C., Fumasoni, I., Satoh, K., et al, 2009. Characterization of WRKY co-regulatory networks in rice and Arabidopsis. BMC Plant Biol. 9:120

CE

Bevilacqua, C. B., Basu, S., Pereira, A., Tseng, T. M., Zimmer, P. D., and Burgos, N. R., 2015. Analysis of Stress-Responsive Gene Expression in Cultivated and Weedy Rice Differing in Cold Stress Tolerance. PLoS One. 10:e0132100

AC

Blanvillain, R., Kim, J. H., Wu, S., Lima, A., and Ow, D. W., 2009. OXIDATIVE STRESS 3 is a chromatin-associated factor involved in tolerance to heavy metals and oxidative stress. Plant J. 57, 654-665

Chinnusamy, V., Ohta, M., Kanrar, S., Lee, B. H., Hong, X., Agarwal, M., et al, 2003. ICE1: a regulator of cold-induced transcriptome and freezing tolerance in Arabidopsis. Genes Dev. 17, 1043-1054 Colebrook, E. H., Thomas, S. G., Phillips, A. L., and Hedden, P., 2014. The role of gibberellin signalling in plant responses to abiotic stress. J Exp Biol. 217, 67-75 Cushman, J. C., and Bohnert, H. J., 2000. Genomic approaches to plant stress tolerance. Curr Opin Plant Biol. 3, 117-124 Doebley, J. F., Gaut, B. S., and Smith, B. D., 2006. The molecular genetics of crop domestication. Cell 127, 1309-1321 Dubouzet, J. G., Sakuma, Y., Ito, Y., Kasuga, M., Dubouzet, E. G., Miura, S., et al., 2003. OsDREB

19

ACCEPTED MANUSCRIPT genes in rice, Oryza sativa L., encode transcription activators that function in drought-, high-saltand cold-responsive gene expression. Plant J. 33, 751-763 Feng, X., Wang, C., Nan, J., Zhang, X., Wang, R., Jiang, G., et al., 2017. Updating the elite rice variety Kongyu 131 by improving the Gn1a locus. Rice 10, 35 Frary, A., Nesbitt, T. C., Grandillo, S., Knaap, E., Cong, B., Liu, J., et al., 2000. fw2.2: a quantitative trait locus key to the evolution of tomato fruit size. Science 289, 85-88 Gegas, V. C., Nazari, A., Griffiths, S., Simmonds, J., Fish, L., Orford, S., et al., 2010. A genetic framework for grain size and shape variation in wheat. Plant Cell 22, 1046-1056 Hegedus, D., Yu, M., Baldwin, D., Gruber, M., Sharpe, A., Parkin, I., et al., 2003. Molecular

PT

characterization of Brassica napus NAC domain transcriptional activators induced in response to biotic and abiotic stress. Plant Mol. Biol. 53, 383–397

RI

Heino, P., Sandman, G., Lång, V., Nordin, K., and Palva, E. T., 1990. Abscisic acid deficiency prevents development of freezing tolerance in Arabidopsis thaliana (L.) Heynh. Theor Appl Genet. 79,

SC

801-806

Hoagland, R.E., Paul, R.N., 1978. A comparative SEM study of red rice and several commercial rice (Oryza sativa) varieties. Weed Sci. 26, 619–625

NU

Hossain, M. A., Cho, J. I., Han, M., Ahn, C. H., Jeon, J. S, An, G., et al., 2010. The ABRE-binding bZIP transcription factor OsABF2 is a positive regulator of abiotic stress and ABA signaling in rice. J Plant Physiol. 167, 1512-1520

MA

Huang, C. H., Kuo, W. Y., and Jinn, T. L., 2012. Models for the mechanism for activating copper-zinc superoxide dismutase in the absence of the CCS Cu chaperone in Arabidopsis. Plant Signal Behav. 7, 428-430

Hufford, M. B., Xu, X., Van, H. J., Pyhäjärvi, T., Chia, J. M., Cartwright, R. A., et al., 2012.

D

Comparative population genomics of maize domestication and improvement. Nat Genet. 44, 808-811

PT E

Ishikawa R., Toki N., Imai K., Sato Y.I., Yamagishi H, Shimamoto Y, Ueno K, Morishima H, Sato T., 2005. Origin of weedy rice grown in Bhutan and the force of genetic diversity. Genetic Resources and Crop Evolution 52, 395-403.

Izawa, T., Konishi, S., Shomura, A., and Yano, M., 2009. DNA changes tell us about rice domestication.

CE

Curr Opin Plant Biol. 12, 185-192

Jain, M., and Khurana, J. P., 2009. Transcript profiling reveals diverse roles of auxin-responsive genes during reproductive development and abiotic stress in rice. FEBS J. 276, 3148-3162

AC

Jaspers, P., and Kangasjärvi, J., 2010. Reactive oxygen species in abiotic stress signaling. Physiol Plant 138, 405-413

Kanneganti, V., and Gupta, A. K., 2008. Overexpression of OsiSAP8, a member of stress associated protein (SAP) gene family of rice confers tolerance to salt, drought and cold stress in transgenic tobacco and rice. Plant Mol Biol. 66, 445-462 Kobayashi, Y., Yamamoto, S., Minami, H., Kagaya, Y., and Hattori, T., 2004. Differential activation of the rice sucrose nonfermenting1-related protein kinase2 family by hyperosmotic stress and abscisic acid. Plant Cell 16, 1163-1177 Kong, Z., Li, M., Yang, W., Xu, W., and Xue, Y., 2006. A novel nuclear-localized CCCH-type zinc finger protein, OsDOS, is involved in delaying leaf senescence in rice. Plant Physiol. 141, 1376-1388 Kotak, S., Larkindale, J., Lee, U., von Koskull-Döring, P., Vierling, E., and Scharf, K. D., 2007.

20

ACCEPTED MANUSCRIPT Complexity of the heat stress response in plants. Curr Opin Plant Biol. 10, 310-316 Kuo, W. Y., Huang, C. H., and Jinn, T. L., 2013. Chaperonin 20 might be an iron chaperone for superoxide dismutase in activating iron superoxide dismutase (FeSOD). Plant Signal Behav. 8:e23074 Lang, V., Mantyla, E., Welin, B., Sundberg, B., and Palva, E. T., 1994. Alterations in Water Status, Endogenous Abscisic Acid Content, and Expression of rab18 Gene during the Development of Freezing Tolerance in Arabidopsis thaliana. Plant Physiol. 104, 1341-1349 Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L., 2009. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25

PT

Li, L. F., Li, Y. L., Jia, Y., Caicedo, A. L., and Olsen, K. M., 2017. Signatures of adaptation in the weedy rice genome. Nat Genet. 49, 811-814

RI

Li, Y. H., Zhao, S. C., Ma, J. X., Li, D., Yan, L., Li, J., et al., 2013. Molecular footprints of domestication and improvement in soybean revealed by whole genome re-sequencing. BMC

SC

Genomics. 14:579

Li, Z., Hu, G., Liu, X., Zhou, Y., Li, Y., Zhang, X., et al., 2016. Transcriptome Sequencing Identified Genes and Gene Ontologies Associated with Early Freezing Tolerance in Maize. Front Plant Sci.

NU

7:1477

Lin, R., Zhao, W., Meng, X., Min, W., and Peng, Y., 2007. Rice gene OsNAC19 encodes a novel NAC-domain transcription factor and responds to infection by Magnaporthe grisea. Plant Sci. 172,

MA

120-130

Liu, C., Wu, Y., and Wang, X., 2012. bZIP transcription factor OsbZIP52/RISBZ5: a potential negative regulator of cold and drought stress response in rice. Planta 235, 1157-1169 Lu, B.R., Ellstrand, N., 2014. World food security and the tribe Triticeae (Poaceae): Genetic resources

D

of cultivated, wild, and weedy taxa for crop improvement. Journal of Systematics and Evolution 52, 661-666

PT E

Ma, Q., Dai, X., Xu, Y., Guo, J., Liu, Y., Chen, N., et al., 2009. Enhanced tolerance to chilling stress in OsMYB3R-2 transgenic rice is mediated by alteration in cell cycle and ectopic expression of stress genes. Plant Physiol. 150, 244-256 Maia, L. C. D., Cadore, P. R. B., Benitez, L. C., Danielowski, R., Braga, E. J. B., and Fagundes, P. R. 419-429

CE

R., et al., 2016. Transcriptome profiling of rice seedlings under cold stress. Funct Plant Biol. 44, Mao, X., Cai, T., Olyarchuk, J. G., and Wei, L., 2005. Automated genome annotation and pathway

AC

identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21, 3787-3793

Miura, K., Lee, J., Gong, Q., Ma, S., Jin, J. B., Yoo, C. Y., et al., 2011. SIZ1 regulation of phosphate starvation-induced root architecture remodeling involves the control of auxin accumulation. Plant Physiol. 155, 1000-1012 Mizoi, J., Shinozaki, K., and Yamaguchi-Shinozaki, K., 2012. AP2/ERF family transcription factors in plant abiotic stress responses. Biochim Biophys Acta. 1819, 86-96 Nakashima, K., Tran, L. S., Van Nguyen, D., Fujita, M., Maruyama, K., Todaka, D., et al., 2007. Functional analysis of a NAC-type transcription factor OsNAC6 involved in abiotic and biotic stress-responsive gene expression in rice. Plant J. 51, 617-630 Noctor, G., Veljovic-Jovanovic, S., Driscoll, S., Novitskaya, L., and Foyer, C. H., 2002. Drought and oxidative load in the leaves of C3 plants: a predominant role for photorespiration? Ann Bot. 89,

21

ACCEPTED MANUSCRIPT 841-850 Oh, C. S., Choi, Y. H., Lee, S. J., Yoon, D. B., Moon, H. P., and Ahn, S. N., 2004. Mapping of Quantitative Trait Loci for Cold Tolerance in Weedy Rice. Breeding Sci. 54, 373-380 Parker, H. G., Shearin, A. L., and Ostrander, E. A., 2010. Man's best friend becomes biology's best in show: genome analyses in the domestic dog. Annu Rev Genet. 44, 309-336 Peleg, Z., and Blumwald, E., 2011. Hormone balance and abiotic stress tolerance in crop plants. Curr Opin Plant Biol. 14, 290-295 Qi, J., Liu, X., Shen, D., Miao, H., Xie, B., Li, X., et al., 2013. A genomic variation map provides insights into the genetic basis of cucumber domestication and diversity. Nat Genet. 45, 1510-1515

PT

Ramamoorthy, R., Jiang, S. Y., Kumar, N., Venkatesh, P. N., and Ramachandran, S., 2008. A comprehensive transcriptional profiling of the WRKY gene family in rice under various abiotic

RI

and phytohormone treatments. Plant Cell Physiol. 49, 865-879

Sanghera, G. S., Wani, S. H., Hussain, W., and Singh, N. B., 2011. Engineering cold stress tolerance in

SC

crop plants. Curr Genomics 12, 30-43

Sharoni, A. M., Nuruzzaman, M., Satoh, K., Shimizu, T., Kondoh, H., Sasaya, T., et al., 2011. Gene structures, classification and expression models of the AP2/EREBP transcription factor family in

NU

rice. Plant Cell Physiol. 52, 344-360

Shen, C., Li, D., He, R., Fang, Z., Xia, Y., Gao, J., et al., 2014. Comparative transcriptome analysis of rna-seq data for cold-tolerant and cold-sensitive rice genotypes under cold stress. J. Plant Biol. 57,

MA

337-348

Shibasaki, K., Uemura, M., Tsurumi, S., and Rahman, A., 2009. Auxin response in Arabidopsis under cold stress: underlying molecular mechanisms. Plant Cell. 21, 3823-3838 Subudhi, P. K., Singh, P. K., DeLeon, T., Parco, A., Karan, R., Biradar, H., 2014. Mapping of seed

D

shattering loci provides insights into origin of weedy rice and rice domestication. J Hered. 105, 276-287

PT E

Suh, H. S., Sato, Y. I. , and H. Morishima, H., 1997. Genetic characterization of weedy rice (Oryza sativa L.) based on morphophysiology, isozymes and RAPD markers. Theor. Appl. Genet. 94: 316–321

Suh, J. P., Jeung, J. U., Lee, J. I., Choi, Y. H., Yea, J. D., Virk, P. S., 2010. Identification and analysis of

CE

QTLs controlling cold tolerance at the reproductive stage and validation of effective QTLs in cold-tolerant genotypes of rice (Oryza sativa L.). Theor Appl Genet. 120, 985-999 Suzuki, N., and Mittler, R., 2006. Reactive oxygen species and temperature stresses: A delicate balance

AC

between signaling and destruction. Physiol Plant. 126, 45-51 Tähtiharju, S., and Palva, T., 2001. Antisense inhibition of protein phosphatase 2C accelerates cold acclimation in Arabidopsis thaliana. Plant J. 26, 461-470 Thomashow, M. F., 1999. PLANT COLD ACCLIMATION: Freezing Tolerance Genes and Regulatory Mechanisms. Annu Rev Plant Physiol Plant Mol Biol. 50, 571-599 Thurber, C. S., Reagon, M., Gross, B. L., Olsen, K. M., Jia, Y., and Caicedo, A. L., 2010. Molecular evolution of shattering loci in U.S. weedy rice. Mol Ecol. 19, 3271-3284 Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D. R, et al., 2012. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 7, 562-578 Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al., 2010. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform

22

ACCEPTED MANUSCRIPT switching during cell differentiation. Nat Biotechnol. 28, 511-515 Tuteja, N., 2007. Abscisic Acid and abiotic stress signaling. Plant Signal Behav. 2, 135-138 Vaughan, L. K., Ottis, B. V., Prazak-Havey, A. M., Bormans, C. A., Sneller, C., Chandler, J. M., 2001. Is all red rice found in commercial rice really Oryza sativa? Weed Science 49, 468-476 Vinocur, B., and Altman, A. (2005). Recent advances in engineering plant tolerance to abiotic stress: achievements and limitations. Curr Opin Biotechnol. 16, 123-132 Wang, D.Z., Jin, Y.N., Ding, X.H., Wang, W.J., Zhai, S.S., Bai, L.P., and Guo, Z.F., 2017. Gene Regulation and Signal Transduction in the ICE–CBF–COR Signaling Pathway during Cold Stress in Plants. Biochemistry (Moscow) 82(10), 1103-1117

PT

Wang, G. J., Miao, W., Wang, J. Y., Ma, D. R., Li, J. Q., and Chen, W. F., 2013. Effects of exogenous abscisic acid on antioxidant system in weedy and cultivated rice with different chilling sensitivity

RI

under chilling stress. J Agro Crop Sci. 199, 200-208

Wang, J., Yang, Y., Liu, X., Huang, J., Wang, Q., Gu, J., et al., 2014. Transcriptome profiling of the

SC

cold response and signaling pathways in Lilium lancifolium. BMC Genomics 15:203 Wang, R. L., Stec, A., Hey, J., Lukens, L., and Doebley, J., 1999. The limits of selection during maize domestication. Nature 398, 236-239

NU

Xie, Q., Frugis, G., Colgan, D., and Chua, N. H., 2000. Arabidopsis NAC1 transduces auxin signal downstream of TIR1 to promote lateral root development. Genes Dev. 14, 3024-3036 Yang, Y. W., Chen, H. C., Jen, W. F., Liu, L. Y., and Chang, M. C., 2015. Comparative Transcriptome

MA

Analysis of Shoots and Roots of TNG67 and TCN1 Rice Seedlings under Cold Stress and Following Subsequent Recovery: Insights into Metabolic Pathways, Phytohormones, and Transcription Factors. PLoS One. 10:e0131391

Young, M. D., Wakefield, M. J., Smyth, G.K., and Oshlack, A., 2010. Gene ontology analysis for

D

RNA-seq: accounting for selection bias. Genome Biol. 11,R14 Yun, K. Y., Park, M. R., Mohanty, B., Herath, V., Xu, F., Mauleon, R., et al., 2010. Transcriptional

PT E

regulatory network triggered by oxidative signals configures the early response mechanisms of japonica rice to chilling stress. BMC Plant Biol. 10,16 Zentella, R., Zhang, Z. L., Park, M., Thomas, S. G., Endo, A., Murase, K., et al., 2007. Global analysis of della direct targets in early gibberellin signaling in Arabidopsis. Plant Cell 19, 3037-3057

CE

Zhang, T., Huang, L., Wang, Y., Wang, W., Zhao, X., Zhang, S., 2017. Differential transcriptome profiling of chilling stress response between shoots and rhizomes of Oryza longistaminata using RNA sequencing. PLoS One. 12, e0188625

AC

Zhang, T., Zhao, X., Wang, W., Pan, Y., Huang, L., Liu, X., 2012. Comparative transcriptome profiling of chilling stress responsiveness in two contrasting rice genotypes. PLoS One. 7, e43274 Zhang, X., Xiong, Y., Defraia, C., Schmelz, E., and Mou, Z., 2008. The Arabidopsis MAP kinase kinase 7: A crosstalk point between auxin signaling and defense responses? Plant Signal Behav. 3, 272-274 Zheng, X., Chen, B., Lu, G., and Han, B., 2009. Overexpression of a NAC transcription factor enhances rice drought and salt tolerance. Biochem Biophys Res Commun. 379, 985-989

23

ACCEPTED MANUSCRIPT

Table 1. Summary of RNA-Seq data and sequence assembly Raw reads

Clean reads Number

%

Clean Error Q20 (%) bases (Gb) rate (%)

Q30 (%)

GC content(%)

6.96

0.03

95.33

90.68

55.79

CT2_9311 57,580,104 53,884,698 93.58

6.74

0.03

95.26

90.53

55.41

CK1_9311 56,805,258 54,220,860 95.45

6.78

0.03

95.26

90.58

53.11

CK2_9311 54,656,844 52,235,434 95.57

6.53

0.03

95.34

90.73

53.69

CT1_0326 62,170,152 57,838,054 93.03

7.23

0.03

94.91

90.00

54.06

CT2_0326 59,474,334 55,982,920 94.13

7.00

0.03

95.17

90.54

52.32

CK1_0326 54,942,754 52,492,070 95.54

6.56

0.03

95.43

90.95

51.96

CK2_0326 56,402,514 53,946,658 95.65

6.74

0.03

95.46

90.98

52.31

CT1_0335 61,178,026 57,218,358 93.53

7.15

0.03

94.09

88.61

50.71

CT2_0335 62,130,758 58,198,270 93.67

7.27

0.03

94.11

88.70

49.78

CK1_0335 52,926,568 50,324,602 95.08

6.29

0.03

94.18

88.53

52.64

CK2_0335 52,813,114 49,962,074 94.60

6.25

0.03

94.22

88.65

52.72

7.81

0.03

94.17

88.54

53.41

6.60

0.03

94.16

88.49

54.92

7.73

0.03

94.17

88.60

52.27

7.02

0.03

94.19

88.62

52.21

CT2_131 57,256,734 52,772,638 92.17

MA

CK1_131 65,553,124 61,870,192 94.38

AC

CE

PT E

D

CK2_131 59,363,864 56,124,972 94.54

RI

NU

CT1_131 67,392,984 62,452,166 92.67

PT

CT1_9311 60,041,148 55,656,884 92.70

SC

Sample

24

PT

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

SC

RI

Figure 1. The REC (relative electrical conductivity) of SCR, CTWR2, CTWR1 and CTR treated at 4°C for 48 h. Vertical bars represent standard error (P<0.05).

25

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

Figure 2. Venn map of the differentially expressed genes in the four rice genotypes under cold stress.

26

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

Figure 3. Hierarchical cluster analysis of the differentially expressed genes based on the log2 (Fold change) values. The red bands represent high gene expression levels; the blue bands represent low gene expression levels.

27

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC

CE

PT E

Figure 4. GO terms of the functional categorization of differentially expressed genes (DEGs) from SCR (A), CTWR2 (B), CTWR1 (C) and CTR (D). Significantly enriched GO terms are indicated by * for corrected a P-value of <0.05.

28

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

Figure 5. Verification of RNA-seq expression via qRT-PCR. The relative expression level of a selected gene was normalized against the rice actin 1 gene. The bars represent standard deviation.

29

ACCEPTED MANUSCRIPT Highlights

1. Some common and special DEGs (differentially expressed genes) were identified in the cold tolerance and cold sensitivity of weedy rice lines and cultivated rice lines. 2. The gene ontology (GO) enrichment analyses showed significantly enriched terms

PT

for biological processes, cellular components and molecular functions.

3. Some genes related to abscisic acid (ABA), gibberellic acid (GA), auxin and

RI

ethylene were identified.

SC

4. Gene resources related to cold stress from weedy rice could be valuable for understanding the mechanisms involved in cold stress and rice breeding for improving

AC

CE

PT E

D

MA

NU

cold tolerance.

30

ACCEPTED MANUSCRIPT Abbreviations DEGs, differentially expressed genes; bHLH, basic helix-loop-helix; LRR, leucine rich repeats; GO, gene ontology; ABA, abscisic acid; GA, gibberellic acid; RNA-Seq, RNA sequencing technology; ROS, reactive oxygen species; CBF, C-repeat binding factor; ABRE, ABA-responsive element; REC, relative electrical conductivity;KEGG, Kyoto Encyclopedia of Genes and

PT

Genomes Pathway; qPCR, quantitative real-time PCR; CTWR2, weedy rice 03-35; CTWR1, weedy rice 03-26; CTR, cultivated rice Kong-Yu 131; SRT, cultivated rice 9311;BUD1, MAP

RI

kinase kinase; NAC1, NAC family transcription factor; OXS3, OXIDATIVE STRESS 3; HSF,

SC

heat shock transcription factors; DREBs, dehydration responsive element binding factors; TFs,

AC

CE

PT E

D

MA

NU

transcription factors.

31

ACCEPTED MANUSCRIPT Authors' contributions S.G., M.Z. and Z.G. conceived and designed the experiment. Z.X., W.Z. and Z.G. conducted the experiments. M.Z. and S.G. analyzed the data. S.G. and Z.G. wrote the manuscript. Z.X. and D.M. revised the final version of the paper. All authors approved

AC

CE

PT E

D

MA

NU

SC

RI

PT

the final version of the manuscript.

32