当前位置：文档库 › The RNA-seq approach to discriminate gene expression profiles in

The RNA-seq approach to discriminate gene expression profiles in

The RNA-seq approach to discriminate gene expression pro ?les in response to melatonin on cucumber lateral root formation

Abstract:Cucumber is a model cucurbitaceous plant with a known genome sequence which is important for studying molecular mechanisms of root development.In this study,RNA sequencing was employed to explore the mechanism of melatonin-induced lateral root formation in cucumber under salt stress.Three groups of seeds were examined,that is,seeds primed without melatonin (CK),seeds primed in a solution containing 10or 500l mol/L

melatonin (M10and M500,respectively).These seeds were then germinated in NaCl solution.The RNA-seq analysis generated 16,866,670sequence reads aligned with 17,920genes,which provided abundant data for the analysis of lateral root formation.A total of 17,552,17,450,and 17,393genes were identi?ed from roots of the three treatments (CK,M10and M500,

respectively).The expression of 121genes was signi?cantly up-regulated,and 196genes were signi?cantly down-regulated in M500which showed an

obvious increase on the number of lateral roots.These genes were signi?cantly enriched in 57KEGG pathways and 16GO terms (M500versus CK).Based on their expression pattern,peroxidase-related genes were selected as the candidates to be involved in the melatonin response.Several transcription factor families might play important roles in lateral root formation processes.A number of genes related to cell wall formation,carbohydrate metabolic processes,oxidation/reduction processes,and catalytic activity also showed di?erent expression patterns as a result of melatonin treatments.This RNA-sequencing study will enable the scienti?c community to better de?ne the

molecular processes that a?ect lateral root formation in response to melatonin treatment.

Na Zhang 1,Hai-Jun Zhang 1,Bing Zhao 1,Qian-Qian Sun 1,Yun-Yun Cao 1,Ren Li 1,Xin-Xin Wu 1,Sarah Weeda 2,Li Li 3,Shuxin Ren 2,Russel J.Reiter 4and Yang-Dong Guo 1

College of Agriculture and Biotechnology,China Agricultural University,Beijing,China;2

School of Agriculture,Virginia State University,Petersburg,VA,USA;3Novogene

Bioinformatics Institute,Beijing,China;4

Department of Cellular and Structural Biology,The University of Texas Health Science Center,San Antonio,TX,USA

Key words:cucumber,gene expression,lateral root formation,melatonin,RNA sequencing Address reprint requests to Yang-Dong Guo,College of Agriculture and Biotechnology,China Agricultural University,No.2,Yuanmingyuan West Rd,Beijing 100193,China.E-mail:yaguo@https://www.wendangku.net/doc/094852731.html,

Shuxin Ren,School of Agriculture,Virginia

State University,PO Box 9061,Petersburg,VA 23806,USA.

E-mail:sren@https://www.wendangku.net/doc/094852731.html,

The ?rst two authors contribute equally to this paper.

Received September 3,2013;Accepted September 10,2013.

Introduction

Cucumber (Cucumis sativus L.),which belongs to the cuc-urbitaceous family,together with squash,loofa,bitter gourd,and muskmelon,are important horticultural crops.Cucumber o?ers various advantages as an experimental plant compared with other cucurbitaceous plants,includ-ing a small genome size and known genome sequence [1,2].The degree of root branching impacts the e?ciency of water uptake,acquisition of nutrients and anchorage by plants.A well-developed root system is essential for the vegetative growth and fruit development of plant.In recent years,the plant root system has emerged as a cen-tral focus of research in many laboratories across the https://www.wendangku.net/doc/094852731.html,teral roots (LRs)are the most dynamic and physiologically active part of the root system.To adapt to a very heterogeneous environment,root architecture is extremely plastic,responding to the soil nutrient concen-trations,soil matrix heterogeneity,and biotic interactions [3].Lateral root formation is critical for the development of root architecture.LR formation,the production of new roots from parent roots,is a complicated developmental process in higher plants.In most dicots,LRs are formed from root pericycle cells adjacent to the protoxylem poles of the parent root.A mature LR primordium is formed

through several developmental stages with well-ordered cell divisions and cell di?erentiation [4].LR formation is regulated by an intrinsic developmental program and envi-ronmental signals.Previous research primarily focused on auxin’s role on LR initiation and development [3,5,6].Recently,some researchers successively reported the e?ects of melatonin on root architecture in di?erent crops includ-ing St.John’s wort,wild leaf mustard,sweet cherry root-stocks,lupin,etc.[7–10].

Melatonin has been found to be a ubiquitous and highly conserved molecule in the plant kingdom [11].There are reports demonstrating the ability of melatonin to alleviate the e?ects of abiotic stresses such as low temperature,osmotic stress,copper stress,and light conditions during seed germination [12–15].An increase in the endogenous concentration of melatonin was correlated with a rise in de novo root formation in St.John’s wort (Hypericum perforatum cv.Anthos)[7].In many plants,melatonin has been shown to promote root growth.When exogenous melatonin was applied to etiolated seedlings of wild leaf mustard (Brassica juncea ),the results showed that melatonin had a stimulatory e?ect on root growth in a concentration dependent manner.Application of 0.1m M melatonin also raised the endogenous levels of free indole acetic acid in mustard roots [8].Similar e?ects were noted

J.Pineal Res.2014;56:39–50

Doi:10.1111/jpi.12095

?2013John Wiley &Sons A/S.Published by John Wiley &Sons Ltd

Journal of Pineal Research

M o l e c u l a r ,B i o l o g i c a l ,P h y s i o l o g i c a l a n d C l i n i c a l A s p e c t s o f M e l a t o n i n

in an experiment on sweet cherry rootstocks,where mela-tonin at a low concentration caused an auxinic response concerning the number and length of roots,but at high concentration,it was inhibitory to rooting in all the tested rootstocks[9,16].The mechanisms of this di?erential response are still unknown.

Melatonin modulates root system architecture by stimu-lating lateral root formation but minimally a?ecting pri-mary root growth or root hair development[17].LRs originate exclusively from pericycle founder cells located opposite xylem poles[18].Melatonin can induced the appearance of root primordia from pericycle cells in etio-lated hypocotyls of Lupinus albus L.,modifying the pat-tern of distribution of adventitious or lateral roots,the time course,the number and length of adventitious roots, and the number of LRs[10].

