Skip to main content

Whole-transcriptome RNA sequencing reveals the global molecular responses and ceRNA regulatory network of mRNAs, lncRNAs, miRNAs and circRNAs in response to copper toxicity in Ziyang Xiangcheng (Citrus junos Sieb. Ex Tanaka)

Abstract

Background

Copper (Cu) toxicity has become a potential threat for citrus production, but little is known about related mechanisms. This study aims to uncover the global landscape of mRNAs, long non-coding RNAs (lncRNAs), circular RNAs (circRNAs) and microRNAs (miRNAs) in response to Cu toxicity so as to construct a regulatory network of competing endogenous RNAs (ceRNAs) and to provide valuable knowledge pertinent to Cu response in citrus.

Results

Tolerance of four commonly used rootstocks to Cu toxicity was evaluated, and ‘Ziyang Xiangcheng’ (Citrus junos) was found to be the most tolerant genotype. Then the roots and leaves sampled from ‘Ziyang Xiangcheng’ with or without Cu treatment were used for whole-transcriptome sequencing. In total, 5734 and 222 mRNAs, 164 and 5 lncRNAs, 45 and 17 circRNAs, and 147 and 130 miRNAs were identified to be differentially expressed (DE) in Cu-treated roots and leaves, respectively, in comparison with the control. Gene ontology enrichment analysis showed that most of the DEmRNAs and targets of DElncRNAs and DEmiRNAs were annotated to the categories of ‘oxidation-reduction’, ‘phosphorylation’, ‘membrane’, and ‘ion binding’. The ceRNA network was then constructed with the predicted pairs of DEmRNAs-DEmiRNAs and DElncRNAs-DEmiRNAs, which further revealed regulatory roles of these DERNAs in Cu toxicity.

Conclusions

A large number of mRNAs, lncRNAs, circRNAs, and miRNAs in ‘Ziyang Xiangcheng’ were altered in response to Cu toxicity, which may play crucial roles in mitigation of Cu toxicity through the ceRNA regulatory network in this Cu-tolerant rootstock.

Background

Copper (Cu) is an essential micronutrient for plant growth and development. As a redox-active transition element, Cu plays key roles in photosynthesis, respiration, C and N metabolism, oxidative stress protection, lignification, pollen fertility, and ethylene perception [1,2,3,4]. Most of the functions rendered by Cu are ascribed to enzymatically bound Cu, which catalyzes redox reactions [1]. In plants, there are more than 100 Cu-containing proteins, such as plastocyanin (PC), copper/zinc superoxide dismutase (CSD), cytochrome c oxidase (COX), laccase (LAC), diamine oxidases (DAO), and polyphenol oxidases [1,2,3]. Despite being essential, Cu is easily toxic to plants, even at a supra-optimal level [5, 6]. Excessive Cu inhibits plant growth and uptake of other mineral nutrients and alters enzyme systems, membrane integrity and many other biochemical and physiological processes [5]. Unfortunately, in the last decades, Cu contamination has become a global issue due to the long-term use of Cu-containing fungicides and bactericides, wastewater irrigation, and unconscionable Cu mining [6, 7]. In China, Cu is ranked as the fourth most contaminating heavy metal of agricultural lands [7, 8]. Thus, it is important and pressing to gain better understanding of physiological and molecular responses to Cu excess.

It is worth mentioning that plants have evolutionarily developed a tightly regulated system to balance the uptake, efflux, chelation, distribution, and utilization of Cu [9, 10]. In this system, a number of functional proteins (such as Cu transporters and chaperones), transcription factors (TFs), as well as non-coding RNAs (ncRNAs) may be involved in regulation of Cu homeostasis. For instance, Cu transporters, including CTR-like Cu transporters (COPTs), P-type heavy metal ATPases (HMAs), yellow stripe-like (YSL) proteins, zinc/iron-regulated transporter (ZRT/IRT)-related proteins (ZIPs), and cation diffusion facilitators (CDFs), were shown to directly participate in Cu uptake, transport, and distribution [2, 3, 9, 11]. Among them, COPTs (such as COPT1, COPT2, and COPT6) are mainly responsible for Cu uptake from soil and redistribution to reproductive organs [12, 13]. The HMAs such as AtHMA5, AtHMA7/RAN1, AtHMA6/PAA1, and AtHMA8/PAA2 of Arabidopsis thaliana are involved in Cu transport into the xylem, chloroplast, or thylakoid [7, 14,15,16]. YSL16 protein functions in recycling of Cu from older tissues to young tissues and grains [17]. Apart from these functional proteins, a conserved transcription factor called SPL7 (Squamosa Promoter binding-Like 7) has been shown to be a central regulator of Cu homeostasis [18, 19]. SPL7 regulates the expression of multiple targets, such as COPT1, COPT2, COPT6, cupric reductases FRO4 and FRO5, and Cu-regulated microRNAs (miR397, miR398, miR408 and miR857), that contain reiterative Cu-response elements (CuREs) with a GTAC motif in their promoters under Cu deficiency or excess [18, 20,21,22].

In addition to the protein-coding RNAs, emerging evidence has revealed that ncRNAs also play essential regulatory roles in plant responses to abiotic stress [23]. MicroRNAs (miRNAs), a major class of ncRNA with a length of 19 to 24 nucleotides, participate in Cu homeostasis by repressing translation or directly degrading Cu-related proteins in plants [22, 24]. Previous studies indicate that Cu deficiency induced the expression of miR397, miR408, and miR857, which repressed a number of Cu-containing proteins, including CSD1, CSD2, LAC, COX subunit 5b (COX5b-1), and Cu chaperone for SOD (CCS1) [20, 25, 26], while Cu excess downregulated miR398 to induce the expression of CSD1 and CSD2 [27]. Recently, two types of ncRNA, long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs), were discovered in plants. LncRNAs belong to a group of ncRNA longer than 200 nucleotides and can regulate the expression of genes through cis−/trans-acting or miRNA sponges [28, 29]. CircRNAs are endogenous covalently closed circular RNAs that are generated by alternative circularization [30]. A competing endogenous RNA (ceRNA) hypothesis has been proposed that the lncRNAs, circRNAs, mRNAs and pseudogenes can act as ceRNAs to competitively bind common miRNA response elements (MREs), and thus regulate a wide range of biological and developmental processes [31,32,33,34,35,36,37]. These ceRNAs are also named miRNA sponges to sequester specific miRNAs and inhibit their functions [30, 31, 34, 35]. In several studies, the ceRNA regulatory theory has been used to uncover molecular mechanisms of plant biology. For example, a ceRNA network was constructed to dissect function of lncRNAs in phosphate starvation of rice [32]. A complex ceRNA network consisting of lncRNAs, mRNAs, and miRNAs was established for maize seed development [33]. A large number of circRNAs, possibly acting as ceRNAs, were found in the ethylene pathway of tomato [38]. Recently, the ceRNA networks were also reported to be related to A. thaliana leaf development, tomato flowering, and cucumber heat stress response [35, 37, 39, 40]. However, whether lncRNAs and circRNAs participate in Cu homeostasis by acting as ceRNAs in plants remains to be investigated.

Citrus is widely grown worldwide. With extensive application of Cu-containing bactericides for controlling citrus canker disease, Cu toxicity has become a potential threat for citrus. However, relevant researches to understand Cu toxicity response are greatly limited. In this study, whole-transcriptome RNA sequencing (RNAseq) of Ziyang Xiangcheng (Citrus junos Sieb. ex Tanaka), a widely used citrus rootstock in China, was performed so as to uncover the global molecular responses to Cu toxicity at both protein-coding RNAs (mRNAs) and ncRNAs (lncRNAs, miRNAs, and circRNAs) levels. Moreover, the ceRNA networks were constructed for further revealing the underlying regulatory mechanisms in response to Cu toxicity.

Results

Evaluation of tolerance to Cu toxicity in the citrus rootstocks

