Revealing of the microRNA involved regulatory gene networks on

Subscriber access provided by Bibliothèque de l'Université Paris-Sud

Omics Technologies Applied to Agriculture and Food

Revealing of the microRNA involved regulatory gene networks on terpenoid biosynthesis in Camellia sinensis in different growing time points Shiqi Zhao, Xuewen Wang, Xiaomei Yan, Lingxiao Guo, Xiaozeng Mi, Qingshan Xu, Junyan Zhu, Ailin Wu, Linlin Liu, and Chaoling Wei J. Agric. Food Chem., Just Accepted Manuscript • DOI: 10.1021/acs.jafc.8b05345 • Publication Date (Web): 06 Nov 2018 Downloaded from http://pubs.acs.org on November 8, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 38

Journal of Agricultural and Food Chemistry

1 2

Revealing of the microRNA involved regulatory gene networks on

3

terpenoid biosynthesis in Camellia sinensis in different growing time points

4 5

Shiqi Zhao1+, Xuewen Wang1,2+, Xiaomei Yan1, Lingxiao Guo1, Xiaozeng Mi1, Qingshan Xu1,

6

Junyan Zhu1, Ailin Wu1, Linlin Liu1, Chaoling Wei1*

7 8

*Corresponding author; E-mail: [email protected]; Tel/fax: +86 551 65786765

9

1State

Key Laboratory of Tea Plant Biology and Utilization, Anhui Agricultural University,

10

West 130 Changjiang Road, Hefei, 230036, Anhui, People’s Republic of China

11

2Department

12

+

13

Email of authors: Shiqi Zhao, [email protected], Xuewen Wang, [email protected],

14

Xiaomei Yan, [email protected], Lingxiao Guo, [email protected], Xiaozeng

15

Mi,

16

[email protected], Ailin Wu, [email protected], Linlin Liu, [email protected]

of Genetics, University of Georgia, Athens, USA

These authors contributed equally to this work.

[email protected],

Qingshan

Xu,

[email protected],

17 18 19 20 21 22

ACS Paragon Plus Environment

Junyan

Zhu,


23

Abstract

24

Tea, made from leaves of Camellia sinensis, has long been consumed worldwide for its

25

unique taste and aroma. Terpenoids play important roles not only in tea beverage aroma

26

formation, but also in the productivity and quality of tea plantation due to their significant

27

contribution to light harvesting pigments and phytohormones. To date, however, the

28

regulation of terpenoid synthase genes remains unclear. Herein, the analyses of metabolomics,

29

sRNAs, degradome and transcriptomics were performed and integrated for identifying key

30

regulatory miRNA-target circuits on terpenoid biosynthesis in leaf tissues over five different

31

months in which the amount of terpenoids in tea leaves varies greatly. Four classes of

32

miRNA-TF pairs which might play a central role in the regulation of terpenoid biosynthesis

33

were also uncovered. Ultimately, a hypothetical model was proposed that mature miRNAs

34

maintained by light regulator at both the transcriptional and posttranscriptional levels

35

negatively regulate the targets to control terpenoid biosynthesis.

36

Keywords: Tea plant, miRNA-target pair, co-expression network, terpenoid volatiles, light

37 38 39 40 41 42 43 44


Page 2 of 38

Page 3 of 38

45


Introduction

46

Tea is the most widely consumed, refreshing, non-alcoholic and thirst-quenching beverage

47

made from the leaves of Camellia sinensis. It has a variety of positive health benefits1. The

48

quality of tea is important to its market value and is mostly depended on its taste and aroma.

49

In general, caffeine, tea polyphenols, amino acids and tea polysaccharides play important

50

roles in tea taste, while the volatile compounds generate the aroma2-4. Due to seasonal climate

51

changes and nutritional status variations of tea plants, the quality of “spring tea”, “summer

52

tea”, and “autumn tea” varies in China5, 6. The abundances of terpenoid volatiles such as

53

monoterpenes and sesquiterpenes, key aroma compounds for tea products’ quality7, varied

54

over growing seasons.

55

In model plants, the terpene biosynthesis pathway has been extensively studied. The MEP

56

pathway in plastids is mainly responsible for producing monoterpenes and diterpenes, while

57

the MVA pathway is largely associated with sesquiterpenes and triterpenes production8. TPS

58

of different classes catalyze the rate-limiting step of converting terpenoid precursors into

59

monoterpenes, sesquiterpenes, and diterpenes, respectively. Many genes (for instance, DXR,

60

HMGR, and TPS) in related to terpenoid biosynthesis in tea trees have been studied, and

61

some genes has been functionally characterized in heterologous model plants or in vitro, as a

62

genetic transformation system for tea plants has not yet been established6, 9-11. At present, the

63

majority of the studies have focused on volatile terpene synthases in tea plants, including

64

monoterpene synthases and sesquiterpene synthases. Several identified terpene synthase genes

65

such as those encoding linalool synthase, nerolidol synthase, and ocimene synthase6, 12-15 have

66

been expressed in Escherichia coli or yeast expression systems and functionally characterized



67

in vitro. However, a few genes, such as (E)-nerolidol/linalool synthase has been further shown

68

to catalyze the expected biochemical reactions in vivo14.

69

MiRNAs play predominant roles in regulating the expression of enzyme genes involved in

70

various physiological processes responding to various stresses (such as light, temperature,

71

diseases and insect pests). Several studies have shown that terpenoid volatile production is

72

tightly regulated through TFs by miRNAs, while direct regulation of enzyme genes in

73

terpenoid biosynthesis pathway through miRNA has not been proven yet. In Pinus taeda,

74

PtMYB14, positively regulates the production of multiple sesquiterpenes, while MsMYB

75

negatively regulates the accumulation of monoterpene and perturbs the sesquiterpene and

76

diterpene derivatives in Mentha spicata16,

77

terpenoid volatiles, probably not solely by activating the transcription of TPS genes, but also

78

probably by enhancing the metabolic flux towards the isoprenoid pathway18. A regulatory

79

network formed by auxin-responsive factors (ARF6/8) and MYBs (MYB21/24) is involved in

80

regulating floral production of volatile sesquiterpenes in the synchronized flower

81

development process19,

82

sesquiterpene production through integrating both GA and JA signals into the transcriptional

83

regulation of sesquiterpene synthase genes21. Similarly, in three JA-responsive AP2/ERF TFs,

84

ORCA3 increases the accumulation of terpenoid indole alkaloids; AaERF1/2, are involved in

85

positively regulating expression of sesquiterpene synthase genes22-24. In Gossypium hirsutum,

86

GaWRKY1 activates the transcription of the sesquiterpene synthase gene encoding CAD1,

87

consequently leading to the biosynthesis of sesquiterpene gossypol25. A spearmint YABBY

88

TF, MsYABBY5, negatively regulates monoterpenes production26. Moreover, limited reports

20.

17.

The Arabidopsis MYB TF, PAP1, increases

In Arabidopsis thaliana inflorescences, MYC2 promotes


Page 4 of 38

Page 5 of 38


89

also revealed that miRNA is also involved in the regulation of terpenoid volatile production in

90

plants. In Arabidopsis thaliana, miR156 negatively regulates (E)--caryophyllene

91

biosynthesis through its target SPL9, which directly binds to TPS21 promoter and activates its

92

expression27. Further investigations are required for a better understanding the regulatory

93

roles of miRNAs in terpenoid volatile production in plants. In recent years, although a large

94

number of genes including those encoding miRNAs have been obtained from tea trees by

95

cloning or RNA sequencing, few studies about miRNAs and TFs on the regulation of terpenes