Serotonin N–acetyl transferase(NAT)is believed to be the rate-limiting enzyme in melatonin biosynthesis.Trans-genic rice seedlings expressing sheep NAT showed enhanced seminal root growth compared with wild-type rice seeds.Enhanced melatonin levels were detected in T3 homozygous seedlings[19].By microarray analysis,gene expression pro?les were altered in the transgenic rice[20]. Soil salinity is one of the most signi?cant abiotic stresses limiting plant growth,productivity,and geographical dis-tribution[21].Application of0.1l M melatonin noticeably alleviated the growth inhibition of M.hupehensis seedlings caused by high salinity[22].In our previous study,we found that melatonin had an e?ect on improving osmic stress tolerance.Under osmic stress circumstances,melato-nin-stimulated lateral root generation and showed a detectable ability on strengthening the root system[15]. Although many researches demonstrated a similar phe-nomenon,the mechanism of this e?ect remains unknown. The current study provides insight into the molecular and cellular events associated with the action of melatonin on lateral root formation.

In recent years,several di?erential screening techniques (e.g.,di?erential display,subtraction libraries,and di?er-ential hybridization)have made it possible to characterize genes that are di?erentially expressed after certain treat-ments[23–25].Recently,methods,such as serial analysis of gene expression(SAGE)and microarray analysis,have allowed us to visualize global changes in transcript abun-dance in a spatial,temporal,or conditional ways[26,27]. More recently developed RNA deep-sequencing technol-ogy provides a platform for measuring di?erences in gene expression in a manner that is more sensitive than that of traditional microarray hybridization techniques[28].This new method dramatically changes the procedures used to identify the di?erentially expressed genes[29,30].The RNA-seq method generates absolute information,rather than relative gene expression measurements,and it is more sensitive for low-expressed transcripts.Thus,the RNA-seq method avoids many of the inherent limitations of micro-array analysis.Here,we use deep RNA sequencing com-bined with digital gene expression pro?le(DGE)analysis to rapidly identify and analyze the e?ects of melatonin on cucumber LR formation.In this work,three sequencing libraries prepared from control and melatonin-treated samples were sequenced using Illumina/Solexa platform.The digital gene expression analysis was used to discrimi-nate gene expression di?erences in response to melatonin under NaCl stress in cucumber.

Materials and methods

Plant materials

Cucumber seeds were pretreated with melatonin by prim-ing in a melatonin solution at4°C for24h(0,10, 500l mol/L),that is,CK,M10,and M500,respectively. Then,the seeds were germinated in Petri dishes with three layers of?lter paper(Whatman International Ltd,Maid-stone,UK)moistened with10mL150mmol/L NaCl solution for2days.The Petri dishes were placed at25°C in darkness.Totally,100simultaneously germinated seeds were selected in each group to develop the root system. The experiments were repeated three times.In the sam-pling process,the roots and hypocotyls were separated. The experimental materials for sequencing were obtained by mixing100roots together in each group.

Staining of lateral root primordium and peroxidase activity measurement

Staining of LR primordium was processed according to De Tomasi[31].Peroxidase activity was measured accord-ing to Scebba et al.[32].The method is based on monitor-ing H2O2decomposition rate by peroxidase,using guaiacol as a hydrogen donor.The reaction was initiated by adding50l L enzyme extracts to1950l L phosphate bu?er(65m M,pH 5.5)containing11m M H2O2and 2.25m M guaiacol.The rate of color development was determined by recording the absorbance of the reaction solution at470nm every0.1s.One unit of peroxidase activity was de?ned as an absorbance change of 0.01units/s;activity was expressed as units/g protein x s. RNA preparation

Total RNA was extracted from the roots using TRIzol reagent according to the manufacturer’s protocol(Invitro-gen,Burlington,ON,Canada).Then,RNA degradation and contamination was monitored on1%agarose gels. After that,RNA purity was checked using the Nano Pho-tometer spectrophotometer(IMPLEN,Westlake Village, CA,USA).Additionally,we measured RNA concentra-tion using Qubit RNA Assay Kit in Qubit2.0Flurometer (Life Technologies,Carlsbad,CA,USA).Lastly,RNA integrity was assessed using the RNA Nano6000Assay Kit of the Bioanalyzer2100system(Agilent Technologies, Santa Clara,CA,USA).

Transcriptome sample preparation for sequencing

The total amount of3l g RNA per sample was used as input material for the RNA sample preparations.All three samples had RIN values above8.0.Sequencing libraries were generated using Illumina TruSeq RNA Sample Preparation Kit(Illumina,San Diego,CA,USA)follow-ing manufacturer’s recommendations,and three index

Zhang et al.

codes were added to attribute sequences to each sample. Brie?y,mRNA was puri?ed from total RNA using poly-T oligo-attached magnetic beads.Fragmentation was carried out using divalent cations under elevated temperature in Illumina proprietary fragmentation bu?er.First-strand cDNA was synthesized using random oligonucleotides and SuperScript II.Second-strand cDNA synthesis was subse-quently performed using DNA polymerase I and RNase H.Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities,and enzymes were removed.After adenylation of3′ends of DNA fragments, Illumina PE adapter oligonucleotides were ligated to pre-pare for hybridization.To select cDNA fragments of pref-erentially200bp in length,the library fragments were puri?ed with AMPure XP system(Beckman Coulter, Beverly,MA,USA).DNA fragments with ligated adap-tor molecules on both ends were selectively enriched using Illumina PCR Primer Cocktail in a10cycle PCR.Prod-ucts were puri?ed(AMPure XP system)and quanti?ed using the Agilent high-sensitivity DNA assay on the Agi-lent Bioanalyzer2100system.

Clustering and sequencing

The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS(Illumia,San Diego,USA) according to the manufacturer’s instructions.After cluster generation,the library preparations were sequenced on an Illumina Hiseq2000platform and100bp single-end reads were generated.

Validation of RNA-seq data by real-time quantitative PCR