To compare tolerance to Cu toxicity, four widely used citrus rootstocks, trifoliate orange (TO), ‘Ziyang Xiangcheng’ (XC), red tangerine (RT), and ‘Shatian’ pummelo (ST), were subjected to excessive Cu treatment, followed by evaluation of plant phenotypic and physiological parameters. As shown in Fig. 1, after 25 d of treatment, top leaves of TO and ST showed a yellow color, while those of XC and RT were almost normal (similar to the phenotype of CK). Although the relative increase rate of plant height was significantly suppressed upon Cu toxicity, XC had a minimal impact among the four rootstocks. Chlorophyll contents were also significantly reduced by Cu treatment, but XC had a higher level than the other three rootstocks at 25 and 40 d of treatment. In addition, excessive Cu treatment resulted in a drastic increase of MDA in TO, RT, and ST relative to CK, but increase of MDA in XC was not conspicuous. These results indicated that XC was the most tolerant genotype to Cu toxicity among the tested rootstocks. Therefore, XC was selected for high-throughput RNAseq to reveal global transcriptome of mRNA, lncRNA, circRNA and miRNA in response to Cu toxicity.

Fig. 1
figure1

Phenotype and physiological changes of four rootstocks under normal (CK) and Cu toxicity. Four common citrus rootstocks, trifoliate orange (TO), ‘Ziyang Xiangcheng’ (XC), red tangerine (RT), and ‘Shatian’ pummelo (ST), were grown under normal and excess of Cu (187.5 μM, 125×) conditions. After 25 d of treatment, representative pictures were photographed. Plant height increase rate and malonaldehyde (MDA) content were determined at 40 d, while chlorophyll a (chla) and b (chlb) contents were measured at 0 d, 25 d, and 40 d. Data are means ± SE (n = 3). Different letters indicate significant differences at P < 0.05 by Tukey test

Global responses of mRNA to Cu toxicity

From RNAseq data, we identified 30,123 protein-coding genes in leaves and roots of XC by using pummelo as reference genome, and their log2FC values are presented as Volcano Plot pictures in Fig. 2a and b, which ranged from − 8.2 to 7.9 in roots and from − 5.3 to 4.8 in leaves. Among all of these genes, 5734 (2162 up-regulated, 3572 down-regulated) and 222 (132 up-regulated, 90 down-regulated) DEmRNAs were identified in Cu-treated root (CuR/CKR) and leaf (CuL/CKL) groups, respectively (Fig. 2c and Additional file 1: Table S1). The number of DEmRNAs in the root was significantly higher than those in the leaf, suggesting that the root had more predominant responses to Cu toxicity. Moreover, 137 DEmRNAs were common between CuR/CKR and CuL/CKL, implying that they might participate in the basic response process under Cu toxicity. A heat map of DEmRNAs showed the general expression profiles of DEmRNAs in each treatment and also showed that the three repeats of each treatment always clustered together while the Cu-treated group and the CK group were clustered separately (Fig. 2d).

Fig. 2
figure2

Identification and analysis of differentially expressed mRNAs (DEmRNAs) under Cu toxicity. a, b Volcano Plot pictures showing log2FC values and FDR of mRNAs in CuL/CKL and CuR/CKR. c Venn diagram showing the number of DEmRNAs in CuL/CKL and CuR/CKR. d Heat map of all DEmRNAs. e GO enrichment of DEmRNAs in the root

To explore the functions of the DEmRNAs, GO annotation and GO enrichment analysis were performed. Base on GO annotation, we found that most DEmRNAs in the leaf and root were annotated to the metabolic process, single-organism process, localization process and response to stimulus under biological process (BP), to membrane, cell, organelle, and extracellular region under cellular component (CC), and to binding, catalytic activity, transporter activity, antioxidant activity, transcription factor activity and nutrient reservoir activity under molecular function (MF) (Additional file 10: Figure S1a and b). In particular, there were 317 DEmRNAs in response to stimulus, 596 in the membrane, 45 involved in antioxidant activity, 1 with metallochaperone activity, 32 with molecular transducer activity, 126 nucleic acid binding transcription factors, 18 with nutrient reservoir activity, 19 with receptor activity, and 188 with transporter activity in the root (Additional file 8: Table S8). GO enrichment analysis showed that the oxidation-reduction process, phosphorylation process, photosynthesis process, lignin catabolic process, phenylpropanoid catabolic process, membrane, dehydrogenase activity, peroxidase activity, protein kinase activity, iron ion binding, and heme binding were significantly enriched in the leaf and root (Fig. 2e and Additional file 10: Figure S1c), and these GO terms of DEmRNAs possible have primarily participated in mitigation of Cu toxicity.

Global responses of lncRNA to Cu toxicity

Apart from mRNAs, we also identified 243 known lncRNAs and 1033 novel lncRNAs in XC from the RNAseq data by blasting to known lncRNAs of citrus in the GREENC database and performing CNCI, CPC, CPAT and PfamScan analysis (Fig. 3a). Comparison of the genomic characterizations of the lncRNAs with mRNAs showed that their transcripts were similar in length distribution, except lncRNA had relative higher numbers of long transcripts (> 4500 bp) than mRNA; for exon number, a higher percentage of lncRNAs had 2 to 4 exons; in addition, lncRNAs had a shorter ORF length and lower FPKM value (Fig. 3b-e). At a cutoff with an absolute value of log2FC > 1 and FDR < 0.05, 164 (103 up-regulated, 61 down-regulated) and 5 (1 up-regulated, 4 down-regulated) DElncRNAs were identified in CuR/CKR and CuL/CKL groups, respectively (Fig. 3f and Additional file 2: Table S2). The log2FC values of DElncRNAs ranged from − 10.0 to 13.2 in the root, and − 11.1 to 8.8 in the leaf. The general expression profiles of DElncRNAs are shown in Fig. 3g. Similar to DEmRNAs, the DElncRNAs in the Cu-treated group and CK group were clustered separately, while their three repeats were clustered together.

Fig. 3
figure3

Identification and analysis of differentially expressed lncRNAs (DElncRNAs) under Cu toxicity. a Venn diagram showing the number of lncRNAs identified by CNCI, CPC, PfamScan and CPAT methods. be Comparison of lncRNA with mRNA with respect to the transcript length, exon number, ORF length and FPKM value. f Venn diagram showing the number of DElncRNAs in CuL/CKL and CuR/CKR. g Heat map of all DElncRNAs. h GO enrichment of targets of DElncRNAs in the root

To explore the potential functions of these DElncRNAs, their cis- and trans-targeted mRNAs were predicted with bioinformatics methods (Additional file 3: Table S3). GO annotation of the targets of DElncRNAs in the root showed that they were annotated in 16 terms under BP (mainly involved in metabolic process, cellular process, single-organism process, localization and response to stimulus), 13 terms under CC (mainly involved in membrane, cell, organelle and macromolecular complex), and 10 terms under MF (mainly involved in binding, catalytic activity, transporter activity, electron carrier activity, transcription factor activity, and antioxidant activity) (Additional file 11: Figure S2). GO enrichment analysis of these targets showed that significantly enriched GO terms were the photosynthesis process, phosphorylation process, oxidation-reduction process, lignin catabolic process, phenylpropanoid catabolic process, MCM complex, photosystem I reaction center, extracellular region, membrane, dehydrogenase activity, peroxidase activity, protein kinase activity, iron ion binding, and heme binding (Fig. 3h).

Global responses of circRNA to Cu toxicity

In total, 2404 circRNAs were identified in the leaf and root of XC, and 60.48, 28.62, and 10.90% of them belong to the intergenic region type, exon type, and intron type, respectively (Fig. 4c). The sequence length distribution of circRNAs is shown in Fig. 4a, and most of them were 10,000 bp to 50,000 bp, or shorter than 1200 bp. Chromosome 2 (chr2) included maximum circRNAs, followed by chr3, chr5, and chr8 (Fig. 4b). Log2FC values of circRNAs are displayed in Fig. 4d and e, and 45 (28 up-regulated, 17 down-regulated) and 17 (11 up-regulated, 6 down-regulated) DEcircRNAs were identified in CuR/CKR and CuL/CKL groups, respectively (Fig. 4f and Additional file 4: Table S4), among which, 1 DEcircRNAs were common between CuR/CKR and CuL/CKL. A heat map showed the general expression profiles of DEcircRNAs in each treatment, and most DEcircRNAs were highly expressed in the CuR group (Fig. 4g). These results suggest that the circRNAs are also involved in the responses to Cu toxicity.

