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)

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.


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 supraoptimal 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 longterm 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.
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, wholetranscriptome 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. 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

Plant materials and treatments
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 CuSO 4 (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 [41].

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-

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 [43]. R statistical package software EdgeR (version 3.14.0) [44] was utilized for differential expression analysis with an absolute value of log 2 FC > 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 [45,46]. 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 log 2 FC > 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 [47]. The trans-targets were predicted as described by Ou et al [48].

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 [49]. CIRCexplorer obtains junction reads via a two-step TopHat and TopHat-Fusion mapping strategy as described by Zhang et al [50]. 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 log 2 FC > 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 lowquality 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 log 2 FC > 1 and FDR < 0.005 by DEGseq. Target prediction of DEmiRNAs was performed with psRobot local version 1.01 [51].

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 [51] 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) [52]. 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.

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.

Global responses of mRNA to Cu toxicity
From RNAseq data, we identified 30123 protein-coding genes in leaves and roots of XC by using pummelo as reference genome, and their log 2 FC values are presented as Volcano Plot pictures in Fig.   2a and 2b, which ranged from -8.2 to 7.9 in roots and from -5.  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).
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: Fig. 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: Fig. 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 log 2 FC > 1 and FDR < 0.05, 164 (103 upregulated, 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 log 2 FC 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 Cutreated group and CK group were clustered separately, while their three repeats were clustered together.
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: Fig. 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 10000 bp to 50000 bp, or shorter than 1200 bp. Chromosome 2 (chr2) included maximum circRNAs, followed by chr3, chr5, and chr8 (Fig. 4b). Log 2 FC values of circRNAs are displayed in Fig. 4d and 4e, 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.

Global responses of miRNA to Cu toxicity
In the present study, a total of 23333512, 24526156, 21822295, and 25871923 raw reads were generated from CKL, CKR, CuL and CuR small RNA libraries, respectively. Of these raw reads, we obtained 15050388, 15173260, 13474928 and 13748074 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 5d, and miR166a-3p and nov-m0105-3p exhibited the highest expression abundance. From known miRNAs 12 (10 upregulated, 2 down-regulated) and 3 (all up-regulated) DEmiRNAs, and from novel miRNAs 135 (26 upregulated, 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 5f and Additional file 5: Table S5).
The general expression profiles of these DEmiRNAs are shown in Fig. 5g and 5h. Their expression levels exhibited obvious differences between CK and Cu-treated samples and between root and leaf samples.
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: Fig. 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 6b).

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 upregulated; 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.

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 log 2 FC 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.

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 [53], 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 [22,27,54,55].
MiR395b and miR395c of A. thaliana were found to be up-regulated under Cu and cadmium (Cd) toxicity [56]. In addition, miR166 and miR399 were shown to play important roles in manganese (Mn), Cd, arsenic (As) or aluminum (Al) toxicity [57]. 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 [57,58]. 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 [59,60]. The OPT proteins consisting of YSL and PT clades have been demonstrated to participate in metal homeostasis through the translocation of metal-chelates [61]. The MCOs, which belong to the blue Cu proteins containing one to six Cu atoms per molecule, are suggested to function in redox reactions [62]. 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 [63,64]. It is worth mentioning that except the miR398targeted 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 [65]. An autophagy promoting factor (APF) lncRNA was shown to interact with miR188-3p, and thus affected ATG7 expression and autophagic cell death in heart [66]. 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 [67]. 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 [68]. 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 [69]. 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.

Conclusions
Tolerance evaluation showed that XC was most tolerant to Cu toxicity. Whole-transcriptome RNAseq

Consent for publication
Not applicable.

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

Competing interests
The authors declare that they have no competing interests

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. Additional file 9: Table S9 List of primers used in qRT-PCR.     Heat map of all known (g) and novel (h) DEmiRNAs.  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. Figure 9 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.

Figure 10
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 mRNA.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.
Supplementary materials.rar