96

in tea plants have been reported.

97

The purpose of this study is to comprehensively explore the function of miRNA-target

98

pairs in regulating terpenoid biosynthesis through identifying the potential month-responsive

99

miRNAs and their targets in the tea plant in different growing times (April, June, August,

100

September, and October). The data of metabolomics, sRNAs, degradome and transcriptomics

101

in tea leaves excised in different time points were generated and integrated to identify the key

102

regulatory miRNA-targeted circuits of terpenoid biosynthesis. Our results revealed the

103

regulatory gene networks containing miRNA-target pairs on terpenoid biosynthesis pathway

104

in tea.

105

Materials and methods

106

Plant Materials and Total RNA Extraction

107

8-year-old cloned tea plants (Camellia sinensis cv. Shuchazao) were grown in the

108

experimental nursery under natural daylight conditions at the 916 Tea Plantation in Shucheng

109

County. The cultivation and management were described previously6. Twenty tea plants were

110

grown at 33 cm×120 cm in an experimental plot (4 rows×5 columns). In total, three



111

experimental plots were setup in this study. The second fully expanded leaves (position with

112

reference to the apical bud) excised from tea shoots in an experimental plot were pooled and

113

used as an independent biological sample in a sampling month. Three biological samples were

114

obtained from all the plots for each sampling month (April 10, June 6, August 5, September

115

15, and October 18 in 2015). All the samples were immediately frozen in liquid nitrogen and

116

then maintained at -80 °C.

117

RNA was extracted using the Total RNA Purification kit (NorgenBiotek Corporation,

118

Canada) according to the manufacturer’s protocol. The quality and quantity of the RNA

119

extracts were examined using Agilent Technologies 2100 Bioanalyzer (Agilent Technologies,

120

Palo Alto, CA, USA). The RNA of tea plants in response to shading were from our previously

121

study28.

122

Small RNA Library Construction, Sequencing and Sequencing Data Analysis

123

One microgram of total RNA for each sample was used for small RNA library construction

124

with TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA) according to

125

manufacturer’s instructions. The libraries were used to for single-end sequencing on an

126

Illumina Hiseq2500 at the LC-BIO (Hangzhou, China) following the vendor’s recommended

127

protocol.

128

MiRNA sequences were identified by ACGT101-miR v4.2. Firstly, adapter dimers, junk,

129

low complexity, mRNA, RFam RNAs (rRNA, tRNA, snRNA, snoRNA) and repeats were

130

removed from the raw reads. Subsequently, unique sequences with length in 18- ~ 25-nt were

131

mapped to the precursors in miRBase 21.0 by Bowtie search to find those miRNAs, which are

132

known. Length variation at both 3’ and 5’ ends and one mismatch inside of the sequence were


Page 6 of 38

Page 7 of 38


133

allowed in the alignment. The mapped pre-miRNAs were further aligned against tea genome

134

(www.plantkingdomgdb.com/tea_tree/) to detect their genomic locations by Megablast. Then,

135

the unmapped unique sequences were also searched against tea genome by Bowtie, and the

136

hairpin RNA structure-containing sequences were predicated from the flank 120 nt sequences

137

using RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi). The criteria for

138

secondary structure prediction were applied according to Ambady and Dominko29.

139

Differentially expressed miRNAs between different time points were identified using

140

Student’s t test method based on normalized read counts. In this test, differences in transcript

141

level were considered statistically significant at p-value  0.0530, 31.

142

Degradome Library Construction, Sequencing, and Sequencing Data Analysis

143

Approximately 20 g of total RNA was used to construct the degradome library. The

144

method differed considerably from previous reports32-34 with the following modifications: (1)

145

Approximately 150 ng of poly (A)+ RNA was used as input RNA and annealing with

146

Biotinylated Random Primers; (2) Strapavidin capture of RNA fragments through

147

Biotinylated Random Primers; (3) 5’ adaptor ligation to those RNAs containing 5’-

148

monophosphates only; (4) Reverse transcription and PCR for amplification; (5) Libraries

149

were sequenced using the 5’ adapter only, to obtain the sequences of the first 36 nucleotides

150

of the inserts which presented at the 5’ ends of the original RNAs. Then the single-end

151

sequencing (50 bp) was performed on an Illumina Hiseq2500 at LC-BIO (Hangzhou, China)

152

following the vendor’s recommended protocol.

153

The CleaveLand 3.0 software package33, 35 and the ACGT101-DEG program (LC Sciences,

154

Houston, TX, USA) were used to identify candidate targets of the miRNAs. Briefly, the



Page 8 of 38

155

degradome reads were mapped to the appropriate transcriptome, followed by integrating the

156

mapped

157

'CleaveLand3_map2dd.pl'. Small RNA/mRNA alignments were then generated using

158

'targetfinder.pl'. The degradome density file was compared to the target predictions and

159

significant hits were presented using script 'CleaveLand3_analysis.pl'. "T—plots" of the

160

targets were generated using the 'CleaveLand3_t—plotter.pl' script. All targets were grouped

161

into five categories based on the abundance of the resulting mRNA tag relative to the overall

162

profile of degradome reads that matched the target36. Small RNA and degradome datasets

163

were all deposited in the NCBI SRA database (BioProject accession number: PRJNA429080).

164