Fig. 4
figure4

Identification and analysis of differentially expressed circRNAs (DEcircRNAs) under Cu toxicity. ac Sequence length, chromosome, and type distribution of all identified circRNAs. d, e Volcano Plot pictures showing log2FC values and FDR of circRNAs in CuL/CKL and CuR/CKR. f Venn diagram showing the number of DEcircRNAs in CuL/CKL and CuR/CKR. g Heat map of all DEcircRNAs

Global responses of miRNA to Cu toxicity

In the present study, a total of 23,333,512, 24,526,156, 21,822,295, and 25,871,923 raw reads were generated from CKL, CKR, CuL and CuR small RNA libraries, respectively. Of these raw reads, we obtained 15,050,388, 15,173,260, 13,474,928 and 13,748,074 clean reads after removing adaptor sequences, low-quality sequences, and reads shorter than 18 nt and longer than 32 nt. The lengths of most clean reads were 20–24 nt (Fig. 5a). Small RNA classification showed that 81% of clean reads were rRNA (42%) and unmatched (39%), and there were also 14% miRNA, 5% tRNA, and 1% of other types (Fig. 5b). From 14% of clean reads, we finally identified 149 known miRNAs and 336 novel miRNAs. The top 10 expressed miRNAs in each sample are shown in Fig. 5c and d, and miR166a-3p and nov-m0105-3p exhibited the highest expression abundance. From known miRNAs 12 (10 up-regulated, 2 down-regulated) and 3 (all up-regulated) DEmiRNAs, and from novel miRNAs 135 (26 up-regulated, 109 down-regulated) and 127 (42 up-regulated, 85 down-regulated) DEmiRNAs were identified in CuR/CKR and CuL/CKL groups, respectively (Fig. 5e and f and Additional file 5: Table S5). The general expression profiles of these DEmiRNAs are shown in Fig. 5g and h. Their expression levels exhibited obvious differences between CK and Cu-treated samples and between root and leaf samples.

Fig. 5
figure5

Identification and analysis of differentially expressed miRNAs (DEmiRNAs) under Cu toxicity. a Length distribution of all identified small RNAs. b Percentage of different types of small RNAs. c, d Top 10 expressed known and novel miRNAs in each sample. e, f Venn diagram showing the number of known (e) and novel (f) DEmiRNAs in CuL/CKL and CuR/CKR. g, h Heat map of all known (g) and novel (h) DEmiRNAs

Targeted mRNAs of these DEmiRNAs are listed in Additional file 6: Table S6. We found that 84.7% of DEmRNAs in the leaf (188/222) and 81.0% of DEmRNAs in the root (4642/5734) were targeted by one or multiple DEmiRNAs. GO annotation of the targets of DEmiRNAs in the root showed that most of them were annotated to the metabolic process, cellular process, single-organism process, localization process, and response to stimulus under BP, to membrane, cell, organelle and macromolecular complex under CC, and to binding, catalytic activity and transporter activity under MF (Additional file 12: Figure S3a and b). GO enrichment analysis showed that the significantly enriched GO terms of the targets of known DEmiRNAs in root were photosynthesis, microtubule-based movement, carbohydrate metabolic process, DNA polymerase complex, microtubule, membrane, cellulose synthase activity, microtubule binding, and catalytic activity, while those of novel DEmiRNAs were microtubule-based movement process, phosphorylation process, protein modification process, oxidation-reduction process, DNA polymerase complex, kinesin complex, extracellular region, membrane, protein kinase activity, transferase activity, anion binding, and catalytic activity (Fig. 6a and b).

Fig. 6
figure6

GO enrichment of targets of known (a) and novel (b) DEmiRNAs in the root

CeRNA regulatory network in response to Cu toxicity

To reveal the global regulatory network of protein-coding RNAs and ncRNAs under Cu toxicity, a ceRNA network was constructed using DEmiRNAs, DEmRNAs, DElncRNAs, and DEcircRNAs based on ceRNA theory. In total, 5739 DEmRNAs, 64 DElncRNAs, and 5 DEcircRNAs were predicted as targets of 251 miRNAs in the root and leaf. When their correlation was further filtered with SCC < − 0.5, we obtained 3819 DEmiRNA-DEmRNA and 10 DEmiRNA-DElncRNA interactions in the root and 12 DEmiRNA-DEmRNA interactions in the leaf (Fig. 7). In this ceRNA network, Nov-m0238-3p, Nov-m0074-5p, Nov-m0183-3p, miR166c-5p, Nov-m0128-3p, Nov-m0328-5p, miR165a-5p, and miR535c are involved in more than 100 nodes, suggesting that they may act as core regulators. In addition, lncRNAs including TCONS_00012501, TCONS_00012960, TCONS_00025983, TCONS_00027274, TCONS_00034874, TCONS_00036810, TCONS_00042747, TCONS_00051884 and TCONS_00062474 participated in the network.

Fig. 7
figure7

CeRNA network constructed with all DEmRNAs, DElncRNAs and DEmiRNAs in the root (a) and leaf (b). Color represents the up-regulated (red color) and down-regulated (blue color) levels

Considering that the network contains enormous information and each one cannot be displayed in the figure, we constructed a mini-ceRNA network by reducing the mRNA amount. We identified 284 known Cu-related mRNAs that were reported directly or indirectly in previous literatures from all DEmRNAs of the root; they included Cu-related regulators (SPL, YSL, CDPK, MAPK, and SUMO E3 Ligase, etc.), transporters (COPT, HMA, PPA, CDF, ZIP, OPT, and ABC transporter, etc.) and enzymes (LAC, CSD, CCS, PC, and COX, etc.) according to their functional description (Additional file 7: Table S7 and Additional file 8: Table S8). The mini-ceRNA network was constructed with these 284 DEmRNAs and all DEmiRNAs, DElncRNAs, and DEcircRNAs. Finally, only 261 DEmiRNA-DEmRNA and 10 DEmiRNA-DElncRNA interactions (SCC < − 0.5) containing 18 DEmiRNAs, 129 DEmRNAs, and 9 DElncRNAs were obtained and are displayed in Fig. 8. DEmiRNAs including miR166a-5p, miR395c, miR535c, miR395k, miR166c-5p, miR165a-5p, and miR399a were involved in more than 20 nodes, and all of them were up-regulated in the CuR. A known Cu-related key miRNA, miR398b, was identified in the network, which down-regulated and interacted with 5 DEmRNAs and 2 DElncRNAs. In addition, many known Cu-related key mRNAs such as SPL (Cg5g011720, Cg6g012520, Cg7g016770), YSL (Cg7g013630), HMA (Cg5g002920, Cg5g002930, Cg4g021370), ABC transporter (Cg5g018290, Cg5g027620, Cg3g011050, Cg3g009290, and Cg5g021160 etc.), LAC (Cg6g004840, Cg7g002640, Cg6g004880) and ZIP (Cg8g019240, Cg9g029160, Cg2g037280) were down-regulated in the network.

Fig. 8
figure8

Mini-ceRNA network constructed with 284 known Cu-related DEmRNAs and all DElncRNAs and DEmiRNAs in the root. Color represents the up-regulated (red color) and down-regulated (blue color) levels

qRT-PCR validated expression correlation between miRNAs and ceRNAs under Cu toxicity

To confirm the results of RNAseq and validate expression correlation of miRNA and their targets, six miRNAs (miR398b, miR8175, miR157c-3p, miR166a-5p, miR165a-5p, and Nov-m0284-5p) and their targeted mRNAs and lncRNAs were selected from the ceRNA network for qRT-PCR. As shown in Table 1, the qRT-PCR results agreed well with RNAseq data except several low-abundance mRNAs were undetectable. In addition, the miRNA and its targets showed quite correct up- or down-regulated relationships. For example, miR398 was down-regulated and all of its predicted targets were up-regulated; miR8175 was up-regulated and all of its targets were down-regulated. This result not only suggests reliability of RNAseq data in this study, but also validates the negatively correlated expression between miRNAs and ceRNAs.

Table 1 Confirmation of expression levels of miRNAs and their targeted lncRNAs and mRNAs based on qRT-PCR

Expression patterns of candidate RNAs between XC and TO