Validation of RNA-seq data for12di?erent genes was per-formed by real-time quantitative PCR(Q-PCR).The prim-ers of selected genes were designed using Primer-priemer5software(PREMIER Biosoft,Palo Alto,CA, USA)and synthesized by Sangon.The primer pairs are sum-marized in Table S1.The cDNA was synthesized from1l g of total RNA using PrimeScript RT reagent Kit(Takara, Da Lian,China)in20l L of reaction mixture.The amount of the ampli?ed DNA was monitored by?uorescence at the end of each cycle using7500Real-Time PCR System (Applied Biosystems).Each plate was repeated three times in independent runs for all reference and selected genes. Gene expression was evaluated by the2àDD C t method[33]. Data analysis

Quality control

Raw data(raw reads)of fastq format were?rst processed through in-house Perl scripts.In this step,clean data (clean reads)were obtained by removing reads containing adapter,reads containing ploy-N and low-quality reads from raw data.At the same time,Q20,Q30,GC content, and sequence duplication level of the clean data were calculated.All the downstream analyses were based on the clean data with high quality.

Reads mapping to the reference genome

Reference genome and gene model annotation?les were downloaded from genome website(https://www.wendangku.net/doc/094852731.html,) directly.An index of the reference genome was built using Bowtie v0.12.8(Broad Institute,Cambridge,MA,USA) and paired-end clean reads were aligned to the reference genome using TopHat v1.4.0(Broad Institute).We selected TopHat as the mapping tool for that TopHat can generate a database of splice junctions based on the gene model annotation?le and thus a better mapping result than other nonsplice mapping tools.

Quanti?cation of gene expression level

HTSeq v0.5.3(EMBL,Heidelberg,Germany)was used to count the reads numbers mapped to each gene.And then RPKM of each gene was calculated based on the length of the gene and reads count mapped to this gene.RPKM, Reads Per Kilobase of exon model per Million mapped reads,considers the e?ect of sequencing depth and gene length for the reads count at the same time,and is cur-rently the most commonly used method for estimating gene expression levels[34].

Di?erential expression analysis

Prior to di?erential gene expression analysis,for each sequenced library,the read counts were adjusted by edgeR program package through one scaling normalized factor.

Di?erential expression analysis of two conditions was performed using the DEGSeq R package(1.12.0;

TNLIST,Beijing,China).The P-values were adjusted using the Benjamini and Hochberg method[35].Corrected P-value of0.005and log2(fold change)of?1were set as the threshold for signi?cantly di?erential expression.

GO and KEGG enrichment analysis of di?erentially expressed genes

Gene ontology(GO)enrichment analysis of di?erentially expressed genes was implemented by the GOseq R pack-age,in which gene length bias was corrected.GO terms with corrected P-value<0.05were considered signi?cantly enriched by di?erential expressed genes.

KEGG is a database resource for understanding high-level functions and utilities of the biological system,such as the cell,the organism and the ecosystem,from molecu-lar-level information,especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies(http://www.

genome.jp/kegg/).We used KOBAS software(KOBAS, Surrey,UK)to test the statistical enrichment of di?eren-tial expression genes in KEGG pathways.

Results and discussion

To test the e?ect of melatonin on lateral root formation, we pretreated cucumber seeds with a melatonin solution before germination.Then,we picked the simultaneously germinated seeds to develop the root system.Melatonin

41 DGE analysis of melatonin on lateral root formation in cucumber

pretreatment showed clearly recognizable e?ects in the root.After 2days development,CK group only showed short primary roots and few lateral roots,but in the mela-tonin-treated groups,lateral roots had emerged out of the primary root surface.And the number of lateral roots was in a concentration-dependent manner;thus,the M500plants had more lateral roots than the M10group (Fig.1).To show this e?ect more clearly,we stained the lateral roots primordia when the primary root emerged.The mel-atonin pretreated group showed more lateral roots primor-dia than control (Fig.2).M500showed the highest lateral roots primordia density.It is,therefore,evident that mela-tonin has a positive e?ect on lateral root formation under NaCl stress.

To investigate the expression patterns of genes in cucumber roots in the three samples,the RNA extracted from the roots were sequenced by Illumina Hiseq 2000platform,a total of 16,866,670(CK),15,436,994(M10),16,191,788(M500)sequence reads were generated.

An overview of the sequencing and assembly is outlined in Table 1,Figure S1.After removal of adaptor sequences,duplication sequences,ambiguous reads and low-quality reads,16,717,677high-quality clean reads (16,866,670,and 99.1%of the raw data)of CK remained.Ninety percent of the clean reads data had Phred-like quality scores at the Q30level (an error probability of 0.001).We mapped the sequenced reads to the genome.In the CK sample and melatonin,pretreated samples (M10,M500),respectively,89.5%,91.7%,and 92.4%of the total reads (15million)from RNA-seq data were mapped uniquely to the genome,while small proportions were mapped multiply to genome (Table S2)

Comparison of changes in gene expression between con-trol group and melatonin-treated plants had shown not only similarities but also considerable di?erences.Based on deep sequencing of the three DGE libraries in the current study,17,000genes (65%of reference genes in cucumber)were detected in each library.The samples exhibited transcript complexity with and without melatonin treatment.The Venn diagram shows the distri-bution of expressed genes among the three samples.We can see that not as many genes were detected in melatonin-treated samples as in control samples.Hundreds of genes were detected di?erentially expressed in three samples.Among these genes,16,711were expressed at all three sam-ples,16,994were coexpressed in CK and M10,16,932

were

Fig.1.Melatonin’s e?ect on cucumber lateral root formation under NaCl stress.CK-seeds pretreated with water;M10-seeds pretreated with 10l mol/L melatonin;M500-seeds pretreated with 500l mol/L

melatonin.

Fig.2.Staining of the lateral roots primordia.Cucumber pretreated with melatonin showed more lateral roots primordia.CK-seeds pretreated with water;M10-seeds pretreated with 10l mol/L melatonin;M500-seeds pretreated with 500l mol/L melatonin.

Zhang et al.

coexpressed in CK and M500,and16,924were coex-pressed in M10and M500.The number of melatonin con-centration speci?cally expressed genes was337(CK),190 (M10),and178(M500),respectively(Fig.3;see Data S1 for the gene list for each category).Although there were 16,711genes expressed in all three samples,many of them were quantitatively regulated.While some of these genes had little variation with or without treatment,which may be related to their housekeeping functions.These dynamic patterns of gene expression reinforce the conclusion that melatonin a?ects the lateral root formation by a highly complex process.

To obtain statistical con?rmation of the di?erences in gene expression among the developmental stages,we used the RPKM to normalize the expression level of genes.All these uniquely mapped reads were used for calculating the genes RPKM values.Genes with RPKMs in the interval0.1–3.57are considered to be expressed at low level;genes with RPKMs in the inter-val 3.57–15is considered to be expressed at medium level;genes with RPKM beyond15is considered to be expressed at high level.The percentage of high level expressed genes in control was a little larger than in melatonin-treated samples(Table2).On the contrary, the percentage of low level expressed genes in control was smaller than melatonin-treated samples.Melatonin showed an inhibitory e?ect on gene expression.To min-imize false positives and negatives,we concluded that a statistical analysis was reliable when applied to genes with an RPKM value>3.57in at least one of the three

samples.

We made a hierarchical clustering of the di?erentially

expressed genes based on the three sample’s log10RPKM,

so we could observe the gene expression pattern overall.

The blue bands identify low gene expression quantity,and

the orange represent the high gene expression quantity

(Fig.4).

The expression pro?les of the di?erentially expressed

genes(DEGs)were illustrated by a cluster analysis based

on the k-means method[36].We did this to identify simi-

lar expression patterns of all DEGs across the set of

stages,and how those patterns were similar or di?erent

between DEGs.Six expression patterns(clusters)of di?er-

entially expressed genes were identi?ed(Fig.5).The most

abundant group was subcluster3,with102genes whose

expression showed a positive slope when treated with mel-

atonin.These genes were expressed at their highest level at

M500.Subcluster5showed a similar pattern.The second

most abundant group was subcluster1,which contained

95genes that began to up-regulate at M500.Subcluster2

and subcluster4consisted of32and59genes,respectively,

whose expression showed a continuous negative slope

when treated with melatonin.Subcluster6was composed

of53genes that were down-regulated at M10and then up-

regulated at M500.

To validate the RNA-seq data,Q-PCR of12randomly

selected genes was performed.As shown in Fig.6,there

was a strong positive correlation(R2=0.97)between RNA-seq data and Q-PCR data.The Q-PCR expression

quantities were basically consistent with their transcript

abundance changes identi?ed by RNA-seq,which means

the RNA-seq data were credible.

To identify genes showing a signi?cant change in

expression during di?erent developmental stages,the

di?erentially expressed tags between two samples were

identi?ed using software of DEG-Seq and P-values were

corrected by Benjamini and Hochberg method[35].Cor-

Table1.Summary of sequence assembly after illumina sequencing.CK-samples pretreated with water;M10-samples pretreated with 10l mol/L melatonin;M500-samples pretreated with500l mol/L melatonin

Sample name Raw reads Clean reads Clean bases Error rate(%)Q20(%)Q30(%)GC content(%) CK16,866,67016,717,677 1.67G0.0396.9890.9244.15

M1015,436,99415,315,528 1.53G0.0396.9790.8944.25

M50016,191,78816,071,539 1.61G0.0397.0391.04

44.08

Fig.3.Venn diagram showing the overlaps between control sam-ples and samples treated with melatonin.CK-samples pretreated with water;M10-samples pre-treated with10l mol/L melatonin; M500-samples pretreated with500l mol/L melatonin.

Table2.Statistics of genes in di?erent expression-level interval.

CK-samples pretreated with water;M10-samples pretreated with 10l mol/L melatonin;M500-samples pretreated with500l mol/L melatonin

RPKM interval CK(%)M10(%)M500(%) 0–0.13683(15.84)3915(16.84)3930(16.90)

0.1–0.3584(2.51)524(2.25)588(2.53)

0.3–3.574327(18.61)4257(18.31)4300(18.50)

3.57–156309(27.14)6386(27.47)6226(26.78)

15–605657(24.33)5527(23.77)5524(23.76) >602688(11.56)2639(11.35)2680(11.53)

43 DGE analysis of melatonin on lateral root formation in cucumber

rected P -value of 0.005and log2(fold change)of ?1were set as the threshold for signi?cantly di?erential expression.A total of 113signi?cantly changed genes were detected between the CK and M100cucumber root libraries,with 16genes up-regulated and 97genes down-regulated based on CK.Between the CK and M500cucumber root libraries,a total of 317DEGs were detected,with 121up-regulated genes and 196down-regulated genes based on CK (Fig.7).This suggested that the di?erentiation of expressed genes between CK and M500was larger than that between CK and M100.

We adopted a RNA-sequencing approach to identify and analyze the di?erential expression of cucumber root genes associated with melatonin regulation under salt stress.The samples pretreated with 500l M of melatonin showed the most remarkable change in lateral root num-ber;subsequently,expression of the 391di?erentially expressed genes was examined (Fig.3).Among them,77genes were regulated in a melatonin-dependent manner (Fig.8).The one with highest di?erentially expressed ratio was a gene that encodes peroxidase.The melatonin-dependent regulated genes were interesting candidates for further physiological and molecular investigations of melatonin improved salt stress tolerance in cucumber.

To facilitate the global analysis of gene expression,the signi?cantly di?erentially expressed genes were assigned to di?erent functional categories using GOseq R package.The annotations were veri?ed manually and integrated using gene ontology (GO)classi?cation (http://www.gene https://www.wendangku.net/doc/094852731.html,).GO assignments were used to classify the functions of the di?erentially expressed genes.The 217selected signi?cantly di?erentially expressed genes between CK and M500(q -value <0.05)were categorized into 36functional groups that were clustered in three main catego-ries (biological process,cellular component,and molecular function)of the GO classi?cation (Fig.9).Gene ontology (GO)analysis of these genes did not reveal signi?cantly enriched groups of genes belonging to the cellular compo-nent categories.In the two main categories (biological pro-cess and molecular function)of the GO classi?cation,there were 16and 20functional groups,respectively.Oxi-dation/reduction process (GO:0055114)with 43genes was dominant in the main category of biological process.Cata-lytic activity (GO:0003824)consisted of 126genes was dominant in the main categories of molecular function.We also noticed a high percentage of genes from func-tional groups response to stress (GO:0006950)with 26genes and oxidoreductase activity (GO:0016491)with 45genes,in the two main categories,respectively.The genes of the di?erent expression clusters associated with di?erent functional categories clearly indicate the molecular and cellular events involved in cucumber lateral root development.

To identify the biological pathways that are active in melatonin-induced lateral root formation process,we mapped the detected genes to reference canonical path-ways in the Kyoto Encyclopedia of Genes and Genomes (KEGG)(http://www.genome.ad.jp/kegg/)[37]and com-pared these with the whole transcriptome background,with a view of searching genes involved in metabolic

Fig.4.Hierarchical clustering of the di?erentially expressed genes,using the RNA-seq data derived from three samples (CK,M10,and M500)based on log 10RPKM values.

Zhang et al.

signal transduction pathways that were signi?cantly enriched.The di?erentially expressed genes between M500and CK were assigned to 57KEGG pathways (shown in Data S2).Those pathways with the greatest representation by unique genes were the metabolic pathways (ko01100)with 21members and biosynthesis of secondary metabo-lites (ko01110)with 13members.These annotations provide a valuable resource for investigating speci?c pro-cesses,functions and pathways involved in melatonin-induced lateral root formation process.

RNA-seq analysis created a comprehensive view of the participation of several multigene families in melatonin’s e?ect on cucumber.Melatonin has a positive e?ect in cucumber lateral root formation process under NaCl stress.Some members of multigene family were identi?ed,and their expression pro?les were characterized in detail.These results indicated that these genes are likely to partic-ipate in the regulation of melatonin-related lateral root formation.

All major processes of life depend on di?erential gene expression,which is largely controlled by the activity of transcription factors.At the process of transcription,a number of di?erent factors recognize DNA in a sequence-speci?c manner to regulate the frequency of initiation of transcription.Transcription factors can be activators,repressors,or both.In the present study,we detected many transcription factors that were di?erentially expressed between the melatonin-treated and untreated

L o g 2 (r a t i o )

1.5

0.0

–0.5–1.0–1.5–2.0–2.5–3.01.00.5

0.0

–0.5L o g 2 (r a t i o )

L o g 2 (r a t i o )

1.0

Subcluster_1, 95 genes

Subcluster_3, 102 genes

Subcluster_4, 59 genes

Subcluster_5, 22 genes

Subcluster_6, 53 genes

Subcluster_2, 32 genes

–1

–2–3–4

0.50.0–0.5–1.0–1.5–2.00.0

0.00.51.01.52.02.53.0

–0.5

–1.0–1.5–2.0–2.5C K

M 10

M 500

C K

M 10

M 500

C K

M 10

M 500

C K

M 10

M 500

C K

M 10

M 500

C K

M 10

M 500

Fig.5.The clustering of di?erentially expressed genes.The six major clusters obtained by K-means algorithm,representing up-regulated (3,5),transient (6),and down-regulated (1,2,4)clusters.Expression ratios are expressed as Log2.

DGE analysis of melatonin on lateral root formation in cucumber

samples under NaCl stress.These transcription factors include MYB,WRKY,NAC,ethylene-responsive tran-scription factor and so on.Most of the transcription fac-tors were down-regulated by melatonin treatment (Table 3).Among all the identi?ed MYB transcription fac-tors,seven showed signi?cantly di?erentially expressed pat-tern.Six of them were down-regulated when treated with melatonin.We detected three signi?cantly di?erentially expressed WRKY transcription factors,two of them were up-regulated by melatonin treatment,and another exhib-ited a lower expression level in melatonin-treated samples.All the di?erentially expressed ethylene-responsive transcription factors and NACs were down-regulated by melatonin (Table 3).Regulation of gene expression at the transcription level controls many crucial biological pro-cesses.Most of these transcription factors have the ability to negatively regulate root related genes and therefore sup-press the root formation process.Seo and Park [6]reported that MYB transcription factor can regulate lateral root meristem activation.The MYB96-de?cient knockout mutant produced more lateral roots and was more suscepti-ble to drought stress.Over-expression of MYB tanscrition factor can reduce the root growth and lateral root density [38–40].The root growth of NAC2-overexpressing seed-lings was also severely inhibited (nearly no root growth)compared with WT [41].Thus,the down-regulation of these transcription factors may promote the growth of lat-eral roots.The WRKY transcription factors modulate the root system in a positive way.WRKY over-expressed plants showed an increase in lateral root number [42].According to our result a majority of WRKY transcription factors showed high expression level when treated with mel-atonin.This is consistent with our conclusion that melato-nin can improve the ability to develop more lateral roots.These observations indicate that transcription factors play an important role in the process of melatonin-induced lat-eral root formation.

Root formation process is accompanied with an exten-sive cell wall formation.Melatonin a?ects the lateral root formation process through cell wall formation.Changes in expression of cell wall –related genes were apparent in the melatonin-treated samples.We detected 12cell wall –related genes signi?cantly di?erentially expressed in the samples of CK and M500.These genes include

transcrip-

Fig.6.Correlation of fold change values from RNA-seq and Q-PCR.The R 2value is

0.97.

Fig.7.Changes in gene expression pro?le between control sam-ples and samples treated with melatonin.CK-samples pretreated with water;M10-samples pre-treated with 10l mol/L melatonin;M500-samples pretreated with 500l mol/L

melatonin.

Fig.8.Venn diagram showing the number of di?erentially expressed genes between every two samples and the number of joint di?erentially expressed genes.

Zhang et al.

Fig.9.Histogram of gene ontology classi?cation.The results are summarized in three main categories:biological process(BP),cellular component(CC),and molecular function(MF).The right Y-axis indicates the number of genes in a category.The left Y-axis indicates the percentage of a speci?c category of genes in that main category.

Table3.Selected genes about transcription factor with altered expression in the three samples

Group Gene ID Gene Annotation RPKM

CK M10M500

MYB Csa1G561370MYB86Transcription factor MYB8626.4756.9470.21 Csa6G046240MY1R1Transcription factor MYB1R1103.8662.3749.77

Csa6G499890MYB4Myb-related protein Myb441.9356.6512.13

Csa1G071840MYB3Transcription factor MYB383.6162.8834.93

Csa7G046120MYB3Transcription factor MYB334.6214.11 3.50

Csa6G538700MYB3Myb-related protein Hv3319.418.700.89

Csa1G024160MYB05Myb-related protein30523.480.25 2.62 WRKY Csa4G051470WRK41Probable WRKY transcription factor4168.4827.1555.11 Csa3G730800WRK33Probable WRKY transcription factor3369.2026.1331.53

Csa3G727990WRK70Probable WRKY transcription factor7019.6544.7252.04 Ethylene-responsive Csa1G423190AIL5AP2-like ethylene-responsive transcription factor AIL567.0863.7231.11 Csa2G177210ERF92Ethylene-responsive transcription factor1B37.228.1910.69

Csa5G165850EF100Ethylene-responsive transcription factor1A52.4410.4515.49

Csa3G135120ERF92Ethylene-responsive transcription factor1B58.1828.8916.57

Csa5G150420ERF61Ethylene-responsive transcription factor ERF06138.1820.988.53 NAC Csa5G606310NAC29NAC domain-containing protein2974.5047.8728.86 Csa3G101810NAC55NAC domain-containing protein55214.4297.8059.64

Csa6G518170NAC8NAC domain-containing protein821.29 1.95 3.25

Csa2G376790NAC98Protein CUP-SHAPED COTYLEDON215.077.89 1.23 Others Csa1G654860MEP1Metalloendoproteinase19.1351.4273.55 Csa4G303700HS17C17.6kDa class I heat shock protein3226.66151.5359.92

Csa6G408800EDL3EID1-like F-box protein379.4935.3928.84

Csa1G039960DOF54Dof zinc?nger protein DOF5.4116.6863.4148.36

Csa1G541390PATL4Patellin-413.4325.3141.56

Csa6G004540GPAT8Probable glycerol-3-phosphate acyltransferase851.7783.8217.19

Csa4G631570RFS1Probable galactinol–sucrose galactosyltransferase134.3052.3968.51

Csa7G434950G3OX3Gibberellin3-beta-dioxygenase322.4157.2067.58

DGE analysis of melatonin on lateral root formation in cucumber

tion factor WRKY33,extensin,basic endochitinase,endochitinase,pectinesterase,xyloglucan endotransglucos-ylase,and some structural constituent of cell wall.Most genes that participated in cell wall biogenesis and organi-zation showed much higher expression level in the melato-nin-treated samples (Table 4),because root formation could not be separated with cell wall formation.

Peroxidase-related genes were identi?ed among the tran-scripts that were di?erentially expressed when treated with melatonin.We detected 12peroxidase-related genes signi?-cantly di?erentially expressed (Table 5).Eight of them were up-regulated in the melatonin-treated samples.We also tested the peroxidase activity in the samples.Melatonin-treated samples showed highly enhanced perox-idase activity (Fig.10).Peroxidases are important enzymes involved in many developmental processes in plants [43,44].Hydroxyl radical (OH ˙)is formed by peroxidases,and it has a loosening e?ect on cell walls and is,therefore,very important for cell elongation [45].Peroxidases are present in all organs and almost all tissues,but they are particu-larly abundant in roots [46].Peroxidases are involved in root growth and di?erentiation [46].Peroxidases control cell elongation through their auxin oxidase activity.Seed-lings lacking AtPrx33transcripts have shorter roots than the wild-type controls,while seedlings over-expressing AtPrx34exhibit signi?cantly longer roots.These modi?ca-tions of root length are accompanied by corresponding changes of cell length [47].

In the traditional view,auxin plays an important role in the rooting process.We were a?ective to the expression pattern of auxin-related genes but found that those genes exhibited minimal expression di?erences (see Data S3).

Table 4.The expression pattern of cell wall related genes

Gene ID Gene name

RPKM

CK M10M500Csa1G154060Structural constituent of cell wall,putative 17.3641.6686.98Csa4G291360Vegetative cell wall protein gp 71.6988.73162.43Csa7G312940Extensin-3

174.30101.5066.10Csa1G538170Putative pectinesterase/pectinesterase inhibitor 2424.3339.1774.29Csa3G175720Pectinesterase 269.9498.34177.02Csa5G167190Pectinesterase 256.55100.66166.64Csa7G405320Pectinesterase 8.5319.1826.93Csa6G084570Unknown

1.17 3.8911.40Csa3G730800Probable WRKY transcription factor 3369.2026.1331.53Csa6G507520Basic endochitinase C

45.5710.2413.00Csa1G202300Probable xyloglucanendotransglucosylase 25.7548.2763.58Csa4G291360Vegetative cell wall protein gp171.6988.73162.43Csa1G534750

Endochitinase PR4

223.97

110.96

205.72

Table 5.The expression pattern of peroxidase

Gene ID Gene name

RPKM

CK M10M500Csa1G051840PER7Peroxidase 7

44.9678.71136.24Csa1G533650PER1Cationic peroxidase 124.82 4.087.18Csa2G421020PER39Peroxidase 39198.08331.79445.44Csa6G504560PER72Peroxidase 72197.60187.9962.60Csa2G235100PER10Peroxidase 1037.4770.0588.65Csa6G134900PER7Peroxidase 733.6576.31127.25Csa4G285760PER15Peroxidase 152573.662467.371204.17Csa4G286310PERN Peroxidase N 51.4675.42104.75Csa5G608020PER11Peroxidase 1146.1454.579.77Csa7G044780PER56Peroxidase 5616.4075.75156.35Csa2G271470PER51Peroxidase 51

7.1517.0529.39Csa6G083540

PER61Probable peroxidase 61

19.18

33.01

64.45

Fig.10.Melatonin has an e?ect on the peroxidase activities in cucumber roots.CK-samples pretreated with water;M10-samples pretreated with 10l mol/L melatonin;M500-samples pretreated with 500l mol/L melatonin.

Zhang et al.

This may mean that melatonin a?ected the root pattern in an auxin-independent manner.This is consistent with the previous study where root developmental changes elicited by melatonin were independent of auxin signaling.In the study in questions melatonin neither activated auxin-inducible gene expression nor induced the degradation of auxin-responsive marker in that research[17].These data may lead to a better understanding of melatonin’s variable e?ect on modulating the root system in the current study. Melatonin modulates root development in an indirect way,mostly by regulating the reactive oxygen species gen-eration and transcription factors.

In conclusion,we report here the?rst dataset that com-prehensively shows the transcriptional changes during the onset of cucumber root branching.RNA-seq data from a collection of cucumber roots with and without melatonin treatment provided new tools for the analysis of melato-nin’s e?ect on plants.RNA deep sequencing generated 16million sequence reads aligned with17,000genes in each sample.Digital gene expression pro?le screened317DEGs between CK and M500,with121up-regulated genes and 196down-regulated genes.They were categorized into36 GO functional groups.The oxidation/reduction process (GO:0055114)and catalytic activity(GO:0003824)were dominant in the main categories of biological processes and molecular function,respectively.Among those di?eren-tially expressed genes,77of them were regulated in a mela-tonin-dependent manner.The one with the highest di?erentially expressed ratio is a gene encodes peroxidase. This study produced abundant data for the analysis of mel-atonin-induced cucumber lateral roots formation.We high-light the link between the peroxidase and the melatonin-induced lateral root formation.In addition to the peroxi-dase,several transcription factors exhibited synergistic reg-ulatory e?ects with melatonin on lateral root formation. We conclude that melatonin modulates root development in an indirect way,mostly by regulating the reactive oxygen species system.To the best of our knowledge,the present study is the?rst to attempt to perform a deep sequencing of the melatonin’s e?ect in cucumber,which may facilitate identi?cation of melatonin’s e?ect on plants. Acknowledgements

This study was partly supported by the grants to Y.-D. Guo(2009CB119000and2012CB113900).

Author contributions

N.Zhang,H.-J.Zhang,S.Ren and Y.-D.Guo designed research;N.Zhang,H.-J.Zhang,Q.-Q.Sun.B.Zhao and Y.-Y.Cao performed research;N.Zhang,H.-J.Zhang, R.Li,X.-X.Wu,L.Li,S.Ren and Y.-D.Guo ana-lyzed data;N.Zhang,H.-J.Zhang,S.Weeda,S.Ren, R.J.Reiter and Y.-D.Guo wrote the paper. References

1.H UANG S,L I R,Z HANG Z et al.The genome of the cucumber,

Cucumis sativus L.Nat Genet2009;41:1275–1281.

2.L I Z,Z HANG Z,Y AN P et al.RNA-Seq improves annotation

of protein-coding genes in the cucumber genome.BMC

Genomics2011;12:540.

3.P ERET B,L ARRIEU A,B ENNETT https://www.wendangku.net/doc/094852731.html,teral root emergence:a

di?cult birth.J Exp Bot2009;60:3637–3643.

4.M ALAMY JE,B ENFEY https://www.wendangku.net/doc/094852731.html,anization and cell di?erentia-

tion in lateral roots of Arabidopsis thaliana.Development

1997;124:33–44.

5.H IMANEN K,B OUCHERON E,V ANNESTE S et al.Auxin-

mediated cell cycle activation during early lateral root

initiation.Plant Cell2002;14:2339–2351.

6.S EO PJ,P ARK C-M.Auxin homeostasis during lateral root

development under drought condition.Plant Signal Behav

2009;4:1002–1004.

7.M URCH SJ,C AMPBELL SSB,S AXENA PK.The role of serotonin

and melatonin in plant morphogenesis:regulation of auxin-

induced root organogenesis in in vitro cultured explants of St.

John’s wort(Hypericum perforatum L.).In Vitro Cell Dev

Biol2001;37:786–793.

8.C HEN Q,Q I W-B,R EITER RJ et al.Exogenously applied mela-

tonin stimulates root growth and raises endogenous indole-

acetic acid in roots of etiolated seedlings of Brassica juncea.

J Plant Physiol2009;166:324–328.

9.S ARROPOULOU VN,T HERIOS IN,D IMASSI-T HERIOU KN.Mela-

tonin promotes adventitious root regeneration in in vitro

shoot tip explants of the commercial sweet cherry rootstocks

CAB-6P(Prunus cerasus L.),Gisela6(P.cerasus9P.canes-

cens),and MxM60(P.avium9P.mahaleb).J Pineal Res

2012;52:38–46.

10.A RNAO M,H ERNANDEZ-R UIZ J.Melatonin promotes adventi-

tious-and lateral root regeneration in etiolated hypocotyls of

Lupinus albus L.J Pineal Res2007;42:147–152.

11.P AREDES SD,K ORKMAZ A,M ANCHESTER LC et al.Phytomela-

tonin:a review.J Exp Bot2009;60:57–69.

12.M ANCHESTER LC,T AN DX,R EITER RJ et al.High levels of

melatonin in the seeds of edible plants–Possible function in

germ tissue protection.Life Sci2000;67:3023–3029.

13.P OSMYK M,K URAN H,M ARCINIAK K et al.Presowing seed

treatment with melatonin protects red cabbage seedlings

against toxic copper ion concentrations.J Pineal Res2008;

45:24–31.

14.P OSMYK M,B ALABUSTA M,W IECZOREK M et al.Melatonin

applied to cucumber(Cucumis sativus L.)seeds improves ger-

mination during chilling stress.J Pineal Res2009;46:214–223.

15.Z HANG N,Z HAO B,Z HANG H-J et al.Melatonin promotes

water-stress tolerance,lateral root formation,and seed germi-

nation in cucumber(Cucumis sativus L.).J Pineal Res2013;

54:15–23.

16.S ARROPOULOU V,D IMASSI-T HERIOU K,T HERIOS I et al.Mela-

tonin enhances root regeneration,photosynthetic pigments,

biomass,total carbohydrates and proline content in the

cherry rootstock PHL-C(Prunus avium9Prunus cerasus).

Plant Physiol Biochem2012;61:162–168.

17.P ELAGIO-F LORES R,M UNOZ-P ARRA E,O RTIZ-C ASTRO R et al.

Melatonin regulates Arabidopsis root system architecture

likely acting independently of auxin signaling.J Pineal Res

2012;53:279–288.

18.P ERET B,D E R YBEL B,C ASIMIRO I et al.Arabidopsis lateral

root development:an emerging story.Trends Plant Sci2009;

14:399–408.

19.P ARK S,B ACK K.Melatonin promotes seminal root elonga-

tion and root growth in transgenic rice after germination.

J Pineal Res2012;53:385–389.

49 DGE analysis of melatonin on lateral root formation in cucumber

20.B YEON Y,P ARK S,K IM YS et al.Microarray analysis of genes

di?erentially expressed in melatonin-rich transgenic rice expressing a sheep serotonin N-acetyltransferase.J Pineal Res 2013;55:357–363.

21.A LLAKHVERDIEV SI,S AKAMOTO A,N ISHIYAMA Y et al.Ionic

and osmotic e?ects of NaCl-induced inactivation of photosys-tems I and II in Synechococcus sp.Plant Physiol2000;

123:1047–1056.

22.A RNAO M,H ERNANDEZ-R UIZ J.Chemical stress by di?erent

agents a?ects the melatonin content of barley roots.J Pineal Res2009;46:295–299.

23.T HOMAS P,S CHIEFELBEIN https://www.wendangku.net/doc/094852731.html,parison of gene expression

under in vitro and ex vitro conditions during rooting of grape cuttings through mRNA di?erential display.Plant Cell Tiss Organ Cult2005;80:105–109.

24.Y AO YY,N I ZF,Z HANG YH et al.Identi?cation of di?eren-

tially expressed genes in leaf and root between wheat hybrid and its parental inbreds using PCR-based cDNA subtraction.

Plant Mol Biol2005;58:367–384.

25.Q I X,X U Q,S HEN L et al.Identi?cation of di?erentially

expressed genes between powdery mildew resistant near-isogenic line and susceptible line of cucumber by suppression subtractive hybridization.Sci Hortic2010;126:27–32.

26.N IELSEN KL.Serial Analysis of Gene Expression(SAGE):

Methods and Protocols,Vol.387.Springer,Secaucus,NJ, USA,2008.

27.S AMOLSKI I,DE L UIS A,V IZCAINO JA et al.Gene expression

analysis of the biocontrol fungus Trichoderma harzianum in the presence of tomato plants,chitin,or glucose using a high-density oligonucleotide microarray.BMC Microbiol2009;

9:217.

28.W ILHELM BT,L ANDRY J-R.RNA-Seq-quantitative measure-

ment of expression through massively parallel RNA-sequencing.

Methods2009;48:249–257.

29.C LOONAN N,F ORREST ARR,K OLLE G et al.Stem cell tran-

scriptome pro?ling via massive-scale mRNA sequencing.Nat Methods2008;5:613–619.

30.G ARBER M,G RABHERR MG,G UTTMAN M et https://www.wendangku.net/doc/094852731.html,puta-

tional methods for transcriptome annotation and quanti?ca-tion using RNA-seq.Nat Methods2011;8:469–477.

31.D E T OMASI JA.Improving the technic of the Feulgen stain.

Stain Technol1936;11:137–144.

32.S CEBBA F,S EBASTIANI L,V ITAGLIANO C.Activities of antioxi-

dant enzymes during senescence of Prunus armeniaca leaves.

Biol Plant2001;44:41–46.

33.L IVAK KJ,S CHMITTGEN TD.Analysis of relative gene expres-

sion data using real-time quantitative PCR and the2-△△Ct method.Methods2001;25:402–408.

34.M ORTAZAVI A,W ILLIAMS BA,M CCUE K et al.Mapping and

quantifying mammalian transcriptomes by RNA-Seq.Nat Methods2008;5:621–628.

35.B ENJAMINI Y,H OCHBERG Y.Controlling the false discovery

rate-a practical and powerful approach to multiple testing.R Stat Soc Series B-Methodol1995;57:289–300.36.H ARTIGAN JA,W ONG MA.A K-means clustering algorithm.

Appl Stat1979;28:100–108.

37.K ANEHISA M,G OTO S,K AWASHIMA S et al.The KEGG

resource for deciphering the genome.Nucleic Acids Res2004;

32:D277–D280.

38.B OMAL C,B EDON F,C ARON S et al.Involvement of Pinus tae-

da MYB1and MYB8in phenylpropanoid metabolism and secondary cell wall biogenesis:a comparative in planta analy-sis.J Exp Bot2008;59:3925–3939.

39.S HIN R,B URCH AY,H UPPERT KA et al.The Arabidopsis

transcription factor MYB77modulates auxin signal transduc-tion.Plant Cell2007;19:2440–2453.

40.V OLPE V,D ELL’A GLIO E,G IOVANNETTI M et al.An AM-

induced,MYB-family gene of Lotus japonicus(LjMAMI) a?ects root growth in an AM-independent manner.Plant J 2013;73:442–455.

41.H U H,Y OU J,F ANG Y et al.Characterization of transcription

factor gene SNAC2conferring cold and salt tolerance in rice.

Plant Mol Biol2008;67:169–181.

42.Z HOU Q-Y,T IAN A-G,Z OU H-F et al.Soybean WRKY-type

transcription factor genes,GmWRKY13,GmWRKY21,and GmWRKY54,confer di?erential tolerance to abiotic stresses in transgenic Arabidopsis plants.Plant Biotechnol J2008;

6:486–503.

43.G APPER C,D OLAN L.Control of plant development by reac-

tive oxygen species.Plant Physiol2006;141:341–345.

44.K WAK JM,N GUYEN V,S CHROEDER JI.The role of reactive

oxygen species in hormonal responses.Plant Physiol2006;

141:323–329.

45.L ISZKAY A,VAN DER Z ALM E,S CHOPFER P.Production of reac-

tive oxygen intermediates(O2á,H2O2,andáOH)by maize roots and their role in wall loosening and elongation growth.

Plant Physiol2004;136:3114–3123.

46.D UNAND C,C REVECOEUR M,P ENEL C.Distribution of super-

oxide and hydrogen peroxide in Arabidopsis root and their in?uence on root development:possible interaction with per-oxidases.New Phytol2007;174:332–341.

47.P ASSARDI F,T OGNOLLI M,D E M EYER M et al.Two cell wall

associated peroxidases from Arabidopsis in?uence root elon-gation.Planta2006;223:965–974.

Supporting Information

Additional Supporting Information may be found in the online version of this article:

Figure S1.The quality of the raw reads.

Table S1.The primers used for qRT-PCR.

Table S2.Summary of RNA-seq data.

Data S1.List of genes expressed in each category.

Data S2.List of KEGG pathways of the di?erentially expressed genes.

Data S3.List of auxin-related genes expression pattern.

Zhang et al.