The GO enrichment analysis was performed by Agrigo (http://bioinfo.cau.edu.cn/agriGO/)

165

and the non-redundant GO enrichment terms were obtained and visualized using Revigo

166

(http://revigo.irb.hr/). The KEGG enrichment analysis was performed by Perl script.

167

Differentially Expressed Target Gene Analysis

degradome

data

into

a

"degradome

density

file"

using

script

168

Transcriptome sequencing reads from tea samples collected over the five different months

169

were generated from our recent study6. Salmon v0.8.2 was used to map all the clean reads to

170

the CDS of all genes and to count the number of reads mapped to CDS of each gene with

171

default parameters. The RPKM of each gene’s CDS were calculated by Perl script based on

172

the length of the gene’s CDS and the number of reads mapped to this gene’s CDS.

173

Differential expression analysis of samples was performed using the DESeq2 package (1.14.1)

174

in R based on clean read counts. The cutoff for pairwise comparisons was set at a fold change

175

greater than 2 and q-values (FDR adjusted p -values) < 0.01.

176

Co-expression analysis


Page 9 of 38


177

The weighted co-expression network analysis of 11,245 DEGs was constructed using the

178

WGCNA (v1.51) package in R (v3.3.2)37. Unsigned co-expression networks were constructed

179

using the blockwise module function with default settings, except that the power was 23 with

180

“unsigned” of TOMType while minModuleSize and mergeCutHeight were set at 30 and 0.2,

181

respectively. To further investigate the relationship between each module and terpenoid

182

volatiles, the eigengene value was calculated for each module and used to test the association

183

with 11 kinds of terpenoid volatiles from our recent report6. The total connectivity and

184

intramodular connectivity (function softConnectivity), kME (for modular membership, also

185

known as eigengene-based connectivity), and kME-P value were calculated for the DEGs,

186

which were clustered into 14 modules. The gene co-expression networks for selected genes

187

were visualized using Cytoscape (version 3.4.0)38.

188

Data validation and gene expression analysis by qRT-PCR

189

To validate the high-throughput sequencing results, qRT-PCR was performed to quantify

190

the transcript levels of nine differentially expressed miRNAs and nine targets. The RNA

191

samples used for qRT-PCR were the same as for the experiments mentioned above. The

192

reverse transcription reactions for miRNAs were carried out using the miRcute miRNA

193

First-Strand cDNA Synthesis Kit (Code No. KR211, TIANGEN, Beijing, China) for tailing

194

qRT-PCR and the PrimeScript RT Reagent Kit (Code No. RR037A, TaKaRa, Dalian, China)

195

for stem-loop qRT-PCR. The reverse transcription reactions for target cDNA synthesis were

196

conducted using the PrimeScript RT Reagent Kit (Code No. RR036A, TaKaRa, Japan). The

197

specific miRNA forward primers were designed based on the sequences of each miRNAs and

198

the primers for the target genes were designed using the Primer Premier 5 software



199

(http://www.premierbiosoft.com/primerdesign/) (Table S1). GAPDH (GenBank: GE651107)

200

and CsmiR222 (UUUCCAAGACCACCCAUGCCGA) were used as reference genes for

201

target transcript and miRNA quantification respectively. The qRT-PCR reactions for miRNAs

202

and their targets were performed in 96-well plates using a miRcute miRNA qPCR Detection

203

Kit (Code No. FP411, TIANGEN, Beijing, China) and SYBR Premix Taq™ II (Code No.

204

RR820A, TaKaRa, Dalian, China) on a CFX96™ real-time PCR system (Bio-Rad, USA). The

205

relative fold change in gene expression was calculated using the 2−ΔCt method39. Each PCR

206

reaction was repeated three times independently.

207

Modified 5’RNA ligase-mediated RACE for the mapping of mRNA cleavage sites

208

In order to validate the cleavage sites of miRNA in target genes identified by degradome, a

209

modified RLM-5’ RACE was performed by using FirstChoice RLM-RACE Kit (Invitrogen,

210

Thermo Fisher Scientific) as previously described40. The primers used to amplify cleavage

211

products of tea miRNA target genes through 5’RLM-RACE are listed in (Table S1).

212

Northern blot of miRNAs

213

Northern blot analysis was performed as previously reported41 with the DIG-High Prime

214

DNA Labeling and Detection Starter Kit II (cat 11745832910, Roche, US) following the

215

manufacturer’s instructions. Total RNA was separated in a 14% polyacrylamide gel

216

containing 7 M urea and electrically transferred to Hybond-N+ membranes (GE Amersham).

217

The 5’ end-modified digoxin-labeled miRNA probes were synthesized by Sangon Biotech

218

(shanghai, China).

219

Results

220

Small RNA sequencing profile


Page 10 of 38

Page 11 of 38


221

To study the miRNA-involved posttranscriptional regulation of aroma compound

222

metabolism, especially terpene in the tea tree during different periods of the growing months,

223

miRNAs and their transcript levels were investigated by small RNA (sRNA) sequencing.

224

Fifteen small RNA libraries established over five different time points were sequenced by

225

Illumina technology. After a series of data filtering, the valid sequences of 18- to 25- nt were

226

obtained (Table S2). The size distributions of the total and unique valid sRNAs were

227

summarized in Figure S1 and showed similar patterns in the five months. The 24-nt RNAs

228

were the most abundant and diverse sRNAs in the fifteen libraries (Figure S1).

229

The

valid

sequences

were

aligned

against

the

tea

genome

230

(www.plantkingdomgdb.com/tea_tree/) and the precursors in miRBase 21.0. The sequences

231

mapped to the miRBase were identified as known miRNAs and those mapped to tea genome

232

only were considered as predicted miRNAs. All detected miRNAs in five months were

233

categorized into five different groups (Figure 1 and Figure 2). In total, 857 mature miRNAs

234

originated from 687 pre-miRNAs and corresponding to 710 unique mature miRNAs were

235

identified. Of these pre-miRNAs and unique mature miRNAs, 420 unique mature miRNAs

236

corresponding to 425 pre-miRNAs were grouped into Group 1, -2 and -3 (gp1, -2, and -3), all

237

sharing a high similarity with 281 known plant miRNAs in miRBase 21.0. In addition, 290

238

unique mature miRNAs, corresponding to 262 pre-miRNAs in Group 4 (gp4), were identified

239

as novel candidate miRNAs, which were not registered in miRBase 21.0. Mature miRNAs

240

with 21-nt occupied 39.67 % of the total unique miRNAs, followed by 25.67% for those with

241

24-nt and 11.32% for 22-nt. The similar pattern could be found in each sampling time point

242

(Figure 1A). Of these mature miRNAs, 58 (38 known miRNAs and 20 novel miRNAs) were



243

expressed only in September; 36 (13 known miRNAs and 23 novel miRNAs) were expressed

244

only in April. While no miRNA expressed specifically in June was found. A total of 454

245

miRNAs, consisting of 165 known miRNAs and 289 novel miRNAs, expressed in all five

246

sampling time points (Figure 1B).

247

Furthermore, the results of first nucleotide bias analysis showed distinct patterns between

248

known and novel miRNAs. In the known miRNAs, uracil (U, 43.18%) was the most common

249

nucleotide at the 5’ end, followed by guanine (G, 23.30%) and adenosine (A, 20.27%),

250

whereas in the novel miRNAs, adenosine was the most common nucleotide (49.24%) at the 5’

251

terminal, followed by U (33.74%) (Table S3).

252

Differential expression of miRNAs in tea tree in different growing months

253

To further identify the growing month-specific miRNAs in the tea tree, the DEmiRNAs in

254

fifteen libraries was analyzed and compared based on the normalized read counts. In total,

255

221 mature miRNAs (P-value  0.05) showed differential expression patterns of miRNAs

256

including 134 known and 87 novel miRNAs, were found at least one pairwise comparison of

257

samples. The heatmap of potential DEmiRNAs is illustrated in Figure 3. According to the

258

clustering analysis, the expression of 221 miRNAs was very different across the five months.

259

In general, miRNAs with high or low expression levels had similar distribution patterns

260

(Figure 3). Over the five time points, the number of miRNAs with high expression levels (79,

261

35.7% of the total miRNAs) in April was more than those in the other time points.

262

Interestingly, the number of miRNAs at low expression levels in April (59, 26.7%) was also

263

more than those in other months. However, the numbers of miRNAs with high and low

264

expression levels reduced to 23 (10.4%), 27 (12.2%) in October and 24 (10.9%), 39(17.6%) in


Page 12 of 38

Page 13 of 38

265


June, respectively.

266

Furthermore, to study the special expressed miRNA in the five months, we calculated the

267

fold change of miRNA expression through comparing the normalized counts of the same

268

miRNA at different time points. Six miRNAs were identified with significantly higher

269

expression levels in one sample than all other four samples (more than 4 fold changes)

270

(Figure S2). These miRNAs were referred to as preferentially expressed miRNAs

271

(PEmiRNAs). Among these miRNAs, four PEmiRNAs were found in April that contains two

272

miR319, one miR858 and a novel miRNA. A novel miRNA was found each in August and

273

September. But none of PEmiRNAs was found in Jun or October.

274

Target prediction of the known and novel miRNAs by degradome sequencing

275

For better understanding the potential functions of the miRNAs identified, we applied the

276

degradome sequencing technology to identify the targets of the identified miRNAs.

277

Degradome sequencing resulted in 45,067,616 raw reads, representing 11,456,050 unique

278

reads. After a series of data filtering and mapping, 27,604,585 (61.25% of 44,725,436 clean

279

reads) clean reads were successfully mapped to the 32,753 unigenes (88.64% of the 36 951

280

input cDNA sequences) of tea plant. All of the potential cleaved target unigenes were

281

classified into five categories, based on the signature abundance at each occupied unigenes

282

position36. There were 209, 14, 530, 70, and 316 targets in categories 0, 1, 2, 3, and 4

283

respectively (Figure S3).

284

A total of 646 targets of 293 miRNAs were identified from the mixed degradome,

285

containing 101 transcripts targeted by 52 novel miRNAs. Majority of the miRNAs (217)

286

likely cleaved a few transcript targets (the number < 5), while 75 miRNAs were possibly to



287

cleave 5 or more different transcript targets. 54 targets were identified to be cleaved by both

288

ppe-MIR169i-p3_2ss1CT19GT and ppe-MIR169i-p5_2ss1CT19GT, respectively. This was

289

the maximal number of transcripts cleaved by the same miRNA. While 11 miRNAs were

290

detected to cleave only one transcript target. Cleavage sites of some of these target-miRNA

291

pairs have been verified by 5’ RLM-RACE (Figure 4).

292

GO and KEGG pathway analyses of targets

293

Analyses of GO and KEGG on miRNA target genes were performed to further investigate

294

the potential biological roles of miRNAs in tea plants. Of all 647 target genes, 410 target

295

genes were classified into three GO categories: biological process, molecular function and

296

cellular component (Figure S4). The biological processes with the highest numbers of targets

297

included “metabolic process”, “cellular process”, “biological regulation”, “localization” and

298

“response to stimulus”. The GO terms of “cell”, “cell part”, and “organelle” in the cellular

299

component category contained the most abundant targets. The GO terms of “binding” and

300

“catalytic activity” in the category of molecular function had more targets than any other

301

terms. The enrichment analysis results were shown in Figure 5 and Table S4. For the

302

biological process, 30 of 42 enriched GO terms were classified as the metabolic or

303

biosynthetic process including “RNA metabolic process”, “nucleic acid metabolic process”,

304

“organic cyclic compound biosynthetic process” and “cellular biosynthetic process”. Some

305

interesting GO terms, for instance “response to endogenous stimulus”, “response to hormone”,

306

“aromatic compound biosynthetic process”, were revealed. “Cell” and “cell part” were the

307

most enrichment terms in the cellular component category. The two enrichment terms only in

308

molecular function category were both related transcription factor.


Page 14 of 38

Page 15 of 38


309

KEGG annotation was carried out for pathway analysis, and 90 third level pathways out of

310

18 second level pathways, which belonged to five first level pathways, were obtained. In the

311

top 20 pathways ranked by the KEGG-annotated gene number, “plant-pathogen interaction”

312

and “plant hormone signal transduction” included the most target genes, followed by “- acid

313

metabolism”, “tyrosine metabolism” and “RNA transport” (Table S5). In addition,

314

“biosynthesis of secondary metabolites”, including “terpenoid”, “phenylpropanoid”,

315

“flavonoid”, “indole alkaloid”, had more genes. According to enrichment analysis,

316

“-Linolenic acid

317

degradation”, “tyrosine metabolism” were the significant enriched third level pathways (FDR

318

< 0.05).

319

Correlation analysis of miRNAs expression profiles and their targets

metabolism”, “plant hormone signal transduction”, “fatty acid

320

Our transcriptomics study revealed that in total, 99 target genes of 72 DEmiRNAs were

321

also expressed differentially (Table 1, Table S6). A negative expression correlation was found

322

between miRNA and its target in 95 pairs (Figure 6). For example, miR160s (ptc-miR160a,

323

mtr-miR160a) had a lower expression in August and September, while their targets ARFs

324

(CSA010141.1, CSA035909.1) were highly expressed in the two time points.

325

qRT-PCR and miRNA northern blot were employed to confirm the expression profile of

326

miRNA-target pairs including 8 known miRNAs and 2 predicted miRNAs (Figure 7). As is

327

shown in Figure 7, the results of both qRT-PCR and northern blot validated sequencing

328

results. In addition, the expected negative correlation between the expression of the 8 selected

329

miRNAs and their corresponding target transcripts were confirmed by qRT-PCR in Figure 7,

330

suggesting that transcriptional repression of the target genes was most likely resulted from



331

their negative miRNAs regulators. However, we detected an unexpected result for 2

332

miRNA-target pairs that presented a positive relationship (Figure 7).

333

Identification of miRNA-targets related to terpenoid volatiles biosynthesis

334

Based on the miRNA-target pairs and the unrooted neighbor-joining trees, seven miRNAs

335

and seven targets were found to be related to terpenoid biosynthesis in tea tree (Figure S5).

336

Two miR858 (ath-miR858b, ath-miR858b_R-3) were found to target MYB TF

337

(CSA020478.1, CSA024950.1, CSA036360.1). MiR3630 (vvi-miR3630-3p_L-3_1ss22CA)

338

was predicted to target a MYC TF (CSA021679.1) and miR156 (stu-miR156a, mtr-miR156e,

339

rco-miR156f) and miR535 (mdm-miR535a) were likely to target a SPL TF (CSA023442.1).

340

Furthermore, two novel miRNAs (PC-3p-81_33418 and ptc-MIR156f-p3_2ss3CT22CT) were

341

predicted to target two ERF TFs (CSA025745.1 and CSA017147.1), respectively.

342

The expressions of MYBs or SPL positively correlated with the contents of most

343

monoterpenes, negatively related with the contents of most sesquiterpenes. However, ERFs

344

were opposite to MYBs and SPL. Furthermore, the expression of MYC negatively correlated

345

with both contents of monoterpenes and almost all sesquiterpenes. (Table S7 and Table S8)

346

Gene co-expression network related to terpenoid volatiles biosynthesis

347

In order to comprehensively understand the miRNA-target pairs responsible for regulation

348

of terpenoids biosynthesis in different growing months, our previous transcriptomics and

349

metabolomics data sets (Table S7)6 were used for co-expression analysis. We identified 14

350

distinct co-expression modules containing 11,245 unigenes through WGCNA package in R

351

(Figure S6). Of all the modules, the size ranged from 36 (cyan module) to 3,494 genes

352

(turquoise module). To explore the network connections for the regulatory genes related to


Page 16 of 38

Page 17 of 38


353

terpene synthases, we focused on this seven target TF genes identified by the neighbor-joining

354

trees. Finally, five were directly connected with 3591 edges (the edge weight higher than 0.4)

355

in this network, suggesting that these five genes were regulated through these gene networks.

356

We also created a sub-network containing 387 genes within

six categories related to

357

signal transduction, transcription factor, terpenoid biosynthesis, light responsive, heat and

358

cold responsive (Figure 8, Table S9). Each of the six categories contained 144, 83, 25, 70, 39

359

and 26 genes, respectively and some genes were multiple functional and thus grouped into

360

multiple categories (Table S10). The five categories (except terpenoid biosynthesis) may

361

trigger the regulation process for terpenoid volatiles biosynthetic pathway in tea plant through

362

different signal pathways. Interestingly, distinct TF family of hub gene connected different

363

edge genes. In the ERF family, the two ERF genes (CSA025745.1 and CSA017147.1)

364

targeted by different miRNAs were involved in all six gene networks shown in Fig 8. While

365

two MYB genes (CSA024950.1 and CSA036360.1) targeted by same miRNA

366

(ath-miR858b-R_3) shared no edge gene, indicating these two targets were involved in very

367

limited networks.

368

When the edge weight higher than 0.4, 23 genes related terpenoid biosynthesis pathway

369

were found relating to two ERF hub genes, and only two connected with one MYB hub gene

370

(CSA036360.1). Two genes (CSA014975.1 and CSA020011.1), which connected with ERF

371

hub genes, were annotated as monoterpene synthase catalyzing the production of pinene and

372

linalool, respectively. MYB hub gene (CSA036360.1) may regulate two genes (CSA018888.1

373

and CSA008212.1), which encoded a sesquiterpene and a monoterpene synthase, respectively.

374

Only when the edge weight higher than 0.3, two sesquiterpene synthase genes (CSA019490.1



375

and CSA003008.1) connected with MYB (CSA020478.1) and MYC (CSA021679.1), three

376

caryophyllene synthase genes (CSA025571.1, CSA025572.1 and CSA005131.1) connected

377

with SPL (CSA023442.1).

378

Discussion

379

In this study, miRNA and their target identification, classification, differential expression,

380

target determination, and functional correlation in tea plants were studied with an emphasis on

381

the different growing time points and tea character specialized metabolite metabolism. Due to

382

the advance of sequencing technology and recently release of tea genome, our results on

383

miRNA analysis were validated and confirmed with tea genome and gene network data as

384

well as metabolic profiling findings.

385

The identified miRNAs and its’ targets involved in terpenoid biosynthesis

386

Based on small RNA sequencing data, degradome and transcriptome analyses, 99

387

differentially expressed target genes of 72 DEmiRNAs were identified. Importantly, most of

388

these targets belong to TFs, which is consistent with the published reports42. In the previous

389

studies, a large number of miRNAs have been found to target enzyme genes in the terpene

390

synthesis pathway in silico43, three miRNA-enzyme gene pairs in related to terpenoid

391

biosynthesis were identified in our degradome data also. However, none of the miRNAs that

392

target the three genes are DEmiRNAs. So, we focused our attention on miRNA-TF pairs in

393

the DEmiRNA-DEG pairs. Our study revealed that three pairs of miRNA-MYB, one of

394

miRNA-MYC, two of miRNA-ERF, and one of miRNA-SPL were likely involved in

395

terpenoid biosynthesis regulation in tea plants (Figure S5). However, miRNA-WRKY pairs,

396

which regulates the sesquiterpene biosynthesis in cotton25, was not identified in tea plants in


Page 18 of 38

Page 19 of 38


397

this study. Furthermore, in this study, some TFs were found being targeted by conserved or

398

novel miRNAs, which had not been reported before, for instance, miR828-MYB,

399

miR3630-MYC, miR535-SPL, and novel-miRNA-ERF. Further studies are required to

400

validate the regulatory function of these pairs in terpenoid biosynthesis in tea plants.

401

Some environmental factors regulate the synthesis of terpenes via miRNA-TF pairs

402

The environmental factors, such as light, temperature, soil and air humidity, can

403

significantly affect the production and emission of terpenoid volatiles in plants44-46. The

404

relative amount of β–myrcene and α-farnesene significantly decreased with increases in light

405

intensity, whereas nerolidol and linalool showed a significant increase in their relative

406

proportion when the light was on45. The synthesis of monoterpenes (mainly α- and β-pinene)

407

is strongly activated by light47. In our results, the abundance of three sesquiterpenes and one

408

monoterpenes was highest in Jun with the longest duration of sunshine (Table S7). Moreover,

409

the proportions of caryophyllene and nerolidol were highest at 37 °C, while relative amounts

410

of other terpenes were highest at 22°C or 27°C

411

°C), emissions of ocimene remarkably increased46. In addition, low temperature (4 °C) stress

412

could induce the increasing of linalool and decreasing of nerol in tea leaves44. Similarly, we

413

noted that the caryophyllene abundance in tea leaves was highest in Jun at the highest

414

temperature, compared with the other four time points. For some monoterpenes, the highest

415

abundances in tea leaves appeared in April at the lowest temperature.

45.

Interestingly, at high temperatures (> 38

416

Some reports revealed that after light exposure, the expression of most genes in the MVA

417

pathway are suppressed, while almost all genes of the MEP pathway are induced6. The PaIspS

418

mRNA expression in Populus alba is strongly induced by heat stress and light irradiation, but



419

substantially decreases in the dark48. The expression of PtMYB14 which can leads to the

420

accumulation of sesquiterpene decreased slightly responding to cold16. MiR156 and its targets

421

SPL9 which bind TPS21 promoter and activate its expression can be regulated by heat stress27,

422

49.

423

the contents of most terpenoid volatiles (Table S8). In silico analysis revealed that the

424

light-specific and temperature-specific cis-elements were found present in the miRNA

425

promoters. Interestingly enough, enhanced miR156 expression and caryophyllene emission

426

occur simultaneously in Arabidopsis thaliana upon high temperature stress, while SPL9, the

427

target of miR156, can activate the expression of TPS. Here, both miR156 expression and

428

caryophyllene abundance in tea leaves were low in April at the lowest temperature, compared

429

with the other four time points, while the expression of SPL (CSA023442.1) is the highest in

430

April. All of these findings indicate that the abiotic environmental factors may regulate the

431

synthesis of terpenes through miRNA-TF pairs.

432

Light may trigger terpenoid biosynthesis through miRNA-TF pairs

The expression of 7 miRNA-TF pairs identified in this study were significantly related to

433

Despite the fact that light affects terpenoid biosynthesis through some light-responsive

434

genes, it is unclear whether and how miRNAs are involved6, 50. Previously we found that

435

light-specific elements do exist on the promoter regions of many transcription factors6. The

436

present study indicated that light-specific elements are also present in the miRNA promoters

437

(Table S11). Light signaling factors, PIF4, and CDF2 regulate the expression of miRNAs

438

during the dark-to-light transition via two major mechanisms: (1) controlling the transcription

439

through binding directly to the promoters of miRNA genes, (2) affecting the processing of

440

pri-miRNAs51, 52. In this study, the expression of some mature miRNAs and pre-miRNAs in


Page 20 of 38

Page 21 of 38


441

tea leaves were promoted or suppressed under shading treatments (Figure S7). The E3 ligase

442

COP1 is essential for light signaling and stabilizing HYL1 to retain miRNA biogenesis under

443

light53. The COP1–HY5 interaction may specifically target HY5 for ubiquitination and

444

subsequent degradation by the proteasome in the nucleus in the dark54,

445

transferred to the cytoplasm and separated with HY5 under light conditions. Consequently,

446

HY5 directly regulates the downstream genes by binding the promoters of these genes,

447

including MIR genes and chalcone synthase gene56-58. The four TF families (PIF, DOF, COP

448

and HY5) in tea plants were closely associated with the expression of some miRNAs which

449

might regulate terpenoid biosynthesis and the abundance of some terpenoid volatiles (Table

450

S8). Interestingly, these four TF families had the opposite effects on mono- and diterpene

451

biosynthesis. On the other hand, less than half genes in this four families and one miRNA

452

related terpenoid biosynthesis were found closely associated with the duration of sunshine

453

(Table S8). This may be related to the transient expression of the genes or other mechanisms.

55.

COP1 can be

454

A model of regulatory network has been proposed for light signal regulated terpenoid

455

metabolism in tea (Figure 9). Firstly, PIF4 and CDF2 would interact with light photoreceptor

456

(PHY and CRY), and then COP1 is transferred to the cytoplasm from nucleus and separated

457

with HY5 under light. Secondly, PIF4 and CDF2 function at both the transcriptional and

458

post-transcriptional levels to regulate miRNA biogenesis. HY5 binds specifically to the

459

promoter of some MIR genes and positively regulates their expression. Thirdly, transcription

460

factors controlled by miRNA regulate the genes (especially TPS) in terpenoid pathways to

461

control terpenoid volatiles biosynthesis.

462

In this study, analyses of metabolomics, sRNAs, degradome and transcriptomics over



Page 22 of 38

463

different growing time points in tea plants identified a number of miRNA-target pairs. Four

464

classes of miRNA-TF pairs were involved in the terpenoid pathway regulation. A model of

465

regulatory network contained miRNA-target circuits was proposed to dissect the regulation of

466

tea terpenoid metabolism via light signaling and miRNAs. These findings will improve our

467

comprehensive understanding of the regulation of terpenoid biosynthesis in tea and could

468

guide for tea genetic improvement.

469 470

Abbreviations

471

AP2/ERF, Apetala2/ethylene response factor; CAD1, (+)-δ-cadinene synthase; CDF2,

472

cycling DOF transcription factor; CDS, coding sequences; COP1, E3 ligase CONSTITUTIVE

473

PHOTOMORPHOGENESIS 1; CRY, cryptochrome; DEG, differentially expressed gene;

474

DEmiRNA, differentially expressed miRNA; DXR, 1-deoxy-D-xylulose-5-phosphate

475

reductoisomerase; GA, gibberellin; GAPDH, Glyceraldehyde-3-phosphate dehydrogenase;

476

GO, gene ontology; gp, Group; HMGR, 3-hydroxy-3-methylglutaryl coenzyme A reductase;

477

HY5, LONG HYPOCOTYL 5; HYL1, HPONASTIC LEAVES1; JA, jasmonate; KEGG,

478

Kyoto Encyclopedia of Genes and Genomes 2; MEP, 2-C-methyl-D-erythritol 4-phosphate;

479

miRNA, microRNA; MVA, mevalonate; PAP1, Production of Anthocyanin Pigment 1;

480

PEmiRNA,

481

phytochrome-interacting factor;4 RPKM, Reads per Kilobase per Million mapped reads;

482

sRNA, small RNA; TF, Transcription factors; TPS, Terpene synthases;

preferentially

expressed

miRNA;

483 484

Authors and contributors


PHY,

phytochrome;

PIF4,

Page 23 of 38


485

CW conceived and designed the study; SZ analyzed the data. SZ and XW wrote the

486

manuscript; SZ, XY, LG and XM performed the experiments; QX helped to do data analysis,

487

JZ and AW provides experimental guidance. LL did the shading experiment. All authors have

488

read and approved the final version of the manuscript.

489

490

Acknowledgment

491

This work was supported by the National Natural Science Foundation of China (31171608),

492

the Special Innovative Province Construction in Anhui Province (15czs08032), the Special

493

Project for Central Guiding Science and Technology Innovatio of Region in Anhui Province

494

(2016080503B024) and the Programme for Changjiang Scholars and Innovative Research

495

Team in University (IRT1101). We would like to thank the 916 Tea Plantation in Shucheng,

496

Anhui Province, China for providing samples of tea plants. We are grateful to Shengrui Liu at

497

Anhui Agriculture University, who helped us to improve this work.

498

The authors declare that they have no competing interest.

499 500

Supporting Information description

501

Table S1. The primers and adaptors list.

502

Table S2. Overview of sRNA sequencing.

503

Table S3. The distribution nucleotide bias in miRNAs.

504

Table S4. GO enrichment of targets.

505

Table S5. KEGG enrichment of targets.

506

Table S6. The annotation of differentially expressed targets of DEmiRNAs



507

Table S7. The mono- and sesqui-terpene contents (relative abundance) of tea leaves in five

508

months.

509

Table S8. The correlation among genes, terpenes and environment factors.

510

Table S9. Edges in co-expression network.

511

Table S10. Nodes in co-expression network.

512

Table S11. Distribution of elements scanned with promoter analysis software Plant CARE in

513

miRNAs related terpenes biosynthesis.

514 515

Figure S1. Length distribution of total (A) and unique (B) small RNA sequences.

516

Figure S2. The heatmap of preferentially expressed miRNAs.

517

Figure S3. Five categories based on target cleavage positions in degradome.

518

Figure S4. Gene ontology classification of target genes for the identified miRNAs.

519

Figure S5. WGCNA co-expression modules based on DEG data.

520

Figure S6. The miRNAs target related to terpene synthases in tea plant.

521

Figure S7. The qRT-PCR of selected mature miRNAs and pre-miRNAs in tea leaves under

522

shading treatments.

523 524

References

525 526 527 528 529 530 531 532

(1) Gramza-Michałowska, A., Functional Aspects of Tea Camellia sinensis as Traditional Beverage. In Functional Properties of Traditional Foods, Springer: 2016; pp 353-369. (2) Gohain, B.; Borchetia, S.; Bhorali, P.; Agarwal, N.; Bhuyan, L. P.; Rahman, A.; Sakata, K.; Mizutani, M.; Shimizu, B.; Gurusubramaniam, G.; Ravindranath, R.; Kalita, M. C.; Hazarika, M.; Das, S., Understanding Darjeeling tea flavour on a molecular basis. Plant Mol Biol 2012, 78, 577-97. (3) Yang, Z. Y.; Baldermann, S.; Watanabe, N., Recent studies of the volatile compounds in tea. Food Res Int 2013, 53, 585-599. (4) Rawat, R.; Gulati, A.; Babu, G. D. K.; Acharya, R.; Kaul, V. K.; Singh, B., Characterization of volatile


Page 24 of 38

Page 25 of 38

533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576


components of Kangra orthodox black tea by gas chromatography-mass spectrometry. Food Chem 2007, 105, 229-235. (5) Xu, W.; Song, Q.; Li, D.; Wan, X., Discrimination of the production season of Chinese green tea by chemical analysis in combination with supervised pattern recognition. J Agric Food Chem 2012, 60, 7064-70. (6) Xu, Q.; He, Y.; Yan, X.; Zhao, S.; Zhu, J.; Wei, C., Unraveling a crosstalk regulatory network of temporal aroma accumulation in tea plant (Camellia sinensis) leaves by integration of metabolomics and transcriptomics. Environ Exp Bot 2018, 149, 81-94. (7) Schuh, C.; Schieberle, P., Characterization of the key aroma compounds in the beverage prepared from Darjeeling black tea: quantitative differences between tea leaves and infusion. J Agric Food Chem 2006, 54, 916-24. (8) Dubey, V. S.; Bhalla, R.; Luthra, R., An overview of the non-mevalonate pathway for terpenoid biosynthesis in plants. J Biosci 2003, 28, 637-46. (9) Zeng, L.; Watanabe, N.; Yang, Z., Understanding the biosyntheses and stress response mechanisms of aroma compounds in tea (Camellia sinensis) to safely and effectively improve tea aroma. Crit Rev Food Sci Nutr 2018, 1-14. (10) Wei, X.; Yan, X.; Fu-Min, W.; Li-Ping, G.; Ming-Jun, G.; Zheng-Zhu, Z.; Xiao-Chun, W.; Shu, W., In Silico Analysis and Feeding Assays of Some Genes in the Early Steps of Terpenoid Biosynthetic Pathway in Camellia Sinensis. Journal of Tea 2013, 39, 191-198. (11) Fu, J. Y., Molecular cloning and expression analysis of a putative sesquiterpene synthase gene from tea plant (Camellia sinensis). Acta Physiologiae Plantarum 2013, 35, 289-293. (12) Fu, X. M.; Chen, Y. Y.; Mei, X.; Katsuno, T.; Kobayashi, E.; Dong, F.; Watanabe, N.; Yang, Z. Y., Regulation of formation of volatile compounds of tea (Camellia sinensis) leaves by single light wavelength. Sci Rep 2015, 5. (13) Mei, X.; Liu, X. Y.; Zhou, Y.; Wang, X. Q.; Zeng, L. T.; Fu, X. M.; Li, J. L.; Tang, J. C.; Dong, F.; Yang, Z. Y., Formation and emission of linalool in tea (Camellia sinensis) leaves infested by tea green leafhopper (Empoasca (Matsumurasca) onukii Matsuda). Food Chem 2017, 237, 356-363. (14) Liu, G. F.; Liu, J. J.; He, Z. R.; Wang, F. M.; Yang, H.; Yan, Y. F.; Gao, M. J.; Gruber, M. Y.; Wan, X. C.; Wei, S., Implementation of CsLIS/NES in linalool biosynthesis involves transcript splicing regulation in Camellia sinensis. Plant Cell and Environment 2018, 41, 176-186. (15) Zhou, Y.; Zeng, L. T.; Liu, X. Y.; Gui, J. D.; Mei, X.; Fu, X. M.; Dong, F.; Tang, J. C.; Zhang, L. Y.; Yang, Z. Y., Formation of (E)-nerolidol in tea (Camellia sinensis) leaves exposed to multiple stresses during tea manufacturing. Food Chem 2017, 231, 78-86. (16) Bedon, F.; Bomal, C.; Caron, S.; Levasseur, C.; Boyle, B.; Mansfield, S. D.; Schmidt, A.; Gershenzon, J.; Grima-Pettenati, J.; Seguin, A.; MacKay, J., Subgroup 4 R2R3-MYBs in conifer trees: gene family expansion and contribution to the isoprenoid- and flavonoid-oriented responses. J Exp Bot 2010, 61, 3847-64. (17) Reddy, V. A.; Wang, Q.; Dhar, N.; Kumar, N.; Venkatesh, P. N.; Rajan, C.; Panicker, D.; Sridhar, V.; Mao, H. Z.; Sarojam, R., Spearmint R2R3-MYB transcription factor MsMYB negatively regulates monoterpene production and suppresses the expression of geranyl diphosphate synthase large subunit (MsGPPS.LSU). Plant Biotechnol J 2017, 15, 1105-1119. (18) Zvi, M. M.; Shklarman, E.; Masci, T.; Kalev, H.; Debener, T.; Shafir, S.; Ovadis, M.; Vainstein, A., PAP1 transcription factor enhances production of phenylpropanoid and terpenoid scent compounds in rose flowers. New Phytol 2012, 195, 335-45. (19) Mandaokar, A.; Thines, B.; Shin, B.; Lange, B. M.; Choi, G.; Koo, Y. J.; Yoo, Y. J.; Choi, Y. D.; Choi, G.; Browse, J., Transcriptional regulators of stamen development in Arabidopsis identified by transcriptional profiling. Plant J 2006, 46, 984-1008.



577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620

(20) Reeves, P. H.; Ellis, C. M.; Ploense, S. E.; Wu, M. F.; Yadav, V.; Tholl, D.; Chetelat, A.; Haupt, I.; Kennerley, B. J.; Hodgens, C.; Farmer, E. E.; Nagpal, P.; Reed, J. W., A regulatory network for coordinated flower maturation. PLoS Genet 2012, 8, e1002506. (21) Hong, G. J.; Xue, X. Y.; Mao, Y. B.; Wang, L. J.; Chen, X. Y., Arabidopsis MYC2 interacts with DELLA proteins in regulating sesquiterpene synthase gene expression. Plant Cell 2012, 24, 2635-48. (22) van der Fits, L.; Memelink, J., ORCA3, a jasmonate-responsive transcriptional regulator of plant primary and secondary metabolism. Science 2000, 289, 295-7. (23) Zhang, H.; Hedhili, S.; Montiel, G.; Zhang, Y.; Chatel, G.; Pre, M.; Gantet, P.; Memelink, J., The basic helix-loop-helix transcription factor CrMYC2 controls the jasmonate-responsive expression of the ORCA genes that regulate alkaloid biosynthesis in Catharanthus roseus. Plant J 2011, 67, 61-71. (24) Yu, Z. X.; Li, J. X.; Yang, C. Q.; Hu, W. L.; Wang, L. J.; Chen, X. Y., The jasmonate-responsive AP2/ERF transcription factors AaERF1 and AaERF2 positively regulate artemisinin biosynthesis in Artemisia annua L. Mol Plant 2012, 5, 353-65. (25) Xu, Y. H.; Wang, J. W.; Wang, S.; Wang, J. Y.; Chen, X. Y., Characterization of GaWRKY1, a cotton transcription factor that regulates the sesquiterpene synthase gene (+)-delta-cadinene synthase-A. Plant Physiol 2004, 135, 507-15. (26) Wang, Q.; Reddy, V. A.; Panicker, D.; Mao, H. Z.; Kumar, N.; Rajan, C.; Venkatesh, P. N.; Chua, N. H.; Sarojam, R., Metabolic engineering of terpene biosynthesis in plants using a trichome-specific transcription factor MsYABBY5 from spearmint (Mentha spicata). Plant Biotechnol J 2016, 14, 1619-32. (27) Yu, Z. X.; Wang, L. J.; Zhao, B.; Shan, C. M.; Zhang, Y. H.; Chen, D. F.; Chen, X. Y., Progressive regulation of sesquiterpene biosynthesis in Arabidopsis and Patchouli (Pogostemon cablin) by the miR156-targeted SPL transcription factors. Mol Plant 2015, 8, 98-110. (28) Liu, L.; Li, Y.; She, G.; Zhang, X.; Jordan, B.; Chen, Q.; Zhao, J.; Wan, X., Metabolite profiling and transcriptomic analyses reveal an essential role of UVR8-mediated signal transduction pathway in regulating flavonoid biosynthesis in tea plants (Camellia sinensis) in response to shading. BMC Plant Biol 2018, 18, 233. (29) Ambady, S.; Wu, Z.; Dominko, T., Identification of novel microRNAs in Xenopus laevis metaphase II arrested eggs. Genesis 2012, 50, 286-99. (30) Cer, R. Z.; Herrera-Galeano, J. E.; Anderson, J. J.; Bishop-Lilly, K. A.; Mokashi, V. P., miRNA Temporal Analyzer (mirnaTA): a bioinformatics tool for identifying differentially expressed microRNAs in temporal studies using normal quantile transformation. Gigascience 2014, 3, 20. (31) Li, X.; Shahid, M. Q.; Wu, J.; Wang, L.; Liu, X.; Lu, Y., Comparative Small RNA Analysis of Pollen Development in Autotetraploid and Diploid Rice. Int J Mol Sci 2016, 17, 499. (32) Addo-Quaye, C.; Eshoo, T. W.; Bartel, D. P.; Axtell, M. J., Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Curr Biol 2008, 18, 758-762. (33) Addo-Quaye, C.; Miller, W.; Axtell, M. J., CleaveLand: a pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics 2009, 25, 130-1. (34) Ma, Z.; Coruh, C.; Axtell, M. J., Arabidopsis lyrata small RNAs: transient MIRNA and small interfering RNA loci within the Arabidopsis genus. Plant Cell 2010, 22, 1090-103. (35) Addo-Quaye, C.; Snyder, J. A.; Park, Y. B.; Li, Y. F.; Sunkar, R.; Axtell, M. J., Sliced microRNA targets and precise loop-first processing of MIR319 hairpins revealed by analysis of the Physcomitrella patens degradome. RNA 2009, 15, 2112-21. (36) Yang, J. H.; Liu, X. Y.; Xu, B. C.; Zhao, N.; Yang, X. D.; Zhang, M. F., Identification of miRNAs and their targets using high-throughput sequencing and degradome analysis in cytoplasmic male-sterile and its maintainer fertile lines of brassica juncea. BMC Genomics 2013, 14.


Page 26 of 38

Page 27 of 38

621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664


(37) Langfelder, P.; Horvath, S., WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008, 9, 559. (38) Smoot, M. E.; Ono, K.; Ruscheinski, J.; Wang, P. L.; Ideker, T., Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics 2011, 27, 431-2. (39) Livak, K. J.; Schmittgen, T. D., Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 2001, 25, 402-8. (40) Jeyaraj, A.; Zhang, X.; Hou, Y.; Shangguan, M.; Gajjeraman, P.; Li, Y.; Wei, C., Genome-wide identification of conserved and novel microRNAs in one bud and two tender leaves of tea plant (Camellia sinensis) by small RNA sequencing, microarray-based hybridization and genome survey scaffold sequences. BMC Plant Biol 2017, 17, 212. (41) Rio, D. C., Northern blots for small RNAs and microRNAs. Cold Spring Harb Protoc 2014, 2014, 793-7. (42) Samad, A. F. A.; Sajad, M.; Nazaruddin, N.; Fauzi, I. A.; Murad, A. M. A.; Zainal, Z.; Ismail, I., MicroRNA and Transcription Factor: Key Players in Plant Regulatory Network. Frontiers in Plant Science 2017, 8. (43) Gupta, O. P.; Karkute, S. G.; Banerjee, S.; Meena, N. L.; Dahuja, A., Contemporary Understanding of miRNA-Based Regulation of Secondary Metabolites Biosynthesis in Plants. Frontiers in Plant Science 2017, 8. (44) CAO, F.-r.; LIU, K.-b.; LIU, C.-y.; WANG, D.-l., Studies on the Induction of Aromatic Constituents in Fresh Leaves of Lingtou Dancong Tea by Low Temperature Stress [J]. Journal of Tea Science 2006, 2, 011. (45) Gouinguene, S. P.; Turlings, T. C. J., The effects of abiotic factors on induced volatile emissions in corn plants. Plant Physiol 2002, 129, 1296-1307. (46) Staudt, M.; Bertin, N., Light and temperature dependence of the emission of cyclic and acyclic monoterpenes from holm oak (Quercus ilex L.) leaves. Plant, Cell Environ 1998, 21, 385-395. (47) GLEIZES, M.; Pauly, G.; Bernard‐Dagan, C.; Jacques, R., Effects of light on terpene hydrocarbon synthesis in Pinus pinaster. Physiol Plant 1980, 50, 16-20. (48) Sasaki, K.; Ohara, K.; Yazaki, K., Gene expression and characterization of isoprene synthase from Populus alba. FEBS Lett 2005, 579, 2514-2518. (49) Stief, A.; Altmann, S.; Hoffmann, K.; Pant, B. D.; Scheible, W. R.; Baurle, I., Arabidopsis miR156 Regulates Tolerance to Recurring Environmental Stress through SPL Transcription Factors. Plant Cell 2014, 26, 1792-1807. (50) Kawoosa, T.; Singh, H.; Kumar, A.; Sharma, S. K.; Devi, K.; Dutt, S.; Vats, S. K.; Sharma, M.; Ahuja, P. S.; Kumar, S., Light and temperature regulated terpene biosynthesis: hepatoprotective monoterpene picroside accumulation in Picrorhiza kurrooa. Funct Integr Genomics 2010, 10, 393-404. (51) Sun, Z.; Li, M.; Zhou, Y.; Guo, T.; Liu, Y.; Zhang, H.; Fang, Y., Coordinated regulation of Arabidopsis microRNA biogenesis and red light signaling through Dicer-like 1 and phytochrome-interacting factor 4. PLoS Genet 2018, 14, e1007247. (52) Sun, Z.; Guo, T.; Liu, Y.; Liu, Q.; Fang, Y., The Roles of Arabidopsis CDF2 in Transcriptional and Posttranscriptional Regulation of Primary MicroRNAs. PLoS Genet 2015, 11, e1005598. (53) Cho, S. K.; Ben Chaabane, S.; Shah, P.; Poulsen, C. P.; Yang, S. W., COP1 E3 ligase protects HYL1 to retain microRNA biogenesis. Nat Commun 2014, 5, 5867. (54) Osterlund, M. T.; Wei, N.; Deng, X. W., The roles of photoreceptor systems and the COP1-targeted destabilization of HY5 in light control of Arabidopsis seedling development. Plant Physiol 2000, 124, 1520-4. (55) Osterlund, M. T.; Hardtke, C. S.; Wei, N.; Deng, X. W., Targeted destabilization of HY5 during light-regulated development of Arabidopsis. Nature 2000, 405, 462-6. (56) Zhang, H. Y.; He, H.; Wang, X. C.; Wang, X. F.; Yang, X. Z.; Li, L.; Deng, X. W., Genome-wide mapping of the HY5-mediated genenetworks in Arabidopsis that involve both transcriptional and post-transcriptional regulation. Plant J 2011, 65, 346-358.



665 666 667 668 669 670

(57) Andersson, C. R.; Kay, S. A., COP1 and HY5 interact to mediate light‐induced gene expression. Bioessays

671

Table and Figure captions

672

Table 1. The differentially expressed targets of DEmiRNAs validated by degradome

673

sequencing which involved in secondary metabolism.

674

*The number in brackets represents the five classes (with 0 as the highest degradome peak) in

675

which the cleaved target transcripts were categorized based on their signature abundance at

676

each occupied transcript position.

677

Figure 1. Distribution of miRNAs in different growing time points.

678

Length distribution of total miRNAs (A) and Venn diagram of detected miRNAs (B) in

679

different growing months by high-throughput sequencing. The numbers outside the brackets

680

refer to sum of known and unknown miRNAs, and the numbers within the brackets refer to

681

the unknown miRNAs.

682

Figure 2. Summary of the distribution of five groups of miRNAs. CountsMIRb: mapped

683

miRNA numbers in miRBase. Expression level: low, 10 but less than average;

684

high, over average.

685

Figure 3. Heat map showing the expression pattern of differentially expressed miRNAs in

686

different growing months of tea tree. The expression values of miRNAs were normalized by

687

Z-score normalization. The high- and low-expression DEmiRNAs at each time were showed

688

in A and B respectively.

689

Figure 4. mRNA cleavage sites revealed by degradome sequencing and verified by 5’

1998, 20, 445-448. (58) Ang, L.-H.; Chattopadhyay, S.; Wei, N.; Oyama, T.; Okada, K.; Batschauer, A.; Deng, X.-W., Molecular interaction between COP1 and HY5 defines a regulatory switch for light control of Arabidopsis development. Mol Cell 1998, 1, 213-222.


Page 28 of 38

Page 29 of 38


690

RLM-RACE.

691

The mRNA target and its corresponding miRNA are shown in each alignment. Watson–Crick

692

base pairs and G–U base pairs are indicated by “:” and “.”, respectively. Red arrows indicate

693

cleavage sites, and the clone frequencies are shown. In the target plots (t-plots). The X-axis

694

indicates the mRNA position of miRNA targets from 5’ to 3’. The red bar shows the number

695

of reads in degradome of mRNA cleavage site.

696

Figure 5. Analysis of enriched GO (A) and KEGG (B) of target genes for the identified

697

miRNAs.

698

Figure 6. A general negative expression correlation between differentially expressed

699

miRNAs (A) and their corresponding differentially expressed target genes (B) (cor

Revealing of the microRNA involved regulatory gene networks on

Recommend Documents