To further understand the Cu tolerance mechanism of XC, comparative analysis of the expression of known Cu-related miRNAs and ceRNAs was performed between Cu-tolerant XC and Cu-sensitive TO. As shown in Fig. 9, miR398b and its targets (Cg1g001620, Cg1g015400, Cg5g003300 and TCONS_0012960) were down-regulated or up-regulated in both roots of XC and TO at 1 d, 3 d and 5 d of Cu toxicity. However, the absolute values of log2FC in XC were significantly higher than those in TO at most time points. The similar changes were found for miR157c-3P and miR535c, as well as their targeted mRNAs and lncRNAs. These results suggest that the Cu-tolerant XC had occurred much more drastic expression of the Cu-related miRNAs, mRNAs and lncRNAs under Cu toxicity, which may be the important reason that leads to higher Cu tolerance in XC.

Fig. 9
figure9

qRT-PCR analysis of the candidate miRNAs, mRNAs and lncRNAs in ‘Ziyang Xiangcheng’ (XC) and trifoliate orange (TO). Three miRNAs and their predicted targets were selected to determine the expression in roots of XC and TO which were treated with excess and normal Cu for 1 d, 3 d and 5 d. FC represents the fold changes of the relative expression levels between CuR and CKR. The data are the means ± SE of three technical replicates. Asterisk (*) indicates the significant difference at P < 0.05 by LSD testing

Discussion

The ceRNA hypothesis is now widely accepted since it was reported several years ago [31]. Despite great progress in understanding human disease using the ceRNA theory [41], relatively fewer studies have been carried out in plants. In the ceRNA network, miRNAs play critical role in connection and regulation of different ceRNAs, such as mRNAs, lncRNAs and circRNAs. The mRNAs can be transcribed into proteins that exhibit direct functions, whereas the lncRNAs and circRNAs, not transcribed, can indirectly influence the expression of mRNAs by competitively binding the common miRNAs [31]. Based on this theory, we try to construct a ceRNA regulatory network so as to elucidate the mechanisms underlying tolerance to excessive Cu in the tolerant genotype.

In this study, 96.3% of DEmRNAs were predicted as the targets of 251 DEmiRNAs, which constitute 3819 DEmiRNA-DEmRNA interactions in the root and 12 in the leaf. This result suggests that most of the DEmRNAs may be regulated by the miRNAs under Cu toxicity. From the constructed ceRNA network, we found several miRNAs, including substantially down-regulated miR398b and up-regulated miR165a-5p, miR166a-5p, miR166c-5p, miR395c, miR395k, miR399a and miR535c, which have been previously reported to play a role in metal stress response. Among them, miR398 was shown to be significantly down-regulated under Cu excess but up-regulated under Cu deficiency [20, 22, 27, 42]. MiR395b and miR395c of A. thaliana were found to be up-regulated under Cu and cadmium (Cd) toxicity [43]. In addition, miR166 and miR399 were shown to play important roles in manganese (Mn), Cd, arsenic (As) or aluminum (Al) toxicity [44]. Although no direct evidence is provided to support the function of miR535 in metal stress, another member sharing same superfamily and high sequence similarity, miR156 has been well documented to function in Cd, Al, Mn, and As toxicity by targeting SPL genes [44, 45]. Interestingly, these metal stress-related miRNAs were predicted to target and regulate a large number of known metal stresses-related genes in the ceRNA network (Fig. 8, Additional file 7: Table S7 and Additional file 8: Table S8). For example, miR398b is predicted to target five DEmRNAs that significantly up-regulated by excessive Cu level, including putative ATP-binding cassette (ABC) transporters (Cg5g003300, Cg6g016060), oligopeptide transporter (OPT, Cg1g001620), multicopper oxidase (MCO, Cg1g015400), and mitogen-activated protein kinase kinase kinase (MAPKKK, Cg8g024350). It is well known that ABC transporters are one of the largest families of plant proteins that are suggested to be involved in metal ion uptake, transport and heavy metals detoxification [46, 47]. The OPT proteins consisting of YSL and PT clades have been demonstrated to participate in metal homeostasis through the translocation of metal-chelates [48]. The MCOs, which belong to the blue Cu proteins containing one to six Cu atoms per molecule, are suggested to function in redox reactions [49]. As an important signal transduction pathway, the MAPK signaling cascade plays key roles in response to different extracellular stimuli, and it has been previously found to be activated by excess of Cu and Cd [50, 51]. It is worth mentioning that except the miR398-targeted genes, other DEmRNAs, including putative SPLs, YSLs, HMAs, ABC transporters, LACs and ZIPs, were targeted by up-regulated miRNAs. These target mRNAs have also been shown to play essential roles in homeostasis of Cu or other metals. Therefore, we speculate that the better tolerance of Cu toxicity in XC is ascribed, at least in part, to miRNA-mediated regulation of the Cu-related genes, which allows it to maintain Cu homeostasis through alteration of signal transduction, Cu uptake, transport and detoxification, and antioxidant capacity. This speculation is corroborated by the differential expression levels of the miRNAs and targeted mRNAs in two genotypes with contrasting tolerance to Cu toxicity (Fig. 9).

Accumulating evidence indicates that the lncRNAs and circRNAs, two classes of new ncRNAs, act as ceRNAs (or named ‘miRNA sponges’, ‘target mimics’) to function in a variety of biological processes. In colorectal cancer, the lncRNA H19 was reported to derepress the endogenous genes, Vimentin, ZEB1, and ZEB2, by targeting miR138 and miR200a [52]. An autophagy promoting factor (APF) lncRNA was shown to interact with miR188-3p, and thus affected ATG7 expression and autophagic cell death in heart [53]. In plants, overexpression of lncRNA IPS1 in A. thaliana led to increased accumulation of PHO2, which is targeted by miR399 and functions in phosphate uptake [54]. In addition, circRNAs have been revealed to function as miRNA sponges. For example, a circHIPK3 could sponge to nine miRNAs with 18 potential binding sites and significantly affected human cell growth [55]. In this study, we identified nine DElncRNAs in the ceRNA network, suggesting that these citrus lncRNAs may also function as miRNA sponges to participate in tolerance to Cu toxicity. However, it should be pointed out that most of DElncRNAs were not included in the ceRNA network. It is thus assumed that these lncRNAs involve in Cu toxicity via other mechanisms, such as transcribing as miRNA precursors, induction of DNA methylation, modulation of chromatin modification, or serving as transcription enhancers [56]. It is surpring and unexpected that all of the DEcircRNAs were not included in the ceRNA network. One of the possible reasons is that the criteria used to predict the ceRNA pairs and interactions between DERNAs is too stringent.

Based on the ceRNA theory and the ceRNA network constructed in this study, we propose a mode of action for the differentially expressed mRNAs, miRNAs, lncRNAs and circRNAs in response to the Cu toxicity (Fig. 10). Under Cu toxicity, miRNAs act as the key regulators to modulate up-regulation or down-regulation of important Cu transporters, Cu proteins and Cu regulators (TFs and kinases). As a consequence, uptake, efflux and distribution of Cu will be activated or suppressed, and cell integrity was protected, leading to enhanced tolerance to Cu toxicity. In this process, some lncRNAs can act as ceRNA to competitively bind the MRE of miRNAs, which may indirectly affect the expression of mRNA. Moreover, circRNAs possibly act as another type of ceRNA to sequester Cu-responsive miRNAs and suppress their function.

Fig. 10
figure10

The proposed model in response to Cu toxicity in citrus. Under Cu toxicity, miRNAs act as the key regulators that directly target important Cu transporters, Cu proteins and Cu regulators (TFs and kinases). In this process, some lncRNAs and circRNAs can act as ceRNA to competitively bind the MRE of miRNAs, which may indirectly affect the expression of mRNAs

Conclusions

