- Email: [email protected]
Identification of drought resistant miRNA in Macleaya cordata by high-throughput sequencing
Journal Pre-proof Identification of drought resistant miRNA in Macleaya cordata by high-throughput sequencing Linlan Yu, Li Zhou, Wei Liu, Peng Huang, Ruolan Jiang, Zhaoshan Tang, Pi Cheng, Jianguo Zeng PII:
S0003-9861(19)30922-1
DOI:
https://doi.org/10.1016/j.abb.2020.108300
Reference:
YABBI 108300
To appear in:
Archives of Biochemistry and Biophysics
Received Date: 14 October 2019 Revised Date:
15 January 2020
Accepted Date: 8 February 2020
Please cite this article as: L. Yu, L. Zhou, W. Liu, P. Huang, R. Jiang, Z. Tang, P. Cheng, J. Zeng, Identification of drought resistant miRNA in Macleaya cordata by high-throughput sequencing, Archives of Biochemistry and Biophysics (2020), doi: https://doi.org/10.1016/j.abb.2020.108300. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Inc.
1
Identification of drought resistant miRNA in Macleaya
2
cordata by high-throughput sequencing
3 4 5 6 7 8 9 10 11 12 13 14 15 16
Linlan Yua 2 Li Zhoua 2 Wei Liua 1,3 Peng Huang1,2 Ruolan Jiang1,2 Zhaoshan Tang5 Pi Cheng*1,2
17 18
Wei Liu#: [email protected] Peng Huang: [email protected]
19 20 21 22 23 24 25
Ruolan Jiang: [email protected]
Jianguo Zeng*1,4 1
Hunan Key Laboratory of Traditional Chinese Veterinary Medicine, Hunan Agricultural University,
Changsha, 410128 Hunan China 2
College of Horticulture and Landscape, Hunan Agricultural University, Changsha, 410128 Hunan
China 3
Center of Analytic Service, Hunan Agriculture University, 410208 Changsha, China
4
National and Local Union Engineering Research Center of Veterinary Herbal Medicine Resource and
Initiative, Hunan Agricultural University, Changsha, 410128 Hunan China 5
Micolta Bioresource Inc., Changsha, 410016 China
Email Linlan Yu#: [email protected] Li Zhou#: [email protected]
Zhaoshan Tang: [email protected] Pi Cheng*: [email protected]
Jianguo Zeng*: [email protected] * Corresponding author a These authors contributed equally to this work.
26
Abstract:
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
Drought is one of the most serious factors affecting crop yields in the world. Macleaya cordata (Willd.) is a draught-tolerant medicinal plant that has been proposed as a pioneer crop to be cultivated in arid areas. However, the exact molecular mechanisms through which M. cordata responds to draught stress remain elusive. In recent years, microRNA (miRNAs) in plants have been associated with stress response. Based on these findings, the current study aimed to shed light on the potential regulatory roles of miRNAs in the draught tolerance of M. cordata by employing high-throughput RNA sequencing and degradation sequencing. Six M. cordata plants were randomly divided into two equal experiment groups, including one draught group and one control group. High-throughput sequencing of the M. cordata samples led to the identification of 895 miRNAs, of which 18 showed significantly different expression levels between the two groups. PsRobot analysis and degradation sequencing predicted the differential miRNAs to target 59 and 36 genes, respectively. Functional analysis showed that 38 of the predicted genes could be implicated in the modulation of stress response. Four miRNAs and eight target genes were selected for quantitative real-time polymerase chain reaction (qRT-PCR) validation. The expression trend of each miRNA analyzed by qRT-PCR was consistent with that determined by sequencing, and was negatively correlated with those of its target genes.The results of our current study supported the involvement of miRNAs in the draught tolerance of M. cordata and could pave the way for further investigation into the related regultory mechanisms.
48
Keywords: Macleaya cordata
49
sequencing
50
miRNAs
degradation sequencing
drought
stress
high-throughput
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
Background Drought is one of most serious type of abiotic stress for crops worldwide [1]. Studies have shown that draught conditions can cause the temperature of plant cells to increase and reactive oxygen species (ROS) to accumulate, leading to disrupted cell metabolism and even apoptosis [2]. Furthermore, extreme draught stress can also trigger the closure of stomata, which significantly hampers photosynthesis [3]. Not surprisingly, plants have evolved a wide range of molecular and physiological mechanisms to help them resist draught stress and mitigate the detrimental effects of water deficit [4]. For example, studies on Arabidopsis thaliana have revealed that a number of kinases, including SPK1 and MAPK, are activated by draught-induced up-regulation of abscisic acid (ABA) and can subsequently promote physiological changes designed to prevent excessive water loss [5, 6]. Additionally, there is also mounting evidence that shifts in the biosynthesis of certain cell wall components, such as cellulose and xyloglucan, constitute another strategy that plants employ in response to draught [7]. Better understanding of these response mechanisms would greatly facilitate the development of novel draught-tolerant crop cultivars with enhanced economic value. Macleaya cordata (Willd.) R. Br. is a herbaceous perennial mainly distributed across Southeast China and Southeast Asia [8, 9]. Extracts prepared from M. cordata are an important ingredient in traditional Chinese medicine with demonstrated antimicrobial, insecticidal, antitumoral and anti-inflammatory effects [10-12]. Based on chemical analysis, the main bioactive constituents of the plant consist of isoquinoline alkaloids such as sanguinarine, chelerythrine, protopine and allocryptopine. Interestingly, M. cordata has been shown to exhibit moderate draught tolerance and thus proposed as a pioneer plant species for soil conditioning in arid lands. We have recently discovered that M. cordata exposed to simulated draught stress with 30% PEG6000 showed no observable decrease in growth rate and produced 170% more sanguinarine over 48 h compared to the controls cultivated under regular conditions [11]. Despite these results, the draught response mechanisms in M. cordata remain elusive and require further studies to decipher. MicroRNAs (miRNA) are a diverse family of endogenous single-stranded RNA molecules consisting of 21 to 24 nucleotides [13]. Although miRNA sequences do not encode proteins themselves, they have been increasingly recognized as critical regulators of gene expression [14]. In general, miRNAs are capable of inhibiting the translation, or even directly inducing the degradation, of the target mRNAs by fully or partially complementing with their 3’ untranslated regions (UTR) [15]. Recently, miRNAs have also been shown to participate in plant response to various abiotic stress signals, such as draught, heat, oxidative stress and heavy-metal toxicity [16]. The rapid advances of high-throughput sequencing technologies and bioinformatic algorithms have provided a boon for miRNA research, enabling the systematic examination of miRNAs expression profiles in the same host under different environmental or physiological conditions. Combined, these developments have served as the foundation for our current study, in which we sought to identify miRNAs associated with the draught tolerant traits of M. cordata. Comparison of the miRNAs expression profiles between the drought-stressed treatment group and the control group allowed us to identify 18 differentially expressed candidates. Further functional analysis suggested that the differentially expressed miRNAs were involved in biological functions such as carboxylic ester hydrolase activity, phosphorus metabolic process, membrane. Together,
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138
these findings could increase our understanding of the complex regulatory network that underlies the draught tolerance of M. cordata, and provide a theory basis for encouraging its cultivation in arid regions for vegetation restoration.
Methods Draught induction M. cordata were collected from Huaihua, China with no permission required, and cultivated in Hunan Agricultural University. Six M. cordata plantlets derived from the same explant were grown at 25 °C in Murashige-Skoog (MS) medium (4.42 g/L) supplemented with 30 g/L sucrose and 3 g/L Phytagel at a pH of 5.8 under a 16/8-h light/dark photoperiod for 30 days. The plantlets were randomly divided into two equal groups, including a drought-stressed treatment group and a control group. The three plantlets in the draught group were subjected to draught stress by transfer to a new MS medium containing 30% (w/v) PEG6000, whereas the control group were cultivated in fresh MS medium. After 3 days, the plantlets were harvested, ground to a fine powder in liquid nitrogen, and then stored at -80 °C.
Library construction Total RNA was extracted from 100 mg of ground plant materials using Plant RNA Purification Kit (Thermo Fisher, USA) according the manufacturer’s instructions, followed by the removal of genomic DNA with RNase-free recombinant DNase I (TaKaRa, Japan). The quality of the extracted RNA was verified on a 2100 Bioanalyzer (Agilent Technologies, USA) and a NanoDrop 2000 spectrometer (Thermo Fisher Scientific, USA). Only high-quality RNA, defined as meeting the criteria of OD260/280 = 1.8 ~ 2.2, OD 260/230 ≥ 2.0, RIN ≥ 6.5 and 28S:18S ≥ 1.0, was used to construct a sequencing library.
Library preparation, and Illumina Hiseq xten Sequencing Stranded RNA-seq library for transcriptome analysis was prepared from 5 μg of high-quality total RNA using a Sample Prep Kit (Illumina, CA). First-stranded cDNA was synthesized with random hexamer primers, followed by template removal and replacement strand synthesis with the incorporation of dUTP in place of dTTP. As the polymerase could not utilize dUTP, the elongation of the complementary strand was aborted prematurely, leading to incomplete blunt-ended double-stranded cDNA fragments of varying lengths. These fragments were subsequently purified by AMPure XP beads. A single “A” nucleotide was added to each 3’ end of these fragments to prevent blunt end ligation. The modified cDNA fragments were then ligated with multiple indexing adapters and ollowed by PCR amplification with Phusion DNA polymerase (New England Biolabs, USA) for 15 cycles. The resultant cDNA library was quantified on a TBS-380 Mini-fluorometer (Turner Biosystems, USA) and sequenced on a HiSeq X Ten system (2 × 150bp reads, Illumina, USA).
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
Identification and expression analysis of identified miRNAs The 3’ ends of all low-quality bases with a quality score below 20 were trimmed using in-house PERL scripts, and the adapters were then removed with FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) [17]. Among the remaining sequences, those with 18 to 32 nucleotides were retained. The assembled unique sequences were used as queries for a BLAST search against the Rfam database (version 12.1) [18] to remove non-miRNA members such as rRNA, tRNA and snoRNA. Bowtie [19] (version 1.2.2) was used to map the predicted miRNAs)against the reference genome and predict their secondary structures. On the one hand, known miRNAs were identified by searching miRbase [20] (version 21.0, http://www.mirbase.org/) and retaining the perfectly matched sequences. On the other hand, novel candidates were predicted by MIREAP [21] based on the characteristic hairpin structure of miRNA precursors. In addition, in-house scripts were used to determine the nucleotide bias at each position of the identified miRNAs. The expression level of each miRNA was calculated by miRDeep2 [22] (version 2.0.0.8) in the form of transcript per million reads (TPM). Differential expression analysis was conducted by using the DESeq2 package [23] (version 1.6.3, http://www.bioconductor.org/packages/release/bioc/html/DESeq2.htmls) based on the criteria of |log2(fold change)| > 1 and false discovery rate-adjusted p value < 0.05.
Annotation and functional analysis The target genes of all identified mRNAs were predicted and annotated with psRobot [24]. Gene ontology (GO) analysis was performed by querying the genes against the GO database (http://www.geneontology.org/) [25] using GOATOOLS (version 0.5.9, https://github.com/tanghaibao/GOatools) [26]. GOATOOLS employed Fisher’s exact test to compute the p value that indicates whether the frequency of the target genes assigned to a specific GO term is significantly different from that of randomly selected background genes. The p value was subsequently corrected using four different methods, including Bonferroni, Holm, Sidak and False Discovery Rate. Meanwhile, Kyoto Encyclopedia of Genes and Genomes (KEGG) [27] pathway analysis was performed by querying the genes against the KEGG database (http://www.genome.jp/kegg/) with KOBAS [28]. P value was calculated by applying Fisher’s exact test and then corrected by the Benjamini-Hochberg method.
Degradation sequencing The mRNA targets of the miRNAs were verified by degradation sequencing. After binding with miRNA, the mRNA is cleaved and the resultant 3'-fragment, which contains a free 5'-monophosphate and a 3' polyA tail, can be extended by RNA ligase and subsequently sequenced. In-depth alignment analysis of the sequenced data can then reveal the candidate miRNA splicing sites in the mRNA sequence. Briefly, mRNA molecules were enriched from total
183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
RNA by Total RNA Purification Kit (LC Science, Houston, USA) and then subjected to 5' adaptor ligation, followed by reverse transcription with mixed (PrimeScript®1st Strand cDNA Synthesis Kit , TaKaRa, Japan) and (PCR Ampification Kit Manual, Takara, Japan) according to the manufacturer’s instructions. The generated library was sequenced on a Hiseq 2500 system (Illumina, USA).
Quantitative real-time PCR (qRT-PCR) validation Selected miRNA candidates were reverse-transcribed with RevertAid Premium Reverse Transcriptase (Thermo Fisher Scientific, USA) following the manufacturer’s instructions. The synthesized cDNAs were diluted ten-fold and used directly as template. Each reaction consisted 10 µM each primer and 1 μg of template in 1 × Fast SG qPCR Master Mix (Sangon Biotech, China). The qRT-PCR amplification was performed on a StepOnePlus Real Time PCR System (Applied Biosystems, USA) programmed as follows: 95 °C for 3 min, followed by 45 cycles of 95 °C for 3 s and 60 °C for 30 s. The target mRNAs were reverse-transcribed into cDNAs with Prime Script RT Reagent Kit with gDNA Eraser (Perfect Real Time, TaKaRa, Japan) following the manufacturer’s instructions. Each reaction consisted of 3 µM of each primer and 1 μg of template in 1× FastStart Universal SYBR Green Master (Roche, USA). The qRT-PCR amplification was performed on a 7300 Fluorescence Quantitative PCR Analyzer (Applied Biosystems, USA) programmed as follows: 95°C for 5 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 60 s. All PCR primers used in this study were summarized in Table 2.
Results
Overview of the high-throughput sequencing results Our high-throughput sequencing generated a total of 101375907 raw reads from the six libraries. All low-quality reads were then removed, and only reads consisting of 18 – 32 nucleotides were retained. Eventually, we obtained 80957233 clean reads. Approximately 44.33% of these reads were successfully matched to Rfam. The matching suggested that roughly 23.39% of the clean reads were likely derived from miRNA sequences, with the rest comprising rRNAs, tRNAs, snRNAs and other unknown RNA molecules. The percentage of miRNAs in the clean reads derived from each sequencing library ranged from 17.3% to 29.3%. Detailed descriptions of the various types of RNAs that we identified from the M. cordata samples were provided in Table 1 and additional file 1.
Identification of miRNAs
227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270
Comparison with miRBase-deposited sequences allowed us to identify 355 miRNA candidates belonging to 68 miRNA families. The expression level of each miRNA was then estimated by miRDeep2. Notably, the 13 most highly expressed miRNAs from each sample were all from the miR166 family. Meanwhile, we also predicted an additional 540 novel miRNAs belonging to 334 families. These novel candidates generally showed much lower expression levels than the known miRNAs described above. It is worth emphasizing that the miRNA expression profiles varied significantly among different samples. A comparative analysis indicated that 337 known and 294 predicted miRNA candidates were detected in all six M. cordata samples. In contrast, 36 known and 112 novel miRNAs were found only in the draught group, whereas 42 known and 136 novel candidates were exclusively detected in the control samples.
Differential expression analysis We then compared the miRNA expression profiles between the draught group and the control group. A miRNA is considered differentially expressed if there was a two-fold difference in its expression between the two experiment groups and if the corresponding FDR-adjusted p value was below 0.05. Based on these criteria, nine known and nine novel miRNAs exhibited differential expression between the draught group and the control group. Among them, six known (ath-miR8175, ath-miR827, cca-miR395b, gma-miR395d, gma-miR395i, mtr-miR395h) and five novel (Nov-m0067-5p, Nov-m0176-5p, Nov-m0231-5p, Nov-m0383-5p, Nov-m0396-3p) miRNAs were up-regulated in the draught group, whereas the rest were down-regulated (known miRNAs: cca-miR396a-3p, cca-miR396c and stu-miR156d-3p; novel miRNAs: Nov-m0038-3p, Nov-m0112-5p, Nov-m0219-5p and Nov-m0411-3p). The expression of all detected miRNAs was summarized in figure 1 and 2. Cluster analysis also provided strong evidence that there were clear distinctions between the miRNA expression patterns of the two experiment groups (Fig. 3 and 4).
Target prediction and functional analysis of the drought-responsive miRNAs
The target genes of the identified differentially expressed miRNAs were predicted by psRobot [24]. The nine known miRNAs were predicted to interact with a total of 11 target genes. Target multiplicity, a hallmark of miRNA-mRNA interaction, was unambiguously demonstrated by the observation that most of the miRNAs were predicted to complement the 3’ UTR of more than one mRNA and nearly half of the identified genes were recognized by two or more miRNAs. Similarly, the nine novel miRNA candidates were shown to target 48 genes. In particular, Nov-m0067-5p was predicted to interact with 40 different mRNAs. Nov-m0411-3p could potentially complement the 3’-UTRs of 6 mRNAs, whereas Nov-m0219-5p and Nov-m0396-3p were found to each target one gene. In sharp contrast to what was observed of the known miRNAs, none of the identified genes was predicted to be the target of more than one of the differentially expressed novel miRNA candidates. Known miRNAs such as miR396, miR395,
271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314
miR156-3 and miR827 were predicted to respond to drought by regulating the target genes UBC, ADH5, RPL8 and ALDH. Meanwhile, Nov-miR-0067, Nov-miR-411, Nov-m0383-5p and Nov-m0396-3p were predicted to respond to drought stress by regulating the target genes PLD, GST, ABCB1, SOD2, and APRT.
Annotation and functional analysis of differentially expressed genes Annotation of the combined 59 target genes revealed that they encode a wide assortment of growth- and/or stress-related proteins, such as phospholipase D, ABC transporter, GST, phosphotransferase, kinase, transcription factor, GDSL esterase, zinc finger protein, and pherophorin. The annotated genes were subsequently subjected to GO enrichment and KEGG pathway analyses. The GO annotations of the genes were predicted by GOATOOLS (version: 0.5.9, https://github.com/tanghaibao/GoOatools) and classified into three sub-categories, including molecular function, biological process or cellular component. Based on these results, the GO terms that could be meaningfully interpreted and associated with stress response and/or draught tolerance included hydrolase activity, acting on ester bonds (GO: 0016788), hydrolase activity (GO:0016787) , primary metabolic process (GO:0044238), organic substance metabolic process (GO:0071704), carboxylic ester hydrolase activity (GO:0052689), small molecule binding (GO:0036094), , phosphate-containing compound metabolic process (GO:0006796) and phosphorus metabolic process (GO:0006793). On the other hand, KEGG analysis indicated that the putative target genes of the differentially expressed miRNAs were significantly enriched in a number of pathways that could contribute to plant stress response and/or draught tolerance, including GnRH signaling pathway (ko04912), ABC transporters (ko02010), Ras signaling pathway (ko04014), Glutathione metabolism (ko00480), Plant hormone signal transduction (ko04075) and Phospholipase D signaling pathway (ko04072).
Target validation by degradation sequencing In total, we predicted that the 895 miRNAs that we detected in the M. cordata samples could target a total of 22691 genes. We then sought to systematically verify these miRNA-mRNA interactions through degradation sequencing, which can selectively amplify the 3’ cleavage fragments of miRNAs. The nucleotide sequences of these fragments were then compared against the corresponding miRNA sequences. In case of perfect or near-perfect complementation, the cleaved mRNA would be considered as a target of the corresponding miRNA. Based on these principles, we identified 8022 mRNAs that produced 3’ fragments with the correct cleavage pattern indicative of miRNA binding. These mRNAs were collectively targeted by a combined 550 miRNAs. Importantly, 6310 genes were predicted only by psRobot, whereas 8022 were captured only by degradation sequencing. In the degradation group sequencing, multiple fragments of SPL (SQUAMOSA promoter binding like), GRF (growth regulating factor), and SCL (Scarecrow-like) that were cut by miR156 / 157, miR396 and miR170 / 171 were detected. SPL, GRF and SCL are the star molecules for drought resistance. In addition, KEGG function enrichment was performed on the detected target genes, among which the ratio related to adaptation environment was the
315 316 317 318 319 320 321 322
highest, KEGG function enrichment was performed on the detected target genes, among which the ratio related to adaptation to the environment was the highest, and the signal transduction was also a high proportion. These pathways are closely related to resistance to stress.
323 324 325
miR156d, aqc-miR166a and aqc-miR171a) and sixteen of their target genes (GST, ABCB1, PHOT2, MPC1, At4g37250, ADH5, APS1, PAP1, GRF3, GRF6, SPL6, SPL7, ATHB-15, REV, SCL6 and SCL15) that have previously been associated with draught tolerance for qRT-PCR validation.
326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358
Among them, four miRNAs (cca-miR395b, aqc-miR166a,aqc-miR171a and Nov-m0067-5p) were up-regulated in drought stress samples, and four (cca-miR396a-3p, cca-miR396a-5p, miR156d and Nov-m0411-3p) showed opposite trends. It is worth mentioning that seven of the chosen target genes, including GRF3, GRF6, SPL6, SPL7, At4g37250, PHOT2 and MPC1, were detected by degradation sequencing but not predicted by bioinformatic analysis. The results indicated that the expression patterns of all six miRNAs determined by qRT-PCR were consistent with those shown by RNA sequencing (Fig.5.). In addition, the expression trends of all target genes were found to be inversely correlated with those of their putative miRNA regulators. Among the target genes that we verified, PHOT2 showed extremely significantly augmented expression levels, GST, ABCB1, GRF3, GRF6, SPL6, SPL7, and MPC1 (shaggy-related protein kinase eta) showed significant increased, whereas At4g37250, ADH5, APS1 and PAP1 were extremely significantly down-regulated in the draught group, the rest down-regulated obviously. These results demonstrated that the RNA sequencing data that we obtained were sufficiently reliable and accurate for functional analysis and mechanistic interpretation.
qRT-PCR validation We selected four differentially expressed miRNAs (Nov-m0411-3p, cca-miR396a-3p, Nov-m0067-5p and cca-miR395b), four miRNAs known to regulate drought (cca-miR396a-5p,
Discussion MiRNAs have emerged as key regulators of abiotic stress in plants. As an example, knockdown of miRNA-166 or stimulation of its target gene homeodomain containing protein 4 (OsHB4) was found to improve draught resistance in rice by inducing the development of smaller bulliform cells and abnormal sclerenchymatous cells [29]. In another study, up-regulation of miR-169a or suppression of its target gene nuclear factor -YA (NFYA5) in Arabidopsis thaliana aggravated water loss through the stroma [30]. Overexpression of miR-319 and Osa-miR-319a in Agrostis stolonifera has been shown to enhance its tolerance for drought and salt stress [31]. These findings prompted researchers to speculate the possibility of manipulating miRNA expression to engineer plant cultivars with enhanced stress resistance. A number of differentially expressed miRNAs that we identified have also been shown in previous studies to be potentially implicated in draught resistance. However, there is clear evidence that the expression trends and the underlying regulatory mechanisms of these miRNAs could vary significantly among different plant species. For example, miR-396 was found to be down-regulated in both M. cordata, as evidenced in our current study, and Oryza sativa [32], in
359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402
response to water deficit, but exhibited enhanced expression in draught-stressed Arabidopsis thaliana [33]. On the other hand, miR-396 has been reported to regulate stromatal density and plant transpiration in Arabidopsis thaliana and tobacco by targeting several growth regulating factor genes [34, 35]. Arabidopsis thaliana overexpressing miR-396 was characterized by linear narrow leaves with reduced stomatal density that showed greater resistance to water deficit [36]. However, the opposite effect was evident in Arabidopsis thaliana root [37]. Once again, our degradation sequencing suggested, GRFS was very likely due to the modulation of miR-396. Similarly, miR-156, another down-regulated miRNA candidate in our study, showed decreased expression in response to drought stress in Oryza sativa] [32] and Populus tomentosa [38], but the opposite pattern was observed in peach [39] and Arabidopsis thaliana [33]. Members from the SQUAMOSA promoter binding protein (SPL) family have been shown to be regulated by miR-156 [40]. SPLs are known to help plants combat a variety of abiotic stresses by modulating the anthocyanin synthesis pathway. Therefore, inhibition of miR-156 would de-repress SPLs and concomitantly enhance the draught tolerance of the host. This was consistent with our degradation sequencing data that SPL2, SPL3, SPL6, SPL7 and SPL12 were targeted by miR-156 in draught-stressed M. cordata samples. Based on the above, we speculate that miR396 and miR156 participate in drought stress response by cutting GRFs and SPLs. We also carefully investigated the predicted target genes of the novel miRNAs that exhibited differential expression between the draught group and the control group. For example, GST, a putative target of NOV-m0411-3p, was significantly up-regulated in draught-exposed M. cordata and is a well-established modulator of oxidative stress. GST was reported to contribute to the development of draught resistance in barley [41]. Consistently, overexpression of OsGSTU4 [42] and introduction of the tomato GST gene[43]] could both significantly enhance the draught tolerance of Arabidopsis thaliana. In the latter case, notably, higher GST level resulted in increased production of proline, malondialdehyde and antioxidative enzymes in the engineered plants [44]. The function of ATP-binding cassette transporters (ABC transporters) is transmembrane transport of ions, macromolecules and other substances through energy-driven [45]. ABCB1 was another target gene that we validated and functions as an auxin transporter. It can be activated with other proteins to participate in ion homeostasis in the cytoplasm in response to stress [46]. Studies have shown that the AtABCC5 mutant regulates ABA transport, regulates stomata, reduces transpiration rate, thereby reducing water consumption and enhancing drought tolerance [47]. AtABCB14 has similar functions [48]. As a prominent osmotic regulator, auxin plays a central role in plant response to abiotic stress signals such as water shortage and high salinity [49]. Therefore, it is unsurprising that auxin transporters have been associated with drought tolerance in both maize and Zea mays [50]. Notably, one of the plausible mechanisms through which auxin transporters might modulate draught resistance pathways was postulated to involve their stimulation of ABA production [51]. Taken together, our results suggested that, Nov-m0411-3p, which targets both GST and ABCB1, could be a key miRNA contributor to draught tolerance in M. cordata.
Conclusions
403 404 405 406 407 408 409 410 411 412 413 414 415 416
417 418 419 420 421 422 423 424 425 426 427
In summary, we provided the first systematic study that compared the miRNA expression profiles between draught-stressed M. cordata plants and controls cultivated under normal growth conditions. With a combination of high-throughput RNA sequencing and bioinformatic analysis, we identified a total of 895 miRNAs that showed differential expression as a result of induced water deficit. PsRobot-based bioinformatic analysis and degradation sequencing suggested that these miRNAs could target 6310 and 8022 genes, respectively. Based on functional analysis, these gene targets could be involved in a variety of pathways associated with plant response to draught, such as glutathione metabolism, plant hormone signal transduction and sulfur metabolism. Combined, our findings could facilitate the elucidation of molecular mechanisms that underlie the draught tolerance traits of M. cordata.
Abbreviations Full name
Abbreviations
microRNA 3’ untranslated regions Sphingosine kinase 1 Mitogen-activated protein kinase Abscisic acid ATP-binding cassette subfamily B member 1 Glutathione S-transferase Phototropin 2 Mitochondrial pyruvate carrier 1 Alcohol dehydrogenase 5 purple acid phosphatase 1 SQUAMOSA promoter binding like Growth regulating factor Arabidopsis thaliana homeobox genes REVOLUTA Scarecrow-like gene
miRNA UTR SPK1 MAPK ABA ABCB1 GST PHOT2 MPC1 ADH5 PAP 1 SPL GRF ATHB REVT SCL
Declarations Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Availability of data and materials
428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471
The data sets supporting the results of this article are included within the article and its additional files. Gut metagenome sequences were deposited in the National Center for Biotechnology Information (NCBI) SRP140359. Funding This work was supported by the innovation guidance program project of Hunan provincial science and Technology Department (2016SK3002), the National key laboratory cultivation base construction project (16KFXM09), and the Natural Science Foundation of Hunan Province (No. 2017JJ2115). Authors’ contributions L.Y., P. C. and J.Z. conceived and designed research. L.Y., L.Z., W.L. and P.H. onducted experiments. R.J. and Z.T. analyzed data. L.Y. wrote the manuscript. All authors read and approved the manuscript.
Competing Interests
The authors declare no competing interests.
References
[1] L. Cattivelli, F. Rizza, F.-W. Badeck, E. Mazzucotelli, A.M. Mastrangelo, E. Francia, C. Marè, A. Tondelli, A.M. Stanca, Drought tolerance improvement in crop plants: An integrated view from breeding to genomics, Field Crops Research 105(1-2)
0-14.
[2] M. Farooq, S. Basra, A. Wahid, Z. Cheema, M. Cheema, A. Khaliq, Physiological role of exogenously applied glycinebetaine to improve drought tolerance in fine grain aromatic rice (Oryza sativa L.), Journal of Agronomy and Crop Science 194(5) (2008) 325-333. [3] G. Cornic, Drought stress inhibits photosynthesis by decreasing stomatal aperture–not by affecting ATP synthesis, Trends in plant science 5(5) (2000) 187-188. [4] B. Valliyodan, H.T. Nguyen, Understanding regulatory networks and engineering for enhanced drought tolerance in plants, Curr Opin Plant Biol 9(2) (2006) 189-195. [5] Y. Osakabe, K. Maruyama, M. Seki, M. Satou, K. Shinozaki, K. Yamaguchi-Shinozaki, Leucine-rich repeat receptor-like kinase1 is a key membrane-bound regulator of abscisic acid early signaling in Arabidopsis, The Plant Cell 17(4) (2005) 1105-1119. [6] C. Jonak, S. Kiegerl, W. Ligterink, P.J. Barker, N.S. Huskisson, H. Hirt, Stress signaling in plants: a mitogen-activated protein kinase pathway is activated by cold and drought, Proceedings of the National Academy of Sciences 93(20) (1996) 11274-11279. [7] H. Le Gall, F. Philippe, J.-M. Domon, F. Gillet, J. Pelloux, C. Rayon, Cell wall metabolism in response to abiotic stress, Plants 4(1) (2015) 112-166. [8] P. Kosina, J. Gregorova, J. Gruz, J. Vacek, M. Kolar, M. Vogel, W. Roos, K. Naumann, V. Simanek, J. Ulrichova, Phytochemical and antimicrobial characterization of Macleaya cordata herb, Fitoterapia 81(8)
1006-1012.
[9] K. Pěn?Íková, J. Urbanová, P. Musil, E. Táborská, J. Gregorová, Seasonal Variation of Bioactive Alkaloid Contents in Macleaya microcarpa (Maxim.) Fedde, Molecules 16(12)
3391-3401.
472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515
[10] Z.X. Qing, P. Cheng, X.B. Liu, Y.S. Liu, J.G. Zeng, W. Wang, Structural speculation and identification of alkaloids in Macleaya cordata fruits by high‐performance liquid chromatography/quadrupole‐time‐of‐flight mass spectrometry combined with a screening procedure, Rapid Communications in Mass Spectrometry 28(9) (2014) 1033-1044. [11] L. Yu, C. Pi, T. Qi, Z. Li, W. Wei, Z. Tang, J. Zeng, Changes of Alkaloids Content and Expression of Related Genes for Macleaya by Treating for Different Stress Conditions, Molecular Plant Breeding (2018). [12] J.Í. Vrba, P. Dole?el, J. Vi?ar, M. Modriansky, J. Ulrichová, Chelerythrine and dihydrochelerythrine induce G1 phase arrest and bimodal cell death in human leukemia HL-60 cells, 22(4)
0-1017.
[13] K. Rogers, X. Chen, Biogenesis, Turnover, and Mode of Action of Plant MicroRNAs, Plant Cell 25(7)
2383-2399.
[14] X. Chen, Small RNAs–secrets and surprises of the genome, The plant journal 61(6) (2010) 941-958. [15] X. Chen, Small RNAs in development – insights from plants, Current Opinion in Genetics & Development 22(4)
0-0.
[16] S. Varsha, K. Vinay, D.R. M., K.T. S., W.S. H., MicroRNAs As Potential Targets for Abiotic Stress Tolerance in Plants, Front Plant Sci 7. [17] X. Wu, P.A. Northcott, A. Dubuc, A.J. Dupuy, D.J.H. Shih, H. Witt, S. Croul, E. Bouffet, D.W. Fults, C.G. Eberhart, Clonal selection drives genetic divergence of metastatic medulloblastoma, Nature 482(7386)
529-533.
[18] P.P. Gardner, D. Jennifer, J.G. Tate, E.P. Nawrocki, D.L. Kolbe, L. Stinus, A.C. Wilkinson, R.D. Finn, G.J. Sam, S.R. Eddy, Rfam: updates to the RNA families database, Nucleic Acids Research (suppl_1) (2008) suppl_1. [19] H. Li, R. Durbin, Fast and accurate short read alignment with Burrows-Wheeler transform, Bioinformatics
(2009).
[20] K. Ana, G.J. Sam, miRBase: annotating high confidence microRNAs using deep sequencing data, Nucleic Acids Research (D1) (2013) D1. [21] L.-C. Wan, F. Wang, X. Guo, S. Lu, Z. Qiu, Y. Zhao, H. Zhang, J. Lin, Identification and characterization of small non-coding RNAs from Chinese fir by high throughput sequencing, Bmc Plant Biol 12(1)
146.
[22] M.R. Friedlander, S.D. Mackowiak, N. Li, W. Chen, N. Rajewsky, miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades, Nucleic Acids Research 40(1) 37-52. [23] M.I. Love, W. Huber, S. Anders, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2, Genome Biol 15(12)
550.
[24] H.J. Wu, Y.-K. Ma, T. Chen, M. Wang, X.-J. Wang, PsRobot: a web-based plant small RNA meta-analysis toolbox, Nucleic Acids Research 40(W1)
W22-W28.
[25] M. Ashburner, C.A. Ball, J.A. Blake, D. Botstein, J.M. Cherry, Gene ontology: Tool for the unification of biology, Nat Genet 25(1) (2000) 25-29. [26] D.V. Klopfenstein, Z. Liangsheng, P.B. S., R. Fidel, W.V. Alex, N. Aurélien, M.C. J., Y.J. M., B. Olga, W. Mark, GOATOOLS: A Python library for Gene Ontology analyses, Sci Rep-Uk 8(1)
10872-.
[27] K. Minoru, G. Susumu, KEGG: Kyoto Encyclopedia of Genes and Genomes, Nucleic Acids Research (1) (2000) 1.
516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 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
[28] C. Xie, X. Mao, J. Huang, Y. Ding, J. Wu, S. Dong, L. Kong, G. Gao, C.Y. Li, L. Wei, KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases, Nucleic Acids Research (suppl_2) (2011) suppl_2. [29] J. Zhang, H. Zhang, A.K. Srivastava, Y. Pan, J. Bai, J. Fang, H. Shi, J.-K. Zhu, Knock-down of rice microRNA166 confers drought resistance by causing leaf rolling and altering stem xylem development, Plant Physiology
pp.01432.2017.
[30] W.X. Li, Y. Oono, J. Zhu, X.-J. He, J.-M. Wu, K. Iida, X.-Y. Lu, X. Cui, H. Jin, J.-K. Zhu, The Arabidopsis NFYA5 Transcription Factor Is Regulated Transcriptionally and Posttranscriptionally to Promote Drought Resistance, Plant Cell 20(8)
2238-2251.
[31] M. Zhou, H. Luo, Role of microRNA319 in creeping bentgrass salinity and drought stress response, Plant signaling & behavior 9(4)
e28700.
[32] L. Zhai, Z. Liu, X. Zou, Y. Jiang, F. Qiu, Y. Zheng, Z. Zhang, Genome-wide identification and analysis of microRNA responding to long-term waterlogging in crown roots of maize seedlings, 61(15) (2010) 4157-4168. [33] Sunkar, R., Novel and Stress-Regulated MicroRNAs and Other Small RNAs from Arabidopsis, Plant Cell 16(8)
2001-2019.
[34] H.H. Liu, X. Tian, Y.-J. Li, C.-A. Wu, C.-C. Zheng, Microarray-based analysis of stress-regulated microRNAs in Arabidopsis thaliana, Rna 14(5)
836-843.
[35] T.P. Frazier, G. Sun, C.E. Burklew, B. Zhang, Salt and Drought Stresses Induce the Aberrant Expression of microRNA Genes in Tobacco, 49(2)
159-165.
[36] Y.S. Dongmei Liu, Zhixiang Chen*, Diqiu Yu*, Ectopic expression of miR396 suppresses GRF target gene expression and alters leaf growth in Arabidopsis, Physiol Plantarum 136(2) (2009) 223-236. [37] W. Liu, Z. Yonggang, L. Xiaowei, W. Xingchao, D. Yuanyuan, W. Nan, L. Xiuming, C. Huan, Y. Na, C. Xiyan, Tissue-Specific Regulation of Gma-miR396 Family on Coordinating Development and Low Water Availability Responses, Front Plant Sci 8
1112-.
[38] Y. Ren, L. Chen, Y. Zhang, X. Kang, Z. Zhang, Y. Wang, Identification of novel and conservedPopulus tomentosamicroRNA as components of a response to water stress, Funct Integr Genomics 12(2)
327-339.
[39] V. Eldem, U.C. Akcay, E. Ozhuner, Y. Bakır, S. Uranbey, T. Unver, Genome-wide identification of miRNAs responsive to drought in peach (Prunus persica) by high-throughput deep sequencing, Plos One 7(12) (2012) e50298. [40] M. Arshad, B.A. Feyissa, L. Amyot, B. Aung, A. Hannoufa, MicroRNA156 improves drought stress tolerance in alfalfa ( Medicago sativa ) by silencing SPL13, Plant Science 258
122-136.
[41] M.K. Rezaei, Z.-S. Shobbar, M. Shahbazi, R. Abedini, S. Zare, Glutathione S-transferase (GST) family in barley: Identification of members, enzyme activity, and gene expression pattern, J Plant Physiol 170(14)
1277-1284.
[42] R. Sharma, A. Sahoo, R. Devendran, M. Jain, Over-Expression of a Rice Tau Class Glutathione S-Transferase Gene Improves Tolerance to Salinity and Oxidative Stresses in Arabidopsis, Plos One 9 (2014). [43] X. Jing, X. Xiao-Juan, T. Yong-Sheng, P. Ri-He, X. Yong, Z. Wei, Y. Quan-Hong, Z. Hong, Transgenic Arabidopsis Plants Expressing Tomato Glutathione S-Transferase Showed Enhanced Resistance to Salt and Drought Stress, Plos One 10(9)
e0136960-.
558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582
[44] G. Li, X. Peng, L. Wei, G. Kang, Salicylic acid increases the contents of glutathione and ascorbate
583
Tables list
and temporally regulates the related gene expression in salt-stressed wheat seedlings, Gene 529(2) 321-325. [45] D.C. Rees, E. Johnson, O. Lewinson, ABC transporters: the power to change, Nat Rev Mol Cell Biol 10(3)
218-227.
[46] L. Yang, Y. Jin, W. Huang, Q. Sun, F. Liu, X. Huang, Full-length transcriptome sequences of ephemeral plant Arabidopsis pumila provides insight into gene expression dynamics during continuous salt stress, BMC genomics 19(1) (2018) 717. [47] S.J. Suh, Y.-F. Wang, A. Frelet, N. Leonhardt, M. Klein, C. Forestier, B. Mueller-Roeber, M.H. Cho, E. Martinoia, J.I. Schroeder, The ATP Binding Cassette Transporter AtMRP5 Modulates Anion and Calcium Channel Activities in Arabidopsis Guard Cells, J Biol Chem 282(3)
1916-1924.
[48] M. Lee, Y. Choi, B. Burla, Y.-Y. Kim, B. Jeon, M. Maeshima, J.-Y. Yoo, E. Martinoia, Y. Lee, The ABC transporter AtABCB14 is a malate importer and modulates stomatal response to CO2, Nature Cell Biology 10(10)
1217-1223.
[49] M. Jain, J.P. Khurana, Transcript profiling reveals diverse roles of auxin-responsive genes during reproductive development and abiotic stress in rice , Febs J 276(11) (2009) 3148-3162. [50] R. Yue, S. Tie, T. Sun, L. Zhang, Y. Yang, J. Qi, S. Yan, X. Han, H. Wang, C. Shen, Genome-wide identification and expression profiling analysis of ZmPIN, ZmPILS, ZmLAX and ZmABCB auxin transporter gene families in maize (Zea mays L.) under various abiotic stresses, Plos One 10(3) (2015) e0118751. [51] C. Chai, S.P. K., Comprehensive Analysis and Expression Profiling of the OsLAX and OsABCB Auxin Transporter Gene Families in Rice (Oryza sativa) under Phytohormone Stimuli and Abiotic Stresses, Front Plant Sci 7.
Table 1. Summary of the RNA sequencing results
584 Raw reads Clean reads miRNA
Total_num
s_M_B_1
s_M_B_2
s_M_B_3
s_M_T_1
s_M_T_2
s_M_7_3
101375907 80957233( 79.86%) 18931949 23.39%
18045994 14627819 (81.06%) 4282093 29.27%
15684160 11814745 (75.32%) 3159378 26.74%
19063526 16355079 (85.79%) 3654379 22.34%
15020934 11799819 (78.56%) 2558940 21.69%
17026663 14451133 (84.87%) 3218371 22.27%
16534630 11908638 (72.02%) 2058788 17.29%
585 586 587
Table 2. Primer sequences
Genes 18S -F 18S -R cca-miR396a-3p-RT cca-miR396a-3p-F cca-miR396a-5p-RT cca-miR396a-5p-F
Primer sequence CTTCGGGATCGGAGTAATGA GCGGAGTCCTAGAAGCAACA CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGTTTCCCA ACACTCCAGCTGGGGTTCAAGAAAGCTG CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGCAGTTC ACACTCCAGCTGGGTTCCACAGCTTTCTT
cca-miR395b-RT cca-miR395b-F Nov-m0411-3p-RT Nov-m0411-3p-F Nov-m0067-5p-RT Nov-m0067-5p-F All R ABCB1-F ABCB1-R GST-F GST-R PHOT2-F PHOT2-R MPC1-F MPC1-R AT4G37250-F AT4G37250-R ADH5-F ADH5-R APS1-F APS1-R PAP1-F PAP1-R SPL6-F SPL6-R SPL15-F SPL15-R GRF3-F GRF3-R GRF6-F GRF6-R ATHB-15-F ATHB-15-R REVT-F REVT-R SCL6-F SCL6-R SCL15-F SCL15-R
CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGGAGTTCC ACACTCCAGCTGGGCTGAAGTGTTTGGA CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGAGAGCCT ACACTCCAGCTGGG CAGGCTGTGTGTCC CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGAATCTCT ACACTCCAGCTGGGTTTTTGTCGGGACC TGGTGTCGTGGAGTCG AACCCACTAAAAAGCAAACGATTC GAGGAGCACCACCGTCATTAT CCGGTTCACAAGAAGATTCCG GTTCGACATTAAGAATTAACACCCA GCAGCTTTATGCAGGCTCAC CCTGAAAGGTGTGCGACCA GTGACTGAAGGGCGACTGAA CAGTCGAGGGCAGTGATCTG TCCTTTTCCTGCCCAAGGAC CTTATCCCTTCCTCGCCGTC TTTGCCACACCGATGCCTAT ATCTCGGACTAGGAAGGCAGA GATCGTTTCAGGCTTTCGCC GCCTTTGAAGGCCCATTACAAC CTTGACGTCATTTGCGGTCC TAGCAAATTAAGCACGACAGC CAAAAGAAGTTGCAGAAGGCGA CGCTTACCCGGTAGGAAATAGAA CACCCATTGTGGTGCCATTG GCAAGATCAGTGATGCAGCG GGTGGAGAAGTTCAGATGAAGAGT GTGGTGGTGGTGGTAGTGAG TGGCTACAAGGATTCCATTCACA TGCAGCTAAGTACTTGAACACAA TGACGGCGAATAGAACTCGG CTGGAGGGGGTAGTAGGGTC AGGCGTTGCGGTATATCAGG TCCAACGGCATTCTCGTCAA ATCATCGCCTCCCCTCTCTT ATGCGGTCACCTTGAGGTTC ATGTGATACTGGCACGGCTC GCGTTGGCTTTACGACACTC
588 589
Figure titles list:
590 591
Fig.1. Expression profile of all detected miRNAs Fig.2. Cluster analysis of the differentially expressed known miRNAs
592 593 594 595 596
Fig.3. Cluster analysis of the differentially expressed novel RNAs Fig.4. Enrichment of target Genes in GO pathway Fig.5. Enrichment of target Genes in KEGG pathway Fig.6. Expression of miRNAs and the target genes
S0003-9861(19)30922-1
DOI:
https://doi.org/10.1016/j.abb.2020.108300
Reference:
YABBI 108300
To appear in:
Archives of Biochemistry and Biophysics
Received Date: 14 October 2019 Revised Date:
15 January 2020
Accepted Date: 8 February 2020
Please cite this article as: L. Yu, L. Zhou, W. Liu, P. Huang, R. Jiang, Z. Tang, P. Cheng, J. Zeng, Identification of drought resistant miRNA in Macleaya cordata by high-throughput sequencing, Archives of Biochemistry and Biophysics (2020), doi: https://doi.org/10.1016/j.abb.2020.108300. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Inc.
1
Identification of drought resistant miRNA in Macleaya
2
cordata by high-throughput sequencing
3 4 5 6 7 8 9 10 11 12 13 14 15 16
Linlan Yua 2 Li Zhoua 2 Wei Liua 1,3 Peng Huang1,2 Ruolan Jiang1,2 Zhaoshan Tang5 Pi Cheng*1,2
17 18
Wei Liu#: [email protected] Peng Huang: [email protected]
19 20 21 22 23 24 25
Ruolan Jiang: [email protected]
Jianguo Zeng*1,4 1
Hunan Key Laboratory of Traditional Chinese Veterinary Medicine, Hunan Agricultural University,
Changsha, 410128 Hunan China 2
College of Horticulture and Landscape, Hunan Agricultural University, Changsha, 410128 Hunan
China 3
Center of Analytic Service, Hunan Agriculture University, 410208 Changsha, China
4
National and Local Union Engineering Research Center of Veterinary Herbal Medicine Resource and
Initiative, Hunan Agricultural University, Changsha, 410128 Hunan China 5
Micolta Bioresource Inc., Changsha, 410016 China
Email Linlan Yu#: [email protected] Li Zhou#: [email protected]
Zhaoshan Tang: [email protected] Pi Cheng*: [email protected]
Jianguo Zeng*: [email protected] * Corresponding author a These authors contributed equally to this work.
26
Abstract:
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
Drought is one of the most serious factors affecting crop yields in the world. Macleaya cordata (Willd.) is a draught-tolerant medicinal plant that has been proposed as a pioneer crop to be cultivated in arid areas. However, the exact molecular mechanisms through which M. cordata responds to draught stress remain elusive. In recent years, microRNA (miRNAs) in plants have been associated with stress response. Based on these findings, the current study aimed to shed light on the potential regulatory roles of miRNAs in the draught tolerance of M. cordata by employing high-throughput RNA sequencing and degradation sequencing. Six M. cordata plants were randomly divided into two equal experiment groups, including one draught group and one control group. High-throughput sequencing of the M. cordata samples led to the identification of 895 miRNAs, of which 18 showed significantly different expression levels between the two groups. PsRobot analysis and degradation sequencing predicted the differential miRNAs to target 59 and 36 genes, respectively. Functional analysis showed that 38 of the predicted genes could be implicated in the modulation of stress response. Four miRNAs and eight target genes were selected for quantitative real-time polymerase chain reaction (qRT-PCR) validation. The expression trend of each miRNA analyzed by qRT-PCR was consistent with that determined by sequencing, and was negatively correlated with those of its target genes.The results of our current study supported the involvement of miRNAs in the draught tolerance of M. cordata and could pave the way for further investigation into the related regultory mechanisms.
48
Keywords: Macleaya cordata
49
sequencing
50
miRNAs
degradation sequencing
drought
stress
high-throughput
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
Background Drought is one of most serious type of abiotic stress for crops worldwide [1]. Studies have shown that draught conditions can cause the temperature of plant cells to increase and reactive oxygen species (ROS) to accumulate, leading to disrupted cell metabolism and even apoptosis [2]. Furthermore, extreme draught stress can also trigger the closure of stomata, which significantly hampers photosynthesis [3]. Not surprisingly, plants have evolved a wide range of molecular and physiological mechanisms to help them resist draught stress and mitigate the detrimental effects of water deficit [4]. For example, studies on Arabidopsis thaliana have revealed that a number of kinases, including SPK1 and MAPK, are activated by draught-induced up-regulation of abscisic acid (ABA) and can subsequently promote physiological changes designed to prevent excessive water loss [5, 6]. Additionally, there is also mounting evidence that shifts in the biosynthesis of certain cell wall components, such as cellulose and xyloglucan, constitute another strategy that plants employ in response to draught [7]. Better understanding of these response mechanisms would greatly facilitate the development of novel draught-tolerant crop cultivars with enhanced economic value. Macleaya cordata (Willd.) R. Br. is a herbaceous perennial mainly distributed across Southeast China and Southeast Asia [8, 9]. Extracts prepared from M. cordata are an important ingredient in traditional Chinese medicine with demonstrated antimicrobial, insecticidal, antitumoral and anti-inflammatory effects [10-12]. Based on chemical analysis, the main bioactive constituents of the plant consist of isoquinoline alkaloids such as sanguinarine, chelerythrine, protopine and allocryptopine. Interestingly, M. cordata has been shown to exhibit moderate draught tolerance and thus proposed as a pioneer plant species for soil conditioning in arid lands. We have recently discovered that M. cordata exposed to simulated draught stress with 30% PEG6000 showed no observable decrease in growth rate and produced 170% more sanguinarine over 48 h compared to the controls cultivated under regular conditions [11]. Despite these results, the draught response mechanisms in M. cordata remain elusive and require further studies to decipher. MicroRNAs (miRNA) are a diverse family of endogenous single-stranded RNA molecules consisting of 21 to 24 nucleotides [13]. Although miRNA sequences do not encode proteins themselves, they have been increasingly recognized as critical regulators of gene expression [14]. In general, miRNAs are capable of inhibiting the translation, or even directly inducing the degradation, of the target mRNAs by fully or partially complementing with their 3’ untranslated regions (UTR) [15]. Recently, miRNAs have also been shown to participate in plant response to various abiotic stress signals, such as draught, heat, oxidative stress and heavy-metal toxicity [16]. The rapid advances of high-throughput sequencing technologies and bioinformatic algorithms have provided a boon for miRNA research, enabling the systematic examination of miRNAs expression profiles in the same host under different environmental or physiological conditions. Combined, these developments have served as the foundation for our current study, in which we sought to identify miRNAs associated with the draught tolerant traits of M. cordata. Comparison of the miRNAs expression profiles between the drought-stressed treatment group and the control group allowed us to identify 18 differentially expressed candidates. Further functional analysis suggested that the differentially expressed miRNAs were involved in biological functions such as carboxylic ester hydrolase activity, phosphorus metabolic process, membrane. Together,
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138
these findings could increase our understanding of the complex regulatory network that underlies the draught tolerance of M. cordata, and provide a theory basis for encouraging its cultivation in arid regions for vegetation restoration.
Methods Draught induction M. cordata were collected from Huaihua, China with no permission required, and cultivated in Hunan Agricultural University. Six M. cordata plantlets derived from the same explant were grown at 25 °C in Murashige-Skoog (MS) medium (4.42 g/L) supplemented with 30 g/L sucrose and 3 g/L Phytagel at a pH of 5.8 under a 16/8-h light/dark photoperiod for 30 days. The plantlets were randomly divided into two equal groups, including a drought-stressed treatment group and a control group. The three plantlets in the draught group were subjected to draught stress by transfer to a new MS medium containing 30% (w/v) PEG6000, whereas the control group were cultivated in fresh MS medium. After 3 days, the plantlets were harvested, ground to a fine powder in liquid nitrogen, and then stored at -80 °C.
Library construction Total RNA was extracted from 100 mg of ground plant materials using Plant RNA Purification Kit (Thermo Fisher, USA) according the manufacturer’s instructions, followed by the removal of genomic DNA with RNase-free recombinant DNase I (TaKaRa, Japan). The quality of the extracted RNA was verified on a 2100 Bioanalyzer (Agilent Technologies, USA) and a NanoDrop 2000 spectrometer (Thermo Fisher Scientific, USA). Only high-quality RNA, defined as meeting the criteria of OD260/280 = 1.8 ~ 2.2, OD 260/230 ≥ 2.0, RIN ≥ 6.5 and 28S:18S ≥ 1.0, was used to construct a sequencing library.
Library preparation, and Illumina Hiseq xten Sequencing Stranded RNA-seq library for transcriptome analysis was prepared from 5 μg of high-quality total RNA using a Sample Prep Kit (Illumina, CA). First-stranded cDNA was synthesized with random hexamer primers, followed by template removal and replacement strand synthesis with the incorporation of dUTP in place of dTTP. As the polymerase could not utilize dUTP, the elongation of the complementary strand was aborted prematurely, leading to incomplete blunt-ended double-stranded cDNA fragments of varying lengths. These fragments were subsequently purified by AMPure XP beads. A single “A” nucleotide was added to each 3’ end of these fragments to prevent blunt end ligation. The modified cDNA fragments were then ligated with multiple indexing adapters and ollowed by PCR amplification with Phusion DNA polymerase (New England Biolabs, USA) for 15 cycles. The resultant cDNA library was quantified on a TBS-380 Mini-fluorometer (Turner Biosystems, USA) and sequenced on a HiSeq X Ten system (2 × 150bp reads, Illumina, USA).
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
Identification and expression analysis of identified miRNAs The 3’ ends of all low-quality bases with a quality score below 20 were trimmed using in-house PERL scripts, and the adapters were then removed with FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) [17]. Among the remaining sequences, those with 18 to 32 nucleotides were retained. The assembled unique sequences were used as queries for a BLAST search against the Rfam database (version 12.1) [18] to remove non-miRNA members such as rRNA, tRNA and snoRNA. Bowtie [19] (version 1.2.2) was used to map the predicted miRNAs)against the reference genome and predict their secondary structures. On the one hand, known miRNAs were identified by searching miRbase [20] (version 21.0, http://www.mirbase.org/) and retaining the perfectly matched sequences. On the other hand, novel candidates were predicted by MIREAP [21] based on the characteristic hairpin structure of miRNA precursors. In addition, in-house scripts were used to determine the nucleotide bias at each position of the identified miRNAs. The expression level of each miRNA was calculated by miRDeep2 [22] (version 2.0.0.8) in the form of transcript per million reads (TPM). Differential expression analysis was conducted by using the DESeq2 package [23] (version 1.6.3, http://www.bioconductor.org/packages/release/bioc/html/DESeq2.htmls) based on the criteria of |log2(fold change)| > 1 and false discovery rate-adjusted p value < 0.05.
Annotation and functional analysis The target genes of all identified mRNAs were predicted and annotated with psRobot [24]. Gene ontology (GO) analysis was performed by querying the genes against the GO database (http://www.geneontology.org/) [25] using GOATOOLS (version 0.5.9, https://github.com/tanghaibao/GOatools) [26]. GOATOOLS employed Fisher’s exact test to compute the p value that indicates whether the frequency of the target genes assigned to a specific GO term is significantly different from that of randomly selected background genes. The p value was subsequently corrected using four different methods, including Bonferroni, Holm, Sidak and False Discovery Rate. Meanwhile, Kyoto Encyclopedia of Genes and Genomes (KEGG) [27] pathway analysis was performed by querying the genes against the KEGG database (http://www.genome.jp/kegg/) with KOBAS [28]. P value was calculated by applying Fisher’s exact test and then corrected by the Benjamini-Hochberg method.
Degradation sequencing The mRNA targets of the miRNAs were verified by degradation sequencing. After binding with miRNA, the mRNA is cleaved and the resultant 3'-fragment, which contains a free 5'-monophosphate and a 3' polyA tail, can be extended by RNA ligase and subsequently sequenced. In-depth alignment analysis of the sequenced data can then reveal the candidate miRNA splicing sites in the mRNA sequence. Briefly, mRNA molecules were enriched from total
183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
RNA by Total RNA Purification Kit (LC Science, Houston, USA) and then subjected to 5' adaptor ligation, followed by reverse transcription with mixed (PrimeScript®1st Strand cDNA Synthesis Kit , TaKaRa, Japan) and (PCR Ampification Kit Manual, Takara, Japan) according to the manufacturer’s instructions. The generated library was sequenced on a Hiseq 2500 system (Illumina, USA).
Quantitative real-time PCR (qRT-PCR) validation Selected miRNA candidates were reverse-transcribed with RevertAid Premium Reverse Transcriptase (Thermo Fisher Scientific, USA) following the manufacturer’s instructions. The synthesized cDNAs were diluted ten-fold and used directly as template. Each reaction consisted 10 µM each primer and 1 μg of template in 1 × Fast SG qPCR Master Mix (Sangon Biotech, China). The qRT-PCR amplification was performed on a StepOnePlus Real Time PCR System (Applied Biosystems, USA) programmed as follows: 95 °C for 3 min, followed by 45 cycles of 95 °C for 3 s and 60 °C for 30 s. The target mRNAs were reverse-transcribed into cDNAs with Prime Script RT Reagent Kit with gDNA Eraser (Perfect Real Time, TaKaRa, Japan) following the manufacturer’s instructions. Each reaction consisted of 3 µM of each primer and 1 μg of template in 1× FastStart Universal SYBR Green Master (Roche, USA). The qRT-PCR amplification was performed on a 7300 Fluorescence Quantitative PCR Analyzer (Applied Biosystems, USA) programmed as follows: 95°C for 5 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 60 s. All PCR primers used in this study were summarized in Table 2.
Results
Overview of the high-throughput sequencing results Our high-throughput sequencing generated a total of 101375907 raw reads from the six libraries. All low-quality reads were then removed, and only reads consisting of 18 – 32 nucleotides were retained. Eventually, we obtained 80957233 clean reads. Approximately 44.33% of these reads were successfully matched to Rfam. The matching suggested that roughly 23.39% of the clean reads were likely derived from miRNA sequences, with the rest comprising rRNAs, tRNAs, snRNAs and other unknown RNA molecules. The percentage of miRNAs in the clean reads derived from each sequencing library ranged from 17.3% to 29.3%. Detailed descriptions of the various types of RNAs that we identified from the M. cordata samples were provided in Table 1 and additional file 1.
Identification of miRNAs
227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270
Comparison with miRBase-deposited sequences allowed us to identify 355 miRNA candidates belonging to 68 miRNA families. The expression level of each miRNA was then estimated by miRDeep2. Notably, the 13 most highly expressed miRNAs from each sample were all from the miR166 family. Meanwhile, we also predicted an additional 540 novel miRNAs belonging to 334 families. These novel candidates generally showed much lower expression levels than the known miRNAs described above. It is worth emphasizing that the miRNA expression profiles varied significantly among different samples. A comparative analysis indicated that 337 known and 294 predicted miRNA candidates were detected in all six M. cordata samples. In contrast, 36 known and 112 novel miRNAs were found only in the draught group, whereas 42 known and 136 novel candidates were exclusively detected in the control samples.
Differential expression analysis We then compared the miRNA expression profiles between the draught group and the control group. A miRNA is considered differentially expressed if there was a two-fold difference in its expression between the two experiment groups and if the corresponding FDR-adjusted p value was below 0.05. Based on these criteria, nine known and nine novel miRNAs exhibited differential expression between the draught group and the control group. Among them, six known (ath-miR8175, ath-miR827, cca-miR395b, gma-miR395d, gma-miR395i, mtr-miR395h) and five novel (Nov-m0067-5p, Nov-m0176-5p, Nov-m0231-5p, Nov-m0383-5p, Nov-m0396-3p) miRNAs were up-regulated in the draught group, whereas the rest were down-regulated (known miRNAs: cca-miR396a-3p, cca-miR396c and stu-miR156d-3p; novel miRNAs: Nov-m0038-3p, Nov-m0112-5p, Nov-m0219-5p and Nov-m0411-3p). The expression of all detected miRNAs was summarized in figure 1 and 2. Cluster analysis also provided strong evidence that there were clear distinctions between the miRNA expression patterns of the two experiment groups (Fig. 3 and 4).
Target prediction and functional analysis of the drought-responsive miRNAs
The target genes of the identified differentially expressed miRNAs were predicted by psRobot [24]. The nine known miRNAs were predicted to interact with a total of 11 target genes. Target multiplicity, a hallmark of miRNA-mRNA interaction, was unambiguously demonstrated by the observation that most of the miRNAs were predicted to complement the 3’ UTR of more than one mRNA and nearly half of the identified genes were recognized by two or more miRNAs. Similarly, the nine novel miRNA candidates were shown to target 48 genes. In particular, Nov-m0067-5p was predicted to interact with 40 different mRNAs. Nov-m0411-3p could potentially complement the 3’-UTRs of 6 mRNAs, whereas Nov-m0219-5p and Nov-m0396-3p were found to each target one gene. In sharp contrast to what was observed of the known miRNAs, none of the identified genes was predicted to be the target of more than one of the differentially expressed novel miRNA candidates. Known miRNAs such as miR396, miR395,
271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314
miR156-3 and miR827 were predicted to respond to drought by regulating the target genes UBC, ADH5, RPL8 and ALDH. Meanwhile, Nov-miR-0067, Nov-miR-411, Nov-m0383-5p and Nov-m0396-3p were predicted to respond to drought stress by regulating the target genes PLD, GST, ABCB1, SOD2, and APRT.
Annotation and functional analysis of differentially expressed genes Annotation of the combined 59 target genes revealed that they encode a wide assortment of growth- and/or stress-related proteins, such as phospholipase D, ABC transporter, GST, phosphotransferase, kinase, transcription factor, GDSL esterase, zinc finger protein, and pherophorin. The annotated genes were subsequently subjected to GO enrichment and KEGG pathway analyses. The GO annotations of the genes were predicted by GOATOOLS (version: 0.5.9, https://github.com/tanghaibao/GoOatools) and classified into three sub-categories, including molecular function, biological process or cellular component. Based on these results, the GO terms that could be meaningfully interpreted and associated with stress response and/or draught tolerance included hydrolase activity, acting on ester bonds (GO: 0016788), hydrolase activity (GO:0016787) , primary metabolic process (GO:0044238), organic substance metabolic process (GO:0071704), carboxylic ester hydrolase activity (GO:0052689), small molecule binding (GO:0036094), , phosphate-containing compound metabolic process (GO:0006796) and phosphorus metabolic process (GO:0006793). On the other hand, KEGG analysis indicated that the putative target genes of the differentially expressed miRNAs were significantly enriched in a number of pathways that could contribute to plant stress response and/or draught tolerance, including GnRH signaling pathway (ko04912), ABC transporters (ko02010), Ras signaling pathway (ko04014), Glutathione metabolism (ko00480), Plant hormone signal transduction (ko04075) and Phospholipase D signaling pathway (ko04072).
Target validation by degradation sequencing In total, we predicted that the 895 miRNAs that we detected in the M. cordata samples could target a total of 22691 genes. We then sought to systematically verify these miRNA-mRNA interactions through degradation sequencing, which can selectively amplify the 3’ cleavage fragments of miRNAs. The nucleotide sequences of these fragments were then compared against the corresponding miRNA sequences. In case of perfect or near-perfect complementation, the cleaved mRNA would be considered as a target of the corresponding miRNA. Based on these principles, we identified 8022 mRNAs that produced 3’ fragments with the correct cleavage pattern indicative of miRNA binding. These mRNAs were collectively targeted by a combined 550 miRNAs. Importantly, 6310 genes were predicted only by psRobot, whereas 8022 were captured only by degradation sequencing. In the degradation group sequencing, multiple fragments of SPL (SQUAMOSA promoter binding like), GRF (growth regulating factor), and SCL (Scarecrow-like) that were cut by miR156 / 157, miR396 and miR170 / 171 were detected. SPL, GRF and SCL are the star molecules for drought resistance. In addition, KEGG function enrichment was performed on the detected target genes, among which the ratio related to adaptation environment was the
315 316 317 318 319 320 321 322
highest, KEGG function enrichment was performed on the detected target genes, among which the ratio related to adaptation to the environment was the highest, and the signal transduction was also a high proportion. These pathways are closely related to resistance to stress.
323 324 325
miR156d, aqc-miR166a and aqc-miR171a) and sixteen of their target genes (GST, ABCB1, PHOT2, MPC1, At4g37250, ADH5, APS1, PAP1, GRF3, GRF6, SPL6, SPL7, ATHB-15, REV, SCL6 and SCL15) that have previously been associated with draught tolerance for qRT-PCR validation.
326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358
Among them, four miRNAs (cca-miR395b, aqc-miR166a,aqc-miR171a and Nov-m0067-5p) were up-regulated in drought stress samples, and four (cca-miR396a-3p, cca-miR396a-5p, miR156d and Nov-m0411-3p) showed opposite trends. It is worth mentioning that seven of the chosen target genes, including GRF3, GRF6, SPL6, SPL7, At4g37250, PHOT2 and MPC1, were detected by degradation sequencing but not predicted by bioinformatic analysis. The results indicated that the expression patterns of all six miRNAs determined by qRT-PCR were consistent with those shown by RNA sequencing (Fig.5.). In addition, the expression trends of all target genes were found to be inversely correlated with those of their putative miRNA regulators. Among the target genes that we verified, PHOT2 showed extremely significantly augmented expression levels, GST, ABCB1, GRF3, GRF6, SPL6, SPL7, and MPC1 (shaggy-related protein kinase eta) showed significant increased, whereas At4g37250, ADH5, APS1 and PAP1 were extremely significantly down-regulated in the draught group, the rest down-regulated obviously. These results demonstrated that the RNA sequencing data that we obtained were sufficiently reliable and accurate for functional analysis and mechanistic interpretation.
qRT-PCR validation We selected four differentially expressed miRNAs (Nov-m0411-3p, cca-miR396a-3p, Nov-m0067-5p and cca-miR395b), four miRNAs known to regulate drought (cca-miR396a-5p,
Discussion MiRNAs have emerged as key regulators of abiotic stress in plants. As an example, knockdown of miRNA-166 or stimulation of its target gene homeodomain containing protein 4 (OsHB4) was found to improve draught resistance in rice by inducing the development of smaller bulliform cells and abnormal sclerenchymatous cells [29]. In another study, up-regulation of miR-169a or suppression of its target gene nuclear factor -YA (NFYA5) in Arabidopsis thaliana aggravated water loss through the stroma [30]. Overexpression of miR-319 and Osa-miR-319a in Agrostis stolonifera has been shown to enhance its tolerance for drought and salt stress [31]. These findings prompted researchers to speculate the possibility of manipulating miRNA expression to engineer plant cultivars with enhanced stress resistance. A number of differentially expressed miRNAs that we identified have also been shown in previous studies to be potentially implicated in draught resistance. However, there is clear evidence that the expression trends and the underlying regulatory mechanisms of these miRNAs could vary significantly among different plant species. For example, miR-396 was found to be down-regulated in both M. cordata, as evidenced in our current study, and Oryza sativa [32], in
359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402
response to water deficit, but exhibited enhanced expression in draught-stressed Arabidopsis thaliana [33]. On the other hand, miR-396 has been reported to regulate stromatal density and plant transpiration in Arabidopsis thaliana and tobacco by targeting several growth regulating factor genes [34, 35]. Arabidopsis thaliana overexpressing miR-396 was characterized by linear narrow leaves with reduced stomatal density that showed greater resistance to water deficit [36]. However, the opposite effect was evident in Arabidopsis thaliana root [37]. Once again, our degradation sequencing suggested, GRFS was very likely due to the modulation of miR-396. Similarly, miR-156, another down-regulated miRNA candidate in our study, showed decreased expression in response to drought stress in Oryza sativa] [32] and Populus tomentosa [38], but the opposite pattern was observed in peach [39] and Arabidopsis thaliana [33]. Members from the SQUAMOSA promoter binding protein (SPL) family have been shown to be regulated by miR-156 [40]. SPLs are known to help plants combat a variety of abiotic stresses by modulating the anthocyanin synthesis pathway. Therefore, inhibition of miR-156 would de-repress SPLs and concomitantly enhance the draught tolerance of the host. This was consistent with our degradation sequencing data that SPL2, SPL3, SPL6, SPL7 and SPL12 were targeted by miR-156 in draught-stressed M. cordata samples. Based on the above, we speculate that miR396 and miR156 participate in drought stress response by cutting GRFs and SPLs. We also carefully investigated the predicted target genes of the novel miRNAs that exhibited differential expression between the draught group and the control group. For example, GST, a putative target of NOV-m0411-3p, was significantly up-regulated in draught-exposed M. cordata and is a well-established modulator of oxidative stress. GST was reported to contribute to the development of draught resistance in barley [41]. Consistently, overexpression of OsGSTU4 [42] and introduction of the tomato GST gene[43]] could both significantly enhance the draught tolerance of Arabidopsis thaliana. In the latter case, notably, higher GST level resulted in increased production of proline, malondialdehyde and antioxidative enzymes in the engineered plants [44]. The function of ATP-binding cassette transporters (ABC transporters) is transmembrane transport of ions, macromolecules and other substances through energy-driven [45]. ABCB1 was another target gene that we validated and functions as an auxin transporter. It can be activated with other proteins to participate in ion homeostasis in the cytoplasm in response to stress [46]. Studies have shown that the AtABCC5 mutant regulates ABA transport, regulates stomata, reduces transpiration rate, thereby reducing water consumption and enhancing drought tolerance [47]. AtABCB14 has similar functions [48]. As a prominent osmotic regulator, auxin plays a central role in plant response to abiotic stress signals such as water shortage and high salinity [49]. Therefore, it is unsurprising that auxin transporters have been associated with drought tolerance in both maize and Zea mays [50]. Notably, one of the plausible mechanisms through which auxin transporters might modulate draught resistance pathways was postulated to involve their stimulation of ABA production [51]. Taken together, our results suggested that, Nov-m0411-3p, which targets both GST and ABCB1, could be a key miRNA contributor to draught tolerance in M. cordata.
Conclusions
403 404 405 406 407 408 409 410 411 412 413 414 415 416
417 418 419 420 421 422 423 424 425 426 427
In summary, we provided the first systematic study that compared the miRNA expression profiles between draught-stressed M. cordata plants and controls cultivated under normal growth conditions. With a combination of high-throughput RNA sequencing and bioinformatic analysis, we identified a total of 895 miRNAs that showed differential expression as a result of induced water deficit. PsRobot-based bioinformatic analysis and degradation sequencing suggested that these miRNAs could target 6310 and 8022 genes, respectively. Based on functional analysis, these gene targets could be involved in a variety of pathways associated with plant response to draught, such as glutathione metabolism, plant hormone signal transduction and sulfur metabolism. Combined, our findings could facilitate the elucidation of molecular mechanisms that underlie the draught tolerance traits of M. cordata.
Abbreviations Full name
Abbreviations
microRNA 3’ untranslated regions Sphingosine kinase 1 Mitogen-activated protein kinase Abscisic acid ATP-binding cassette subfamily B member 1 Glutathione S-transferase Phototropin 2 Mitochondrial pyruvate carrier 1 Alcohol dehydrogenase 5 purple acid phosphatase 1 SQUAMOSA promoter binding like Growth regulating factor Arabidopsis thaliana homeobox genes REVOLUTA Scarecrow-like gene
miRNA UTR SPK1 MAPK ABA ABCB1 GST PHOT2 MPC1 ADH5 PAP 1 SPL GRF ATHB REVT SCL
Declarations Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Availability of data and materials
428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471
The data sets supporting the results of this article are included within the article and its additional files. Gut metagenome sequences were deposited in the National Center for Biotechnology Information (NCBI) SRP140359. Funding This work was supported by the innovation guidance program project of Hunan provincial science and Technology Department (2016SK3002), the National key laboratory cultivation base construction project (16KFXM09), and the Natural Science Foundation of Hunan Province (No. 2017JJ2115). Authors’ contributions L.Y., P. C. and J.Z. conceived and designed research. L.Y., L.Z., W.L. and P.H. onducted experiments. R.J. and Z.T. analyzed data. L.Y. wrote the manuscript. All authors read and approved the manuscript.
Competing Interests
The authors declare no competing interests.
References
[1] L. Cattivelli, F. Rizza, F.-W. Badeck, E. Mazzucotelli, A.M. Mastrangelo, E. Francia, C. Marè, A. Tondelli, A.M. Stanca, Drought tolerance improvement in crop plants: An integrated view from breeding to genomics, Field Crops Research 105(1-2)
0-14.
[2] M. Farooq, S. Basra, A. Wahid, Z. Cheema, M. Cheema, A. Khaliq, Physiological role of exogenously applied glycinebetaine to improve drought tolerance in fine grain aromatic rice (Oryza sativa L.), Journal of Agronomy and Crop Science 194(5) (2008) 325-333. [3] G. Cornic, Drought stress inhibits photosynthesis by decreasing stomatal aperture–not by affecting ATP synthesis, Trends in plant science 5(5) (2000) 187-188. [4] B. Valliyodan, H.T. Nguyen, Understanding regulatory networks and engineering for enhanced drought tolerance in plants, Curr Opin Plant Biol 9(2) (2006) 189-195. [5] Y. Osakabe, K. Maruyama, M. Seki, M. Satou, K. Shinozaki, K. Yamaguchi-Shinozaki, Leucine-rich repeat receptor-like kinase1 is a key membrane-bound regulator of abscisic acid early signaling in Arabidopsis, The Plant Cell 17(4) (2005) 1105-1119. [6] C. Jonak, S. Kiegerl, W. Ligterink, P.J. Barker, N.S. Huskisson, H. Hirt, Stress signaling in plants: a mitogen-activated protein kinase pathway is activated by cold and drought, Proceedings of the National Academy of Sciences 93(20) (1996) 11274-11279. [7] H. Le Gall, F. Philippe, J.-M. Domon, F. Gillet, J. Pelloux, C. Rayon, Cell wall metabolism in response to abiotic stress, Plants 4(1) (2015) 112-166. [8] P. Kosina, J. Gregorova, J. Gruz, J. Vacek, M. Kolar, M. Vogel, W. Roos, K. Naumann, V. Simanek, J. Ulrichova, Phytochemical and antimicrobial characterization of Macleaya cordata herb, Fitoterapia 81(8)
1006-1012.
[9] K. Pěn?Íková, J. Urbanová, P. Musil, E. Táborská, J. Gregorová, Seasonal Variation of Bioactive Alkaloid Contents in Macleaya microcarpa (Maxim.) Fedde, Molecules 16(12)
3391-3401.
472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515
[10] Z.X. Qing, P. Cheng, X.B. Liu, Y.S. Liu, J.G. Zeng, W. Wang, Structural speculation and identification of alkaloids in Macleaya cordata fruits by high‐performance liquid chromatography/quadrupole‐time‐of‐flight mass spectrometry combined with a screening procedure, Rapid Communications in Mass Spectrometry 28(9) (2014) 1033-1044. [11] L. Yu, C. Pi, T. Qi, Z. Li, W. Wei, Z. Tang, J. Zeng, Changes of Alkaloids Content and Expression of Related Genes for Macleaya by Treating for Different Stress Conditions, Molecular Plant Breeding (2018). [12] J.Í. Vrba, P. Dole?el, J. Vi?ar, M. Modriansky, J. Ulrichová, Chelerythrine and dihydrochelerythrine induce G1 phase arrest and bimodal cell death in human leukemia HL-60 cells, 22(4)
0-1017.
[13] K. Rogers, X. Chen, Biogenesis, Turnover, and Mode of Action of Plant MicroRNAs, Plant Cell 25(7)
2383-2399.
[14] X. Chen, Small RNAs–secrets and surprises of the genome, The plant journal 61(6) (2010) 941-958. [15] X. Chen, Small RNAs in development – insights from plants, Current Opinion in Genetics & Development 22(4)
0-0.
[16] S. Varsha, K. Vinay, D.R. M., K.T. S., W.S. H., MicroRNAs As Potential Targets for Abiotic Stress Tolerance in Plants, Front Plant Sci 7. [17] X. Wu, P.A. Northcott, A. Dubuc, A.J. Dupuy, D.J.H. Shih, H. Witt, S. Croul, E. Bouffet, D.W. Fults, C.G. Eberhart, Clonal selection drives genetic divergence of metastatic medulloblastoma, Nature 482(7386)
529-533.
[18] P.P. Gardner, D. Jennifer, J.G. Tate, E.P. Nawrocki, D.L. Kolbe, L. Stinus, A.C. Wilkinson, R.D. Finn, G.J. Sam, S.R. Eddy, Rfam: updates to the RNA families database, Nucleic Acids Research (suppl_1) (2008) suppl_1. [19] H. Li, R. Durbin, Fast and accurate short read alignment with Burrows-Wheeler transform, Bioinformatics
(2009).
[20] K. Ana, G.J. Sam, miRBase: annotating high confidence microRNAs using deep sequencing data, Nucleic Acids Research (D1) (2013) D1. [21] L.-C. Wan, F. Wang, X. Guo, S. Lu, Z. Qiu, Y. Zhao, H. Zhang, J. Lin, Identification and characterization of small non-coding RNAs from Chinese fir by high throughput sequencing, Bmc Plant Biol 12(1)
146.
[22] M.R. Friedlander, S.D. Mackowiak, N. Li, W. Chen, N. Rajewsky, miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades, Nucleic Acids Research 40(1) 37-52. [23] M.I. Love, W. Huber, S. Anders, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2, Genome Biol 15(12)
550.
[24] H.J. Wu, Y.-K. Ma, T. Chen, M. Wang, X.-J. Wang, PsRobot: a web-based plant small RNA meta-analysis toolbox, Nucleic Acids Research 40(W1)
W22-W28.
[25] M. Ashburner, C.A. Ball, J.A. Blake, D. Botstein, J.M. Cherry, Gene ontology: Tool for the unification of biology, Nat Genet 25(1) (2000) 25-29. [26] D.V. Klopfenstein, Z. Liangsheng, P.B. S., R. Fidel, W.V. Alex, N. Aurélien, M.C. J., Y.J. M., B. Olga, W. Mark, GOATOOLS: A Python library for Gene Ontology analyses, Sci Rep-Uk 8(1)
10872-.
[27] K. Minoru, G. Susumu, KEGG: Kyoto Encyclopedia of Genes and Genomes, Nucleic Acids Research (1) (2000) 1.
516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 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
[28] C. Xie, X. Mao, J. Huang, Y. Ding, J. Wu, S. Dong, L. Kong, G. Gao, C.Y. Li, L. Wei, KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases, Nucleic Acids Research (suppl_2) (2011) suppl_2. [29] J. Zhang, H. Zhang, A.K. Srivastava, Y. Pan, J. Bai, J. Fang, H. Shi, J.-K. Zhu, Knock-down of rice microRNA166 confers drought resistance by causing leaf rolling and altering stem xylem development, Plant Physiology
pp.01432.2017.
[30] W.X. Li, Y. Oono, J. Zhu, X.-J. He, J.-M. Wu, K. Iida, X.-Y. Lu, X. Cui, H. Jin, J.-K. Zhu, The Arabidopsis NFYA5 Transcription Factor Is Regulated Transcriptionally and Posttranscriptionally to Promote Drought Resistance, Plant Cell 20(8)
2238-2251.
[31] M. Zhou, H. Luo, Role of microRNA319 in creeping bentgrass salinity and drought stress response, Plant signaling & behavior 9(4)
e28700.
[32] L. Zhai, Z. Liu, X. Zou, Y. Jiang, F. Qiu, Y. Zheng, Z. Zhang, Genome-wide identification and analysis of microRNA responding to long-term waterlogging in crown roots of maize seedlings, 61(15) (2010) 4157-4168. [33] Sunkar, R., Novel and Stress-Regulated MicroRNAs and Other Small RNAs from Arabidopsis, Plant Cell 16(8)
2001-2019.
[34] H.H. Liu, X. Tian, Y.-J. Li, C.-A. Wu, C.-C. Zheng, Microarray-based analysis of stress-regulated microRNAs in Arabidopsis thaliana, Rna 14(5)
836-843.
[35] T.P. Frazier, G. Sun, C.E. Burklew, B. Zhang, Salt and Drought Stresses Induce the Aberrant Expression of microRNA Genes in Tobacco, 49(2)
159-165.
[36] Y.S. Dongmei Liu, Zhixiang Chen*, Diqiu Yu*, Ectopic expression of miR396 suppresses GRF target gene expression and alters leaf growth in Arabidopsis, Physiol Plantarum 136(2) (2009) 223-236. [37] W. Liu, Z. Yonggang, L. Xiaowei, W. Xingchao, D. Yuanyuan, W. Nan, L. Xiuming, C. Huan, Y. Na, C. Xiyan, Tissue-Specific Regulation of Gma-miR396 Family on Coordinating Development and Low Water Availability Responses, Front Plant Sci 8
1112-.
[38] Y. Ren, L. Chen, Y. Zhang, X. Kang, Z. Zhang, Y. Wang, Identification of novel and conservedPopulus tomentosamicroRNA as components of a response to water stress, Funct Integr Genomics 12(2)
327-339.
[39] V. Eldem, U.C. Akcay, E. Ozhuner, Y. Bakır, S. Uranbey, T. Unver, Genome-wide identification of miRNAs responsive to drought in peach (Prunus persica) by high-throughput deep sequencing, Plos One 7(12) (2012) e50298. [40] M. Arshad, B.A. Feyissa, L. Amyot, B. Aung, A. Hannoufa, MicroRNA156 improves drought stress tolerance in alfalfa ( Medicago sativa ) by silencing SPL13, Plant Science 258
122-136.
[41] M.K. Rezaei, Z.-S. Shobbar, M. Shahbazi, R. Abedini, S. Zare, Glutathione S-transferase (GST) family in barley: Identification of members, enzyme activity, and gene expression pattern, J Plant Physiol 170(14)
1277-1284.
[42] R. Sharma, A. Sahoo, R. Devendran, M. Jain, Over-Expression of a Rice Tau Class Glutathione S-Transferase Gene Improves Tolerance to Salinity and Oxidative Stresses in Arabidopsis, Plos One 9 (2014). [43] X. Jing, X. Xiao-Juan, T. Yong-Sheng, P. Ri-He, X. Yong, Z. Wei, Y. Quan-Hong, Z. Hong, Transgenic Arabidopsis Plants Expressing Tomato Glutathione S-Transferase Showed Enhanced Resistance to Salt and Drought Stress, Plos One 10(9)
e0136960-.
558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582
[44] G. Li, X. Peng, L. Wei, G. Kang, Salicylic acid increases the contents of glutathione and ascorbate
583
Tables list
and temporally regulates the related gene expression in salt-stressed wheat seedlings, Gene 529(2) 321-325. [45] D.C. Rees, E. Johnson, O. Lewinson, ABC transporters: the power to change, Nat Rev Mol Cell Biol 10(3)
218-227.
[46] L. Yang, Y. Jin, W. Huang, Q. Sun, F. Liu, X. Huang, Full-length transcriptome sequences of ephemeral plant Arabidopsis pumila provides insight into gene expression dynamics during continuous salt stress, BMC genomics 19(1) (2018) 717. [47] S.J. Suh, Y.-F. Wang, A. Frelet, N. Leonhardt, M. Klein, C. Forestier, B. Mueller-Roeber, M.H. Cho, E. Martinoia, J.I. Schroeder, The ATP Binding Cassette Transporter AtMRP5 Modulates Anion and Calcium Channel Activities in Arabidopsis Guard Cells, J Biol Chem 282(3)
1916-1924.
[48] M. Lee, Y. Choi, B. Burla, Y.-Y. Kim, B. Jeon, M. Maeshima, J.-Y. Yoo, E. Martinoia, Y. Lee, The ABC transporter AtABCB14 is a malate importer and modulates stomatal response to CO2, Nature Cell Biology 10(10)
1217-1223.
[49] M. Jain, J.P. Khurana, Transcript profiling reveals diverse roles of auxin-responsive genes during reproductive development and abiotic stress in rice , Febs J 276(11) (2009) 3148-3162. [50] R. Yue, S. Tie, T. Sun, L. Zhang, Y. Yang, J. Qi, S. Yan, X. Han, H. Wang, C. Shen, Genome-wide identification and expression profiling analysis of ZmPIN, ZmPILS, ZmLAX and ZmABCB auxin transporter gene families in maize (Zea mays L.) under various abiotic stresses, Plos One 10(3) (2015) e0118751. [51] C. Chai, S.P. K., Comprehensive Analysis and Expression Profiling of the OsLAX and OsABCB Auxin Transporter Gene Families in Rice (Oryza sativa) under Phytohormone Stimuli and Abiotic Stresses, Front Plant Sci 7.
Table 1. Summary of the RNA sequencing results
584 Raw reads Clean reads miRNA
Total_num
s_M_B_1
s_M_B_2
s_M_B_3
s_M_T_1
s_M_T_2
s_M_7_3
101375907 80957233( 79.86%) 18931949 23.39%
18045994 14627819 (81.06%) 4282093 29.27%
15684160 11814745 (75.32%) 3159378 26.74%
19063526 16355079 (85.79%) 3654379 22.34%
15020934 11799819 (78.56%) 2558940 21.69%
17026663 14451133 (84.87%) 3218371 22.27%
16534630 11908638 (72.02%) 2058788 17.29%
585 586 587
Table 2. Primer sequences
Genes 18S -F 18S -R cca-miR396a-3p-RT cca-miR396a-3p-F cca-miR396a-5p-RT cca-miR396a-5p-F
Primer sequence CTTCGGGATCGGAGTAATGA GCGGAGTCCTAGAAGCAACA CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGTTTCCCA ACACTCCAGCTGGGGTTCAAGAAAGCTG CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGCAGTTC ACACTCCAGCTGGGTTCCACAGCTTTCTT
cca-miR395b-RT cca-miR395b-F Nov-m0411-3p-RT Nov-m0411-3p-F Nov-m0067-5p-RT Nov-m0067-5p-F All R ABCB1-F ABCB1-R GST-F GST-R PHOT2-F PHOT2-R MPC1-F MPC1-R AT4G37250-F AT4G37250-R ADH5-F ADH5-R APS1-F APS1-R PAP1-F PAP1-R SPL6-F SPL6-R SPL15-F SPL15-R GRF3-F GRF3-R GRF6-F GRF6-R ATHB-15-F ATHB-15-R REVT-F REVT-R SCL6-F SCL6-R SCL15-F SCL15-R
CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGGAGTTCC ACACTCCAGCTGGGCTGAAGTGTTTGGA CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGAGAGCCT ACACTCCAGCTGGG CAGGCTGTGTGTCC CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGAATCTCT ACACTCCAGCTGGGTTTTTGTCGGGACC TGGTGTCGTGGAGTCG AACCCACTAAAAAGCAAACGATTC GAGGAGCACCACCGTCATTAT CCGGTTCACAAGAAGATTCCG GTTCGACATTAAGAATTAACACCCA GCAGCTTTATGCAGGCTCAC CCTGAAAGGTGTGCGACCA GTGACTGAAGGGCGACTGAA CAGTCGAGGGCAGTGATCTG TCCTTTTCCTGCCCAAGGAC CTTATCCCTTCCTCGCCGTC TTTGCCACACCGATGCCTAT ATCTCGGACTAGGAAGGCAGA GATCGTTTCAGGCTTTCGCC GCCTTTGAAGGCCCATTACAAC CTTGACGTCATTTGCGGTCC TAGCAAATTAAGCACGACAGC CAAAAGAAGTTGCAGAAGGCGA CGCTTACCCGGTAGGAAATAGAA CACCCATTGTGGTGCCATTG GCAAGATCAGTGATGCAGCG GGTGGAGAAGTTCAGATGAAGAGT GTGGTGGTGGTGGTAGTGAG TGGCTACAAGGATTCCATTCACA TGCAGCTAAGTACTTGAACACAA TGACGGCGAATAGAACTCGG CTGGAGGGGGTAGTAGGGTC AGGCGTTGCGGTATATCAGG TCCAACGGCATTCTCGTCAA ATCATCGCCTCCCCTCTCTT ATGCGGTCACCTTGAGGTTC ATGTGATACTGGCACGGCTC GCGTTGGCTTTACGACACTC
588 589
Figure titles list:
590 591
Fig.1. Expression profile of all detected miRNAs Fig.2. Cluster analysis of the differentially expressed known miRNAs
592 593 594 595 596
Fig.3. Cluster analysis of the differentially expressed novel RNAs Fig.4. Enrichment of target Genes in GO pathway Fig.5. Enrichment of target Genes in KEGG pathway Fig.6. Expression of miRNAs and the target genes