Tolerance evaluation showed that XC was most tolerant to Cu toxicity. Whole-transcriptome RNAseq helped us to identify 5734 (2162 up-regulated, 3572 down-regulated) and 222 (132 up-regulated, 90 down-regulated) DEmRNAs in Cu-treated root and leaf, respectively, in which 1243 were considered to be key candidates in response to excessive Cu. We also identified 243 known and 1033 novel lncRNAs, of which 164 (103 up-regulated, 61 down-regulated) and 5 (1 up-regulated, 4 down-regulated) were significantly differentially expressed in the Cu-treated root and leaf, respectively. From 2404 identified circRNAs, only 45 (28 up-regulated, 17 down-regulated) and 17 (11 up-regulated, 6 down-regulated) DEcircRNAs were identified in root and leaf, respectively, exposed to excessive Cu. In addition, 149 known miRNAs and 336 novel miRNAs were predicted, and 147 and 130 of them were responsive to Cu toxicity in the root and leaf, respectively. GO enrichment analysis showed that most of the DEmRNAs and targets of DElncRNAs and DEmiRNAs were implicated in oxidation-reduction, phosphorylation, membrane, and ion binding. A ceRNA network consisting of differentially expressed mRNAs, miRNAs, and lncRNAs was constructed, which further revealed the critical roles of these DERNAs in tolerance to Cu toxicity.

Methods

Plant materials and treatments

Seeds of four commonly used citrus rootstocks, trifoliate orange (Poncirus trifoliata L. Raf.), ‘Ziyang Xiangcheng’ (Citrus junos Sieb. ex Tanaka), red tangerine (C. reticulata Blanco), and ‘Shatian’ pummelo (C. grandis) were collected from Citrus Research Institute of Southwest University (Chongqing, China). To evaluate the tolerance to Cu toxicity, seeds of the rootstocks were germinated in plastic containers filled with quartz sand at a temperature of 28 °C and a relative humidity of 80%. Uniform seedlings were transplanted into fresh quartz sand washed with distilled water at 25 °C under a 16-h photoperiod (50 μmol·m− 2·s− 1) for sand culture. During sand culture, ½-strength Hoagland’s solution composed of 4 mM Ca (NO3)2, 1.5 mM KNO3, 0.5 mM NH4H2PO4, 1 mM MgSO4, 50 μM Fe-EDTA, 15 μM H3BO3, 10 μM MnSO4, 5 μM ZnSO4, 1.5 μM CuSO4, and 1 μM H2MoO4, was irrigated every 4 d. After 30 d of growth, half of the seedlings were irrigated with the above mentioned solution (used as control, CK), while the rest seedlings were irrigated with ½-strength Hoagland’s solution containing 187.5 μM CuSO4 (125×) for 40 d, which was considered as Cu treatment. Three biological replicates were performed for each treatment.

For RNAseq and quantitative real time PCR (qRT-PCR) validation, seeds of ‘Ziyang Xiangcheng’ and trifoliate orange were first sterilized with 2% sodium hypochlorite for 15 min. After removal of seed coats, the seeds were germinated at a temperature of 28 °C and a relative humidity of 80% for one week. Uniform seedlings were transferred to the above-mentioned normal ½-strength Hoagland’s solution for hydroculture at 25 °C under a 16-h photoperiod (50 μmol·m− 2·s− 1), and the solution was renewed every 10 d. After 30 d of growth, half of the seedlings were renewed with ½-strength Hoagland’s solution containing 75 μM CuSO4 (50×) for excess Cu treatment, while the others were renewed with normal solution (CK). After 1 d, 3 d and 5 d of treatment, the roots and leaves of the seedlings were sampled, frozen in liquid nitrogen, and stored at − 80 °C, respectively. Three biological replicates were performed for each treatment.

Measurement of plant height and contents of chlorophyll and malonaldehyde

Plant height (PH) of the aerial part was measured with a ruler at 0 (PH1) and 40 d of Cu treatment (PH2). Relative increase rate of the plant height was calculated as (PH2-PH1)/PH1 × 100%. Contents of chlorophyll and malonaldehyde (MDA) were determined as described by Fu et al. [57].

RNA extraction, library preparation, and RNA sequencing

Total RNA was extracted from the roots and leaves of CK (abbreviated as CKR and CKL) and Cu-treated samples (abbreviated as CuR and CuL) using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. RNA quality and integrity were estimated with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Only high-quality RNA samples (1.8 < OD260/280 < 2.2, OD260/230 ≥ 2.0, RIN ≥ 7.0, 28S/18S ≥ 1.0) were used to construct the sequencing library.

For mRNA, LncRNA, and circRNA sequencing, 5 μg of total RNA was used to prepare ribosomal RNA (rRNA) removed strand-specific library using a TruSeq Stranded Total RNA Library Prep with the Ribo-Zero Plant Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. There were three biological replicates per treatment, and a total of 12 libraries were prepared. For small RNA sequencing, 4 small RNA libraries were constructed with 3 μg of total RNA from CKL, CKR, CuL and CuR samples and the Truseq Small RNA sample prep Kit (Illumina, San Diego, CA, USA). After libraries were quantified by a TBS-380 Fluorometer (Turner Biosystem, Sunnyvale, CA, USA), deep RNA sequencing was performed using an Illumina HiSeq X Ten platform at Shanghai Majorbio Bio-Pharm Biotechnology Co. Ltd. (Shanghai, China).

Read mapping and transcriptome assembly

The paired-end raw reads were trimmed and quality controlled by SeqPrep (version 1.1, https://github.com/jstjohn/SeqPrep) and Sickle (version 1.33, https://github.com/najoshi/sickle) with default parameters. Then, clean reads were separately aligned to the Pummelo (Citrus grandis) genome [39] in orientation mode using Tophat2 software [58] (version 2.0.13, http://tophat.cbcb.umd.edu/). The mapped reads of each sample were assembled by Cufflink (version 2.2.1, http://cufflinks.cbcb.umd.edu/) in a reference-based approach.

Differentially expressed mRNA and gene ontology (GO) enrichment analysis

To identify differentially expressed mRNAs (DEmRNAs) between CK and Cu-treated samples, the expression level of each transcript was calculated according to the fragments per kilobase of exon per million mapped reads (FRKM) method. RSEM (version 1.2.31, http://deweylab.biostat.wisc.edu/rsem/) was used to quantify gene abundance [59]. R statistical package software EdgeR (version 3.14.0) [60] was utilized for differential expression analysis with an absolute value of log2FC > 1 and FDR < 0.05. GO annotation and functional enrichment analysis of DEmRNAs were carried out using the Omicshare oneline platform (http://www.omicshare.com/tools/).

Identification of lncRNAs and prediction of their target genes

For identification of novel lncRNAs, the transcripts that overlapped with known protein-coding genes on the same strand, transcripts with a fragment count < 3, transcripts shorter than 200 nt, and an open reading frame (ORF) longer than 300 nt were first discarded [61, 62]. Then, we used the Coding Potential Calculator (CPC, version 0.9), Coding-Non-Coding index (CNCI, version 2.0), Coding Potential Accessment Tool (CPAT, version 1.2.4), and Pfam Scan (version 1.6) to filter transcripts with coding potential. The overlapped outputs from CPC, CNCI, CPAT, and Pfam Scan were considered reliably expressed novel lncRNAs. The transcripts were also used to blast the citrus lncRNA sequences that were collected in the GREENC database (http://greenc.sciencedesigners.com/wiki/Main_Page), and the hits with e_value <1E-5 and matched bases ratio > 90% were identified as know lncRNAs. All identified lncRNAs were classified into intergenic, intronic, and antisense lncRNAs using the cuffcompare program in the Cufflinks suite. The expression level of each lncRNA was calculated according to the FRKM method. Differentially expressed lncRNAs (DELncRNAs) were extracted with an absolute value of log2FC > 1 and FDR < 0.05 by EdgeR. The potential cis- and trans-target mRNAs of DELncRNAs were predicted according to the position on the chromosome. The cis-targets were searched within a 10-kb window upstream or downstream of the lncRNA [63]. The trans-targets were predicted as described by Ou et al. [64].

Identification and analysis of circRNAs

CircRNA Identifier (CIRI, version 1.2) and CIRCexplorer (version 1.1.7) tools were used to identify circRNAs in this study. CIRI scans Sequence Alignment/Map (SAM) files and detects junction reads with paired chiastic clipping (PCC) signals and paired-end mapping (PEM) and GT-AG splicing signals as described by Gao et al. [65]. CIRCexplorer obtains junction reads via a two-step TopHat and TopHat-Fusion mapping strategy as described by Zhang et al. [66]. The overlapped outputs from CIRI and CIRCexplorer were used for further analysis. The expression level of each circRNA was calculated according to the Spliced Reads per Billion Mapping (SRPBM) method. Differentially expressed circRNAs (DEcircRNAs) were extracted with an absolute value of log2FC > 1 and FDR < 0.05 by DEGseq (version 1.30.0).

Identification and analysis of miRNAs

The raw data were first quality controlled with Fastx-toolkit software (version. 0.0.13, http://hannonlab.cshl.edu/fastx_toolkit/) to obtain clean small RNA reads by filtering low-quality bases (sanger base quality < 20), sequencing adapters, reads shorter than 18 nt, and reads longer than 32 nt. The assembled unique sequences with clean reads were then BLAST searched against the Rfam database (version 12.1, http://rfam.sanger.ac.uk/) to remove non-miRNA sequences (rRNA, tRNA, snoRNA, etc.). The remaining reads were used to predict known miRNAs through a BLAST search of the miRbase, version 21.0 (http://www.mirbase.org/), and novel miRNAs through analysis of the hairpin structure of the miRNA precursor with Mireap (version 0.2, http://sourceforge.net/projects/mireap) software. The expression level of each miRNA was calculated according to the transcripts per million reads (TPM) method. Differentially expressed miRNAs (DEmiRNAs) were extracted with an absolute value of log2FC > 1 and FDR < 0.005 by DEGseq. Target prediction of DEmiRNAs was performed with psRobot local version 1.01 [67].

Construction and analysis of ceRNAs regulatory network

To reveal the roles and interactions of lncRNAs, circRNAs, miRNAs, and mRNAs, we constructed a lncRNA-circRNA-miRNA-mRNA regulatory network based on the ceRNA hypothesis. psRobot [67] was used to predict the pairs of miRNA-lncRNA, miRNA-mRNA, and miRNA-circRNA. The pairwise correlations of miRNA-lncRNA, miRNA-mRNA and miRNA-circRNA were evaluated using the Spearman correlation coefficient (SCC) and the matched pairs’ expression profile data as described by He et al. [40]. The interaction network was built with RNA pairs of SCC < − 0.5 and visually displayed using Cytoscape software (version 2.8) [68].

qRT-PCR validation of DEmRNAs, DElncRNAs, and DEmiRNAs

qRT-PCR was performed to validate the expression levels of DEmRNAs, DElncRNAs, and DEmiRNAs. Total RNA and small RNA were extracted from samples of CuR and CKR using an RNAprep pure Plant Kit (TIANGEN, Cat#DP432, China) and miRcute miRNA isolation kit (TIANGEN, Cat#DP501, China), respectively. First-strand cDNA was synthesized from 1 μg of total RNA with the HiScript® II Q RT SuperMix (Vazyme, Cat#R223, China) for qPCR of mRNA and with the lnRcute lncRNA First-Strand cDNA Synthesis Kit (TIANGEN, Cat#KR202, China) for qPCR of lncRNA. In addition, 1 μg of small RNA was used for cDNA synthesis using a miRNA 1st Strand cDNA Synthesis Kit (Vazyme, Cat#MR101, China) with the stem-loop primer designed by the stem-loop sequence (GTCGTATCCAGGGTCCGAGGTATTCGCACTGGATACGAC) except for the internal reference U6. qPCR was performed on the Bio-Rad CFX Connect RealTime system using ChamQ™ Universal SYBR® qPCR Master Mix (Vazyme, Cat#Q711), lnRcute lncRNA qPCR Detection Kit (TIANGEN, Cat#FP402) and miRNA Universal SYBR® qPCR Master Mix (Vazyme, Cat#MQ101) following the manufacturer’s instructions. The 2-ΔΔCT method was used to normalize and determine the RNA level relative to an internal reference gene, actin (Cs1g05000.1) or U6. All primers are included in Additional file 9: Table S9.

Availability of data and materials

All supporting data can be found within the manuscript and its additional supporting files.

Abbreviations

ABC transporter:

ATP-binding cassette transporter

BP:

Biological process

CC:

Cellular component

ceRNAs:

competing endogenous RNAs

circRNAs:

circular RNAs

CIRI:

circRNA identifier

CKL:

Cu-untreated leaf

CKR:

Cu-untreated root

CNCI:

Coding-non-coding index

COPTs:

CTR-like Cu transporters

COX:

Cytochrome c oxidase

CPAT:

Coding potential accessment tool

CPC:

Coding potential calculator

CSD:

Copper/zinc superoxide dismutases

Cu:

Copper

CuL:

Cu-treated leaf

CuR:

Cu-treated root

DAO:

Diamine oxidases

DE:

Differentially expressed

FC:

Fold change

FRKM:

The fragments per kilobase of exon per million mapped reads

GO:

Gene ontology

HMAs:

P-type heavy metal ATPases

LAC:

Laccase

lncRNAs:

Long non-coding RNAs

MAPKKK:

Mitogen-activated protein kinase kinase kinase

MCO:

Multicopper oxidase

MDA:

Malonaldehyde

MF:

Molecular function

miRNAs:

microRNAs

MREs:

miRNA response elements

OPT:

Oligopeptide transporter

ORF:

Open reading frame

PC:

Plastocyanin

PCC:

Paired chiastic clipping

PEM:

Paired-end mapping

PH:

Plant height

qRT-PCR:

Quantitative real time PCR

RNAseq:

RNA sequencing

SAM:

Sequence alignment/map

SCC:

Spearman correlation coefficient

SPL7:

Squamosa promoter binding-like 7

SRPBM:

The spliced reads per billion mapping

TFs:

Transcription factors

TO:

Trifoliate orange

TPM:

Transcripts per million reads

XC:

‘Ziyang Xiangcheng’

YSL:

Yellow stripe-like protein

References

  1. 1.

    Broadley M, Brown P, Cakmak I, Rengel Z, Zhao F. Function of nutrients: Micronutrients. p. 191–248. In: P. Marschner (eds.), Mineral Nutrition of Higher Plants, Elsevier, 2012.

  2. 2.

    Burkhead JL, Reynolds KAG, Abdel-Ghany SE, Cohu CM, Pilon M. Copper homeostasis. New Phytol. 2009;182:799–816.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  3. 3.

    Yruela Guerrero I. Copper in plants: acquisition, transport and interactions. Funct Plant Biol. 2009;36:409–30.

    Article  Google Scholar 

  4. 4.

    Yan J, Chia JC, Sheng H, Jung HI, Zavodna TO, Zhang L, Huang R, Jiao C, Craft EJ, Fei Z, et al. Arabidopsis pollen fertility requires the transcription factors CITF1 and SPL7 that regulate copper delivery to anthers and jasmonic acid synthesis. Plant Cell. 2017;29:3012–29.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  5. 5.

    Cambrolle J, Garcia JL, Figueroa ME, Cantos M. Evaluating wild grapevine tolerance to copper toxicity. Chemosphere. 2015;120:171–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  6. 6.

    Leng X, Jia H, Sun X, Shangguan L, Mu Q, Wang B, Fang J. Comparative transcriptome analysis of grapevine in response to copper stress. Sci Rep. 2015;5:17749.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  7. 7.

    Huang XY, Deng F, Yamaji N, Pinson SR, Fujii-Kashino M, Danku J, Douglas A, Guerinot ML, Salt DE, Ma JF. A heavy metal P-type ATPase OsHMA4 prevents copper accumulation in rice grain. Nat Commun. 2016;7:12138.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  8. 8.

    Zhao FJ, Ma Y, Zhu YG, Tang Z, McGrath SP. Soil contamination in China: current status and mitigation strategies. Environ Sci Technol. 2015;49:750–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  9. 9.

    Clemens S. Molecular mechanisms of plant metal tolerance and homeostasis. Planta. 2001;212:475–86.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  10. 10.

    Clemens S, Palmgren MG, Kramer U. A long way ahead: understanding and engineering plant metal accumulation. Trends Plant Sci. 2002;7:309–15.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  11. 11.

    Aguirre G, Pilon M. Copper delivery to chloroplast proteins and its regulation. Front Plant Sci. 2016;6:1250.

    PubMed  PubMed Central  Article  Google Scholar 

  12. 12.

    Sancenon V, Puig S, Mateu-Andres I, Dorcey E, Thiele DJ, Penarrubia L. The Arabidopsis copper transporter COPT1 functions in root elongation and pollen development. J Biol Chem. 2004;279:15348–55.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  13. 13.

    Jung HI, Gayomba SR, Rutzke MA, Craft E, Kochian LV, Vatamaniuk OK. COPT6 is a plasma membrane transporter that functions in copper homeostasis in Arabidopsis and is a novel target of SQUAMOSA promoter-binding protein-like 7. J Biol Chem. 2012;287:33252–67.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  14. 14.

    Andres-Colas N, Sancenon V, Rodriguez-Navarro S, Mayo S, Thiele DJ, Ecker JR, Puig S, Penarrubia L. The Arabidopsis heavy metal P-type ATPase HMA5 interacts with metallochaperones and functions in copper detoxification of roots. Plant J. 2006;45:225–36.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  15. 15.

    Shikanai T, Muller-Moule P, Munekage Y, Niyogi KK, Pilon M. PAA1, a P-type ATPase of Arabidopsis, functions in copper transport in chloroplasts. Plant Cell. 2003;15:1333–46.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  16. 16.

    Abdel-Ghany SE, Muller-Moule P, Niyogi KK, Pilon M, Shikanai T. Two P-type ATPases are required for copper delivery in Arabidopsis thaliana chloroplasts. Plant Cell. 2005;17:1233–51.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  17. 17.

    Zheng L, Yamaji N, Yokosho K, Ma JF. YSL16 is a phloem-localized transporter of the copper-nicotianamine complex that is responsible for copper distribution in rice. Plant Cell. 2012;24:3767–82.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  18. 18.

    Yamasaki H, Hayashi M, Fukazawa M, Kobayashi Y, Shikanai T. SQUAMOSA promoter binding protein-like 7 is a central regulator for copper homeostasis in Arabidopsis. Plant Cell. 2009;21:347–61.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  19. 19.

    Araki R, Mermod M, Yamasaki H, Kamiya T, Fujiwara T, Shikanai T. SPL7 locally regulates copper-homeostasis-related genes in Arabidopsis. J Plant Physiol. 2018;224:137–43.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  20. 20.

    Yamasaki H, Abdel-Ghany SE, Cohu CM, Kobayashi Y, Shikanai T, Pilon M. Regulation of copper homeostasis by micro-RNA in Arabidopsis. J Biol Chem. 2007;282:16369–78.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  21. 21.

    Bernal M, Casero D, Singh V, Wilson GT, Grande A, Yang H, Dodani SC, Pellegrini M, Huijser P, Connolly EL, et al. Transcriptome sequencing identifies SPL7-regulated copper acquisition genes FRO4/FRO5 and the copper dependence of iron homeostasis in Arabidopsis. Plant Cell. 2012;24:738–61.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  22. 22.

    Pilon M. The copper microRNAs. New Phytol. 2017;213:1030–5.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  23. 23.

    Wang J, Meng X, Dobrovolskaya OB, Orlov YL, Chen M. Non-coding RNAs and their roles in stress response in plants. Genom Proteom Bioinf. 2017;15:301–12.

    Article  Google Scholar 

  24. 24.

    Chien PS, Chiang CB, Wang Z, Chiou TJ. MicroRNA-mediated signaling and regulation of nutrient transport and utilization. Curr Opin Plant Biol. 2017;39:73–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  25. 25.

    Abdel-Ghany SE, Pilon M. MicroRNA-mediated systemic down-regulation of copper protein expression in response to low copper availability in Arabidopsis. J Biol Chem. 2008;283:15932–45.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  26. 26.

    Lu S, Li Q, Wei H, Chang M, Tunlaya-Anukit S, Kim H, Liu J, Song J, Sun YH, Yuan L, et al. Ptr-miR397a is a negative regulator of laccase genes affecting lignin content in Populus trichocarpa. P Natl Acad Sci USA. 2013;110:10848–53.

    CAS  Article  Google Scholar 

  27. 27.

    Sunkar R, Kapoor A, Zhu JK. Posttranscriptional induction of two cu/Zn superoxide dismutase genes in Arabidopsis is mediated by downregulation of miR398 and important for oxidative stress tolerance. Plant Cell. 2006;18:2051–65.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  28. 28.

    Chekanova JA. Long non-coding RNAs and their functions in plants. Curr Opin Plant Biol. 2015;27:207–16.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  29. 29.

    Ma K, Shi W, Xu M, Liu J, Zhang F. Genome-wide identification and characterization of long non-coding RNA in wheat roots in response to Ca2+ channel blocker. Front Plant Sci. 2018;9:244.

    PubMed  PubMed Central  Article  Google Scholar 

  30. 30.

    Li QF, Zhang YC, Chen YQ, Yu Y. Circular RNAs roll into the regulatory network of plants. Biochem Bioph Res. 2017;488:382–6.

    CAS  Article  Google Scholar 

  31. 31.

    Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta stone of a hidden RNA language? Cell. 2011;146:353–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  32. 32.

    Xu XW, Zhou XH, Wang RR, Peng WL, An Y, Chen LL. Functional analysis of long intergenic non-coding RNAs in phosphate-starved rice using competing endogenous RNA network. Sci Rep. 2016;6:20715.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  33. 33.

    Zhu M, Zhang M, Xing L, Li W, Jiang H, Wang L, Xu M. Transcriptomic analysis of long non-coding RNAs and coding genes uncovers a complex regulatory network that is involved in maize seed development. Genes. 2017;8:274.

    PubMed Central  Article  CAS  Google Scholar 

  34. 34.

    Ren GJ, Fan XC, Liu TL, Wang SS, Zhao GH. Genome-wide analysis of differentially expressed profiles of mRNAs, lncRNAs and circRNAs during Cryptosporidium baileyi infection. BMC Genomics. 2018;19:356.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  35. 35.

    Meng X, Zhang P, Chen Q, Wang J, Chen M. Identification and characterization of ncRNA-associated ceRNA networks in Arabidopsis leaf development. BMC Genomics. 2018;19:607.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  36. 36.

    Yuan Y, Li J, Xiang W, Liu Y, Shu J, Gou M, Qing M. Analyzing the interactions of mRNAs, miRNAs, lncRNAs and circRNAs to predict competing endogenous RNA networks in glioblastoma. J Neuro-Oncol. 2018;137:493–502.

    Article  CAS  Google Scholar 

  37. 37.

    Yang Z, Yang C, Wang Z, Yang Z, Chen D, Wu Y. LncRNA expression profile and ceRNA analysis in tomato during flowering. PLoS One. 2019;14:e0210650.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  38. 38.

    Wang Y, Wang Q, Gao L, Zhu B, Luo Y, Deng Z, Zuo J. Integrative analysis of circRNAs acting as ceRNAs involved in ethylene pathway in tomato. Physiol Plantarum. 2017;161:311–21.

    CAS  Article  Google Scholar 

  39. 39.

    Wang X, Xu Y, Zhang S, Cao L, Huang Y, Cheng J, Wu G, Tian S, Chen C, Liu Y, et al. Genomic analyses of primitive, wild and cultivated citrus provide insights into asexual reproduction. Nat Genet. 2017;49:765–72.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  40. 40.

    He X, Guo S, Wang Y, Wang L, Shu S, Sun J, et al. Physiol Plantarum. 2019. https://doi.org/10.1111/ppl.12997.

  41. 41.

    Su XQ, Xing JD, Wang ZZ, Chen L, Cui M, Jiang BH. microRNAs and ceRNAs: RNA networks in pathogenesis of cancer. Chinese J Cancer Res. 2013;25:235–9.

    Google Scholar 

  42. 42.

    Leng X, Wang P, Zhao P, Wang M, Cui L, Shangguan L, Wang C. Conservation of microRNA-mediated regulatory networks in response to copper stress in grapevine. Plant Growth Regul. 2017;82:293–304.

    CAS  Article  Google Scholar 

  43. 43.

    Gielen H, Remans T, Vangronsveld J, Cuypers A. Toxicity responses of cu and cd: the involvement of miRNAs and the transcription factor SPL7. BMC Plant Biol. 2016;16:145.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  44. 44.

    Gupta OP, Sharma P, Gupta RK, Sharma I. MicroRNA mediated regulation of metal toxicity in plants: present status and future perspectives. Plant Mol Biol. 2014;84:1–18.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  45. 45.

    Ding Y, Chen Z, Zhu C. Microarray-based analysis of cadmium-responsive microRNAs in rice (Oryza sativa). J Exp Bot. 2011;62:3563–73.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  46. 46.

    Wojas S, Hennig J, Plaza S, Geisler M, Siemianowski O, Sklodowska A, Ruszczynska A, Bulska E, Antosiewicz DM. Ectopic expression of Arabidopsis ABC transporter MRP7 modifies cadmium root-to-shoot transport and accumulation. Environ Pollut. 2009;157:2781–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  47. 47.

    Kim DY, Bovet L, Maeshima M, Martinoia E, Lee Y. The ABC transporter AtPDR8 is a cadmium extrusion pump conferring heavy metal resistance. Plant J. 2007;50:207–18.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  48. 48.

    Lubkowitz M. The oligopeptide transporters: a small gene family with a diverse group of substrates and functions? Mol Plant. 2011;4:407–15.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  49. 49.

    Hoegger PJ, Kilaru S, James TY, Thacker JR, Kues U. Phylogenetic comparison and classification of laccase and related multicopper oxidase protein sequences. FEBS J. 2006;273:2308–26.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  50. 50.

    Jonak C, Nakagami H, Hirt H. Heavy metal stress. Activation of distinct mitogen-activated protein kinase pathways by copper and cadmium. Plant Physiol. 2004;136:3276–83.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  51. 51.

    Yeh CM, Hsiao LJ, Huang HJ. Cadmium activates a mitogen-activated protein kinase gene and MBP kinases in rice. Plant Cell Physiol. 2004;45:1306–12.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  52. 52.

    Liang WC, Fu WM, Wong CW, Wang Y, Wang WM, Hu GX, Zhang L, Xiao LJ, Wan DCC, Zhang JF, et al. The lncRNA H19 promotes epithelial to mesenchymal transition by functioning as miRNA sponges in colorectal cancer. Oncotarget. 2015;6:22513–25.

    PubMed  PubMed Central  Google Scholar 

  53. 53.

    Wang K, Liu CY, Zhou LY, Wang JX, Wang M, Zhao B, Zhao WK, Xu SJ, Fan LH, Zhang XJ, et al. APF lncRNA regulates autophagy and myocardial infarction by targeting miR-188-3p. Nat Commun. 2015;6:6779.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  54. 54.

    Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubio-Somoza I, Leyva A, Weigel D, Garcia JA, Paz-Ares J. Target mimicry provides a new mechanism for regulation of microRNA activity. Nat Genet. 2007;39:1033–7.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  55. 55.

    Zheng QP, Bao CY, Guo WJ, Li SY, Chen J, Chen B, Luo YT, Lyu DB, Li Y, Shi GH, et al. Circular RNA profiling reveals an abundant circHIPK3 that regulates cell growth by sponging multiple miRNAs. Nat Commun. 2016;7:11215.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  56. 56.

    Hou J, Lu D, Mason AS, Li B, Xiao M, An S, Fu D. Non-coding RNAs and transposable elements in plant genomes: emergence, regulatory mechanisms and roles in plant development and stress responses. Planta. 2019;50:23–40.

    Article  CAS  Google Scholar 

  57. 57.

    Fu XZ, Khan EU, Hu SS, Fan QJ, Liu JH. Overexpression of the betaine aldehyde dehydrogenase gene from Atriplex hortensis enhances salt tolerance in the transgenic trifoliate orange (Poncirus trifoliata L. Raf.). Environ Exp Bot. 2011;74:106–13.

    CAS  Article  Google Scholar 

  58. 58.

    Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:R36.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  59. 59.

    Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  60. 60.

    Robinson MD, McCarthy DJ, Smyth GK. EdgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  61. 61.

    Ding ZH, Wu CL, Tie WW, Yan Y, He GY, Hu W. Strand-specific RNA-seq based identification and functional prediction of lncRNAs in response to melatonin and simulated drought stresses in cassava. Plant Physiol Bioch. 2019;140:96–104.

    CAS  Article  Google Scholar 

  62. 62.

    Quinn JJ, Chang HY. Unique features of long non-coding RNA biogenesis and function. Nat Rev Genet. 2016;17:47–62.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  63. 63.

    Jia H, Osak M, Bogu GK, Stanton LW, Johnson R, Lipovich L. Genome-wide computational identification and manual annotation of human long noncoding RNA genes. RNA. 2010;16:1478–87.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  64. 64.

    Ou L, Liu Z, Zhang Z, Wei G, Zhang Y, Kang L, Yang B, Yang S, Lv J, Liu Y, et al. Noncoding and coding transcriptome analysis reveals the regulation roles of long noncoding RNAs in fruit development of hot pepper (Capsicum annuum L.). Plant Growth Regul. 2017;83:1–16.

    Article  CAS  Google Scholar 

  65. 65.

    Gao Y, Wang J, Zhao F. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015;16:4.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  66. 66.

    Zhang XO, Wang HB, Zhang Y, Lu X, Chen LL, Yang L. Complementary sequence-mediated exon circularization. Cell. 2014;159:134–47.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  67. 67.

    Wu HJ, Ma YK, Chen T, Wang M, Wang XJ. PsRobot: a web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res. 2012;40:W22–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  68. 68.

    Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgements

We are grateful to prof. Jihong Liu from Huazhong Agricultural University for his revision of the manuscript.

Funding

This work was financially supported by the National Key Research and Development Program of China (2018YFD1000300, 2017YFD0202000), the National Natural Science Foundation of China (31772280), and the National Citrus Engineering Research Center (NCERC2019001). The funders provided the financial support to the RNAseq, but didn’t play role in in the design of the study, collection, analysis, interpretation of data, and in writing the manuscript.

Author information

Affiliations

Authors

Contributions

XZF and LZP conceived and designed the study. XZF, XYZ, and YZH analyzed the data. JYQ and XZ performed qRT-PCR. MY and LC prepared experimental materials. CPC and LLL measured physiological data. XZF wrote the paper. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Xing-Zheng Fu or Liang-Zhi Peng.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

List of 5734 and 222 DEmRNAs in the root and leaf, respectively.

Additional file 2: Table S2.

List of 164 and 5 DElncRNAs in the root and leaf, respectively.

Additional file 3: Table S3.

Cis and trans targets of DElncRNAs.

Additional file 4: Table S4.

List of 45 and 17 DEcircRNAs in the root and leaf, respectively.

Additional file 5: Table S5.

List of 147 and 130 DEmiRNAs in the root and leaf, respectively.

Additional file 6: Table S6.

DEmiRNAs targeted DEmRNAs in the root and leaf, respectively.

Additional file 7: Table S7.

Searched known Cu-related DEmRNAs by using the keywords reported in previous studies to query their functional descriptions.

Additional file 8: Table S8.

List of 1273 key DEmRNAs in response to Cu toxicity in citrus.

Additional file 9: Table S9.

List of primers used in qRT-PCR.

Additional file 10: Figure S1.

GO annotation of DEmRNAs in the leaf (a) and root (b), and GO enrichment of DEmRNAs in the leaf (c).

Additional file 11: Figure S2.

GO annotation of targets of DElncRNAs in the root.

Additional file 12: Figure S3.

GO annotation of targets of known (a) and novel (b) DEmiRNAs in the root.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Fu, X., Zhang, X., Qiu, J. et al. Whole-transcriptome RNA sequencing reveals the global molecular responses and ceRNA regulatory network of mRNAs, lncRNAs, miRNAs and circRNAs in response to copper toxicity in Ziyang Xiangcheng (Citrus junos Sieb. Ex Tanaka). BMC Plant Biol 19, 509 (2019). https://doi.org/10.1186/s12870-019-2087-1

Download citation

Keywords

  • Citrus
  • Copper
  • CeRNA
  • Transcriptome
  • Non-coding RNA