Cultivar- and tissue-specific transcriptomic changes induced by cadmium in Brassica parachinensis

Background One of the main pathways for cadmium (Cd) transfer from the environment to the human body is through the consumption of leafy vegetables, and Brassica leafy crops tend to be Cd hyper-accumulators. But its response strategies to Cd still lack of systematic study. Results To investigate Brassica response strategies to Cd, we identified two cultivars with different Cd translocation efficiencies and performed mutli-transcriptomic sequencing studies under Cd treatments. Certain transporter families exhibited different temporal expression profiles in the two cultivars and may underlie the different Cd translocation efficiencies. Cd induced a drastic reduction of a 22 nt small RNA, the footprint of a pentatricopeptide repeat (PPR) protein on the chloroplast ndhB transcript and the concomitant down-regulation of the ndhB transcript. A global reduction in the expression of PPR genes was found, revealing previously unknown effects of Cd on organellar gene expression. Analyses of microRNAs (miRNAs) and their target genes by small RNA and degradome sequencing not only revealed Cd-induced changes in miRNAs but also implicated the existence of a regulatory cascade involving bra-miR156, its target gene, and bra-miR397 and bra-miR398 in Cd stress responses. Conclusions The present findings help uncover the impact of Cd stress on the transcriptome of B. parachinensis and provide candidate genes and miRNAs for further investigation. on the post-transcriptional that on the of some chloroplast and nuclear PPR


Introduction
Cadmium (Cd) is a widespread heavy metal that is easily released into the environment by human industrial activities (Sanità Di Toppi and Gabbrielli 1999). It can be transferred to and chronically accumulate in human bodies via the food chain, and one of the main pathways of Cd transfer from plants to humans is through leafy vegetable consumption (Huang et al. 2017). Brassica crop species include various leafy vegetables, including cabbages and pak choi. The hypertolerance and hyper-accumulation of Cd in Brassica species (Lin and Aarts 2012;Rizwan et al. 2018) increases the risk of human exposure to this heavy metal.
Brassica parachinensis L.H.Bailey (Chinese flowering cabbage) is a Brassica leafy vegetable with a short vegetative cycle, large biomass and extensive planting adaptability (Qiu et al. 2011b). It has been recognized as a hyper-accumulator of Cd with different absorption ability depending on the genotype (Qiu et al. 2011a). Its

Cd absorption and accumulation mechanisms warrant investigation.
Cd is a non-essential metal in plants. It is typically absorbed in roots and binds to organic acids and amino acids in the xylem vessels, followed by its transfer to the shoots via xylem (Nakamura et al. 2008;Wu et al. 2015;Khan et al. 2017). Its accumulation in plant cells can damage mitochondria, the photosynthetic apparatus, membranes and other cell parts (Lin and Aarts 2012). Studies over the last two decades have shown that plants avoid or endure Cd stress by chelating Cd (Wójcik and Tukiendorf 2011;Gielen et al. 2017), compartmentalizing and sequestering Cd (Zhang et al. 2018a), promoting Cd efflux (Verret et al. 2004;Zhang et al. 2016) and/or regulating Cd influx (Yao et al. 2018). Several families of metal ion transporters, such as heavy metal ATPases (HMAs) (Zhang et al. 2016), natural resistance-associated macrophage proteins (NRAMPs) (Sasaki et al. 2012;Tang et al. 2017) and ATP-binding cassette subfamily C proteins (ABCCs) (Brunetti et al. 2015), have been implicated in the transport and sequestration of Cd. However, in 5 some cases, the functions of individual members within the same family differ between non-hyper-accumulators and hyper-accumulators. For example, OsHMA3 sequesters Cd into vacuoles in rice roots to reduce its Cd levels in shoots, while SpHMA3 up-regulation in Sedum plumbizincicola enhances Cd tolerance and accumulation in shoots (Sasaki et al. 2014;Liu et al. 2017;Zhang et al. 2018a). In general, individual gene functions depend on the genetic background in a complex way, and transcriptome analyses have shown that Cd induces a wide range of changes in different biological pathways. These include transport (Feng et al. 2018), reactive oxygen species scavenging (Yu et al. 2017a;Xu et al. 2018), secondary metabolite biosynthesis (Chen et al. 2017), cell wall modification (Gao et al. 2013) and DNA or RNA damage (Xu et al. 2018). To understand how plant cope with Cd stress and the related mechanisms, it is essential to clarify the up-and downstream activities and factors related to genes in these biological pathways.
MiRNAs are small non-coding RNAs that post-transcriptionally regulate gene expression by transcript cleavage, translation inhibition or secondary siRNA biogenesis from the target gene transcripts (Yu et al. 2017b). They have been implicated in abiotic stress tolerance (Shriram et al. 2016), including the response to heavy metal stress in plants (Noman and Aqeel, 2017). Some Cd responsive miRNAs have been identified in rice (Wang et al. 2009;Ding et al. 2011), B. napus (Zhou et al. 2012a;Jian et al. 2018), B. parachinensis (Zhou et al. 2017), soybean (Fang et al. 2013) and the hyper-accumulator Sedum alfredii (Han et al. 2016).
Additionally, specific functions of miR166 and miR395 in Cd accumulation and tolerance have been reported in rice and B. napus, respectively (Zhang et al. 2013;Ding et al. 2018). However, knowledge of miRNAs in Cd stress responses in B. parachinensis is rudimentary: only miRNA profiles under prolonged Cd treatments 6 were examined and global detection of miRNA activities by degradome sequencing was not performed (Zhou et al. 2017).
To investigate the effect of Cd stress on B. parachinensis, two cultivars with different Cd translocation efficiency were treated with a relatively high concentration of Cd for a short period, and the mRNA and sRNA transcriptomes as well as degradomes in the roots and leaves were analyzed. Our results showed that roots of the two cultivars largely exhibited similar transcriptional responses to Cd, with many genes related to metal ion transport being up-regulated by Cd. On the other hand, changes in gene expression and sRNA abundance were affected differently in the leaves of the two cultivars. The differences in Cd response in leaves were consistent with differences in their root-to-shoot Cd translocation capacity. The sRNA sequencing revealed a drastic reduction of an RNA footprint of a chloroplast pentatricopeptide repeat (PPR) protein and the concomitant reduction in the levels of the transcript from which the footprint originated. This, and the induction of many PPR genes by Cd, suggested that Cd induces changes in chloroplast RNA metabolism through nucleus-encoded and chloroplast-localized PPR proteins. Our sRNA-seq not only uncovered many novel miRNAs, but also identified many conserved and none-conserved miRNAs being regulated by Cd stress. In addition, degradome sequencing revealed the effects of these miRNAs in Cd response. Three miRNA-target modules seem to have regulatory relationships and form a gene expression cascade in Cd stress response. These genes and miRNAs can help to provide new potential directions for the breeding of food-safe Brassica crops.

Plant materials and Cd treatment
Six B. parachinensis cultivars were used in this experiment. Cultivars YQ, YLCT and JY40 were bought from the Vegetable Research Institute of Guangdong Academy of Agricultural Sciences. YQ80, TC and PT45 were obtained from Guangzhou Academy of Agricultural Sciences, Lianzhou Fengyu Agricultural Technology Co., Ltd., and Guangzhou Yangxin Seedling Co., Ltd., respectively. Seeds were sown in vermiculite and watered with 1/4 Hoagland nutrient solution. After 7 days (d), the seedlings were transferred to a simple hydroponic culture device with 1/2 Hoagland nutrient solution. After 8 d in hydroponic culture, Cd (Cd 2+ ) stress was applied by adding Cd(NO 3 ) 2 to the nutrient solution with final Cd 2+ concentrations of 50 μM or 5 μM.
The nutrient solution was changed every 3 d.

Biomass and Cd content determination
Fifteen-day-old seedlings were used for Cd treatment. Roots and shoots were collected separately at the designated time points to determine the biomass and Cd content. Tissue samples (0.5-1.0 g) were oven-dried then digested with 10 mL 65% HNO 3 using the Microwave Digestion System (EthosONE, Milestone, Italy). Cd contents were analyzed by inductively coupled plasma optical emission spectrometry (ICP-OES, Optima 7000DV, Perkin Elmer, USA). Three biological replicates were detected.

mRNA sequencing and data analysis
Fifteen-day-old TC and YQ seedlings were treated with 50 μM Cd. Roots and the first pair of leaves were separately harvested at 0, 24 and 48 HAC then immediately frozen in liquid nitrogen for total RNA extraction. Three biological replicates were performed for each time point, with each biological replicate consisting of three 8 seedlings grown in different hydroponic growth containers.
Total RNA was extracted using TRI Reagent (Molecular Research Center, Inc.). RNAseq libraries were prepared and sequenced by Beijing Novogene Co., Ltd. using the Illumina HiSeq 2500 platform to obtain 150 bp paired-end reads. Data analysis was performed using the pRNASeqTools pipeline (https://github.com/grubbybio/pRNASeqTools) in its mrna mode. First, the raw reads were trimmed using cutadapt (Martin 2011) with default settings, then the trimmed reads were mapped to the Brassica rapa v2.0 genome ) using STAR (Dobin et al. 2013). Second, the read counts of both introns and exons were combined as the total read count of one transcript using pRNASeqTools. Transcript levels were measured and normalized in TPM (transcripts per million) using TBtools (Chen et al. 2018) (https://github.com/CJ-Chen/TBtools/releases). Finally, DEGs were identified using DEseq2 with a fold change of 1.5 and P-value < 0.01 as the cutoff parameters (Love et al. 2014). Coding sequences (CDSs) of genes from the B. rapa v2.0 genome was compared against the Swiss-Prot protein sequence database (Bairoch and Apweiler 2000) using blastx with e-value < 1e-5 as the cutoff parameters (Camacho et al. 2009) to derive the functional descriptions of genes.
Gene ontology (GO) annotations were performed using GOanna tools available at the AgBase v2.0 website (http://agbase.arizona.edu/index.html) (McCarthy et al. 2011). GO enrichment analysis of DEGs was performed using TBtools (Chen et al. 2018) with Fisher's exact test. REVIGO was used to remove the duplicated GO terms (Supek et al. 2011). Trend analyses of gene expression at 0, 24 and 48 HAC were accomplished using the Short Time-series Expression Miner software (STEM) (Ernst and Bar-Joseph 2006) at the OmicShare tools website (http://www.omicshare.com/tools). The cluster dendrogram and correlation 9 coefficient matrix of all samples were respectively plotted using the 'hclust' and 'heatmap' R functions, and volcano plots of the DEGs were generated using 'ggplot2' in R (Wickham 2016). The above clustering analyses revealed that one biological replicate of each of the samples, YQ roots at 24 HAC, YQ roots at 48 HAC, YQ leaves at 0 HAC, TC leaves at 24 HAC and TC leaves at 48 HAC, was an outlier, only two biological replicates of these samples were retained for DEG analyses.

sRNA-seq library construction and data analysis
To construct sRNA-seq libraries, 20 µg of total RNA was resolved on a 15% urea-PAGE gel, and the ~10-40 nt region was excised from the gel. sRNAs were recovered by shaking the smashed gel in 0.4 M NaCl-DEPC overnight, followed by precipitation in a mixture containing glycogen, sodium acetate and ethanol. sRNA libraries were constructed following the instructions for the NEB Next Multiplex Small RNA Library PrepSet for Illumina kit (NEB E7300). The libraries were sequenced on an Illumina HiSeq 2500 instrument using a 50 bp single-read strategy by Beijing Berry Genomics Co., Ltd. sRNA data analysis was performed using the pRNASeqTools pipeline in its srna mode (Jia et al. 2017). The sequencing data were processed to remove the adapter sequence (AGATCGGAAGAGC) using cutadapt (Martin 2011). The trimmed reads were mapped to the reference genome, which included the B. rapa v2.0 chromosomes, mitochondrial genome sequence (GenBank accession: JF920285) and chloroplast genome sequence (GenBank accession: DQ231548, Wu et al. 2012), using ShortStack with default parameters (Axtell 2013). Normalization was performed by calculating the RPMR (reads per million of 45S rRNA reads) value (Li et al. 2016). The reference genome was tiled into 100 bp windows as individual bins, and reads whose 5' end nucleotides fell within a given bin were assigned to this bin (Li et al. 2016). Annotations for genes and TEs were downloaded from the Brassica Database (http://brassicadb.org/brad/index.php), and rRNAs were located by aligning known full-length rRNAs of plants onto the B. rapa v2.0 genome. As the B. rapa v2.0 genome does not contain MIR gene annotations, we annotated MIR genes by the alignment of MIR/miRNA annotations from the B. rapa v1.5 genome with the v2.0 genome . Prediction of novel miRNAs was performed using the MIRNA analysis mode of ShortStack (Axtell 2013), and only the novel MIR loci identified in at least two libraries were retained. DEMs were identified using DEseq2 with a fold change of 1.5 and P-value < 0.01 as the cutoff parameters (Love et al. 2014).

Degradome/PARE library construction and data analysis
Total RNA was used for degradome/PARE library construction. Equal amounts of total RNA from the samples at 24 HAC and 48 HAC were mixed as the Cd-treated samples. A total of 24 degradome/PARE libraries including three biological replicates were constructed using 75 µg of total RNA (Zhai et al. 2014). The libraries were sequenced using a 50 bp single-end read strategy on an Illumina HiSeq 2500 instrument at Beijing Berry Genomics Co., Ltd. The adapter sequences (TGGAATTCTCGGG) were removed, and reads shorter than 19 nt were filtered out using cutadapt (Martin 2011). Trimmed fastq files were converted to fasta files, and the degradome sites and the corresponding trigger miRNAs were identified using CleaveLand4 (Addo-Quaye et al. 2009). Degradome sites were selected with degradome category = 0 and degradome P-value < 0.05 as the cutoff parameters.
MiRNA cleavage sites identified in at least two biological replicates were retained and were then manually screened using their T-plots. The miRNA cleavage sites with high background noise in the T-plots were removed.

sRNA gel blotting and quantitative RT-PCR
Hydroponic TC and YQ seedlings were prepared as described above, and the samples were harvested at 0, 24 and 48 HAC with three biological replicates. Total RNA was extracted using TRI Reagent, and 20 µg total RNA was resolved in a 15% urea-PAGE gel. sRNA gel blotting assay was performed using a method with enhanced sensitivity (Pall and Hamilton 2008). Biotin-labeled antisense DNA oligonucleotides were used to detect sRNAs. The DNA oligonucleotide probes are listed in Table S1.
Total RNA was reverse-transcribed into single-stranded cDNA using the Brassica rapa subsp. chinensis was used as an internal reference gene , and the primers used for qRT-PCR are listed in Table S1.

Data availability
The sequencing data have been deposited into the NCBI Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/sra) under the BioProject ID PRJNA513391.

Different responses to Cd in two B. parachinensis cultivars
Fifteen-day-old seedlings of six widely grown B. parachinensis cultivars (YQ, YLCT, YQ80, TC, JY40 and PT45) were treated with 50 μM Cd for 6 d then analyzed for Cd accumulation and biomass. No significant difference was observed in the Cd content in roots and shoots among the six cultivars (Fig. S1a). The shoot biomass of all six cultivars was significantly reduced compared to the respective 0 μM Cd controls ( Fig. S1b). In contrast, changes in root biomass varied in the different cultivars: it was significantly reduced in TC and YQ80 and increased in YQ compared to the controls (Fig. S1b). An increase in root biomass in YQ but a decrease in TC after Cd stress was obvious upon visual inspection ( Fig. S1b and c). Based on the different responses to Cd, the YQ and TC cultivars were selected for further analysis.
The Cd content in YQ and TC seedlings during the initial 48 hours after Cd stress (HAC) was measured to determine if there were any differences in Cd absorption and translocation between the two cultivars. In both cultivars, the shoot Cd content continuously increased from 6 to 48 HAC, while the root Cd content was steady from 6 to 24 HAC then rapidly increased at 48 HAC (Fig. 1a). One difference between the cultivars was the significantly higher shoot Cd content in YQ than TC at 12 HAC (Fig.   1a). The Cd content results implied that root-to-shoot Cd translocation mainly occurred within 24 HAC in both cultivars and that this translocation was more rapid in YQ than in TC at 12 HAC (Fig. 1a). To validate the observed differences, we performed another treatment with a lower Cd concentration (5 μM) for a longer time (3 d and 6 d). At both 3 d and 6 d, there were significant differences in the shoot Cd content between YQ and TC (Fig. 1b), indicating more efficient Cd translocation in YQ than TC. Taken together, the findings show that B. parachinensis YQ and TC respond differently to Cd stress in terms of root growth and Cd translocation efficiency.

Cd-induced mRNA transcriptome changes in TC and YQ leaves and roots
To further investigate the response strategies of B. parachinensis to Cd stress, leaves and roots of 15-d-old YQ and TC seedlings treated with 50 μM Cd were 13 harvested at 0, 24 and 48 HAC for mRNA sequencing (mRNA-seq) analysis. As B.
parachinensis is considered a variety of Brassica rapa (Cheng et al. 2016a), the sequencing reads from the mRNA libraries were mapped to the B. rapa genome (version 2.0) , resulting in 89.17% to 93.98% total mapped reads (Table S2). Analysis of differentially expressed genes (DEGs) showed that the DEG numbers and the trend of up-and down-regulation were similar in TC and YQ roots but differed noticeably in the leaves of these two cultivars ( Fig. S2a and b). For example, there were more DEGs in YQ leaves than in TC leaves (Fig. S2b). In contrast to the considerable overlap of root DEGs between the two cultivars and at different Cd treatment time points, in leaves, YQ-specific DEGs were more than the overlapping DEGs between YQ and TC or the TC-specific DEGs (Fig. S2c). Cluster dendrogram and correlation coefficient analyses of the mRNA-seq data showed that YQ leaves at 24 HAC and 48 HAC constituted one cluster, and TC leaves at 0 HAC, 24 HAC and 48 HAC constituted another cluster ( Fig. S3a and b). On the other hand, YQ and TC roots at 0 HAC were in one cluster, and YQ and TC roots at 24 HAC and 48 HAC were in another cluster ( Fig. S3c and d). These data indicate that Cd stress had a more evident impact on gene expression in YQ leaves than TC leaves while roots of the two cultivars responded to Cd similarly.
The DEGs were further grouped according to the patterns of their expression changes over time (i.e., temporal expression profiles). The results showed that the significant profiles (p < 0.05) were similar between TC and YQ roots but different between TC and YQ leaves (Fig. 2a). For TC leaves, the most significant profile was Profile 3, in which expression decreased then increased, followed by Profile 4, in which expression increased then decreased. The enriched GO terms for genes in Profile 3 were mostly related to basal metabolism, such as "cellular amide metabolic 14 process", "cellular nitrogen compound metabolic process" and "gene expression" ( Fig. 2b and Table S3), suggesting that basal metabolism in TC leaves was initially negatively influenced then recovered. The enriched GO terms for Profile 4 genes were related to signaling and stimulus response, such as "gene expression", "signaling" and "regulation of cellular process" (Fig. 2b and Table S3), suggesting that these biological processes were initially enhanced in TC leaves then returned to normal levels. In YQ leaves, basal metabolism related genes were clustered in Profile 2, in which expression decreased then stabilized, and genes related to signaling and stimulus response were mostly found in Profile 5, in which expression increased then stabilized ( Fig. 2b and Table S3). The different expression trends of these similar groups of genes in TC and YQ leaves indicate that the impact of Cd stress on leaves persisted in YQ for a longer time but recovered more rapidly in TC.
In roots, the related genes in GO terms "metal ion transport", "transmembrane transport", "drug transmembrane transport", "divalent metal ion transport" and "transition metal ion transport" were significantly enriched among Profile 5 (increased then stabilized expression) and Profile 6 (increased expression) in both TC and YQ ( Fig. S4 and Table S4). The finding that Cd affected genes related to transmembrane transport suggests that these genes are involved in Cd response.

DEGs in response to Cd were related to metal ion transport and transmembrane transport
Cd absorption and translocation require transporter genes. We found that many of the DEGs with the GO term "metal ion transport" were members of the zinc/iron transporter (ZIP) family, the copper transporter (COPT) family and the HMA family (Table S5). For DEGs with the GO term "transmembrane transport", some of them were found to belong to the ATP-binding cassette transporter B subfamily (ABCB 15 family), ATP-binding cassette transporter C subfamily (ABCC family) and the yellow stripe-like (YSL) family (Table S5). We focused on these DEGs to examine their temporal expression profiles in response to Cd in shoots and roots. Most genes in the ABCC, ABCB, vacuolar cation/proton exchanger (CAX), COPT, HMA, ZIP and YSL families and all of the genes in the ABCG family were more highly expressed in roots than in leaves (Fig. 3). Although most of these genes had similar expression profiles in YQ and TC, several of them had different expression profiles in the two cultivars. For instance, six genes in the ABCC family (ABCC3-3, ABCC3-1, ABCC1-2, ABCC15, ABCC4-1 and ABCC7) were up-regulated in YQ leaves at 24 HAC but not in TC leaves, and HMA2-1 and HMA3 were more highly up-regulated in TC roots than in YQ roots. Two genes in the YSL family (YSL1-3 and YSL1-1) were markedly upregulated in YQ leaves at 24 HAC, and YSL1-2 transcript levels were high in TC leaves at 0 HAC. In all, although most of the transport genes were up-regulated in both cultivar roots, the specific DEGs with different expression profiles between TC and YQ in response to Cd may be related to the differences in Cd translocation efficiency between the two cultivars.

Cd-induced reduction of a chloroplast RNA fragment
We performed small RNA sequencing on the same samples to determine the responses of small RNAs to Cd stress (Table S6). Reads that mapped to the genomes (nuclear, mitochondrial and chloroplast) showed a minor 21 nt peak and a major 24 nt peak in all samples, consistent with observations in Arabidopsis (Li et al. 2016).
Surprisingly, the leaf samples had a peak at 22 nt, and the 22 nt peak was reduced in YQ leaves at 24 HAC and 48 HAC and reduced in TC leaves at 48 HAC (Fig. S5).
Mapping the 22 nt sRNAs to genomic features revealed that nearly half of them were derived from unannotated regions (designated as the "others" category), and the percentage of sRNAs in the "others" category was obviously reduced under Cd stress (Fig. 4a). Further analysis revealed that the "others" category was largely composed of a single sRNA (AGTTACTAATTCATGATCTGGC) located in the 5' UTR of two identical chloroplast genes (ndhB.1 / ndhB.2), so we henceforth refer to it as ndhB sRNA (Fig. 4b). The abundance of this sRNA was reduced by Cd stress in the leaves of both TC and YQ, but the decrease was more rapid in YQ (Fig. 4c). RNA gel blot analysis gave consistent results with the sRNA-seq data (Fig. 4d). The ndhB sRNA was also found in Arabidopsis thaliana at high abundance, and it was proposed to be a 'footprint' of the pentatricopeptide repeat (PPR) protein CRR2 (Ruwe and Schmitz-Linneweber 2012; Ruwe et al. 2016). We did not detect differential expression of CRR2 under Cd stress in RNA-seq (Fig. S6). Among the surrounding genes of this 22 nt sRNA, rps12, rps7 and ndhB belong to one precursor RNA (Hildebrand et al. 1988). RT-qPCR analysis showed that ndhB was down-regulated in YQ leaves at 24 and 48 HAC as well as in TC leaves at 48 HAC, while the levels of rps7 or rps12 transcripts were minimally affected when comparing to ndhB (Fig. 4e).
Therefore, Cd stress reduced the abundance of the ndhB sRNA and, in YQ leaves with higher levels of Cd accumulation (Fig. 1b), also down-regulated the expression of the NADH oxidoreductase genes ndhB.1 / ndhB.2 in leaves.
PPR proteins are encoded in the nucleus but regulate gene expression in mitochondria and chloroplasts in post-transcriptional processes including splicing, RNA maturation, RNA editing and translation (Manna 2015;Cheng et al. 2016b). PPR proteins are RNA-binding proteins that recognize and bind to specific sequences in organellar RNAs (Shikanai and Fujii 2013). The regions bound by PPR proteins are protected from degradation and are found in the small RNA transcriptome (Ruwe et al. 2016). Thus, the ndhB sRNA is a footprint of CRR2 on ndhB transcripts. CRR2 is important for the maturation of the nhdB transcript from the polycistronic RNA rps12-rps7-ndhB transcript (Hashimoto et al. 2003). Loss of function of CRR2 results in the absence of nhdB RNA without much effects on the accumulation of the rps12-rps7 transcript (Hashimoto et al. 2003). The reduction in ndhB RNA levels by Cd treatment but not those of rps12 or rps7, together with the reduction in the levels of the 22 nt ndhB sRNA, suggests that CRR2-mediated ndhB RNA maturation is compromised by Cd stress.
To determine how whether Cd stress might re-program organellar gene expression through nuclear PPR genes in general, we examined the expression profiles of PPR family genes in RNA-seq. Strikingly, 61 and 108 PPR family genes were differentially regulated after Cd stress in TC and YQ leaves, respectively, and nearly all of them were down-regulated at 24 HAC. Perhaps, this suggests that Cd stress led to a compromise in the nuclear control of organellar gene expression. The expression of many PPR genes recovered or even became up-regulated at 48 HAC, indicating a recovery of nucleus-mediated regulation of organellar gene expression. Between YQ and TC, more PPR genes were down-regulated at 24 HAC in YQ, and less recovery of PPR expression was observed at 48 HAC (Fig. 4f), indicating a more severe Cd stress response in YQ leaves.

Differentially expressed miRNAs (DEMs) and their targets in response to Cd stress
As described above, the third most abundant sRNA class by length was the 21 nt sRNA class (Fig. S5), of which 22%-36% were miRNAs (Fig S7). Except in YQ roots, the percentage of 21 nt reads mapping to miRNA regions increased at 48 HAC ( Fig   S7), which may be correlated with the up-regulation of specific 21 nt miRNAs, including bra-miR397-3p, bra-miR397-5p and bra-miR398-3p (Table S7).
Because the miRNA annotations we used were from 2015  and the criteria for plant miRNAs were recently revised by Axtell & Meyers (2018), we performed miRNA discovery with our sRNA-seq data using ShortStack (Axtell 2013) according to the new criteria for plant miRNA annotation. The analysis identified 112 annotated MIR loci meeting the new criteria in at least two independent libraries, and they accounted for approximately 17% of the annotated miRNAs from 2015   (Fig. S8a). Moreover, we identified four new members of known MIR gene families, bra-MIR6028c, bra-MIR403b, bra-MIR6030 and bra-MIR171f S8b). We also identified three novel MIR loci that do not belong to known MIR gene families in miRbase (Fig. S8c). These new MIR family members or novel MIR loci had the potential to form pre-miRNAs with stem-loop structures. In addition, miRNA-5p and miRNA-3p sequences from these loci were both detected in our sRNAseq and together accounted for more than 75% of the reads from the loci (Fig. S8b and c).
As our sRNA-seq only included the leaves and roots of young seedlings and may thus miss miRNAs of low abundance in our annotation, all annotated miRNAs  were used for analysis of differentially expressed miRNAs (DEMs). Except for the leaf samples at 24 HAC, numbers of DEMs were similar between TC and YQ after Cd stress (Fig. S9a). There were more DEMs in roots than in leaves, and the number of overlapping DEMs between TC and YQ was also more in roots than in leaves (Fig. S9b). There were more specific DEMs in YQ leaves at 24 HAC, and this feature was similar to the abundant specific DEGs in YQ leaves at 24 HAC (Fig. S2c), which implied that Cd stress also exerted more evident impact in YQ leaves than TC leaves on miRNAs. We also performed the DEM analysis with only the MIR loci meeting the new criteria. Most of the DEMs were conserved miRNAs (e.g., bra-miR156, bra-miR397 and bra-miR398), but still some non-conserved miRNAs were differentially expressed in response to Cd (Fig. 5 and Fig. 6). For example, bra-miR9560.2-5p (24 nt) was significantly up-regulated in roots and leaves of the two cultivars, and bra-miR5718-3p (22 nt) and bra-miR6202-5p (21 nt) were significantly down-regulated in both TC and YQ roots (Fig. 5 and Fig. 6). Bra-miR403b-3p and -5p from the newly identified MIR403 family member, were identified as up-regulated DEMs in YQ leaves upon Cd stress at 24 HAC, indicating that these non-conserved miRNAs are probably early response factors to Cd stress.
The levels of miRNAs in B. parachinensis seedlings treated by Cd stress were further verified by RNA gel blot analysis ( Fig. 7a and b). Results showed that bra-miR398-3p was strongly up-regulated in leaves and roots under Cd treatment. Cd stress also up-regulated bra-miR397-5p and down-regulated bra-miR156-5p in TC and YQ roots.
In leaves, bra-miR171-3p and bra-miR396-5p were down-regulated by Cd stress and 20 the down-regulation was partially recovered in TC at 48 HAC. The expression profiles of the target genes of the verified miRNAs had changes in the opposite direction compared to their corresponding miRNAs ( Fig. 7c and d). Among them, the target genes of bra-miR398-3p, which encode mavicyanin-like proteins, were obviously down-regulated in both TC and YQ. Laccase (LAC) family members, LAC11-1 and LAC4-2, which are targeted by bra-miR397-5p, were significantly downregulated in YQ roots and slightly down-regulated in TC roots. Three target genes of bra-miR171-3p, encoding the scarecrow-like protein (SCL) family members, were significantly up-regulated in YQ leaves at 24 HAC and in TC leaves at 48 HAC.
Although the target genes of bra-miR156-5p, encoding the SQUAMOSA promoter binding protein-like (SPL) family members, were not identified as DEGs except for SPL9, their transcript levels showed an upward trend under Cd treatment. Regarding the targets of bra-miR396-5p, only one of them was identified as a DEG (BraA01003773, which encodes beta-glucosidase 44, BGL44), and it was downregulated in YQ leaves at 48 HAC. Degradome target-plots (T-plots) further showed that the above-mentioned miRNA target DEGs exhibited a unique and high peak corresponding to the predicted miRNA cleavage site (Fig. S10), thus verifying miRNA-guided cleavage.
To summarize, many miRNAs in B. parachinensis were differentially expressed in response to Cd stress. In particular, bra-miR156-5p, bra-miR398-3p, bra-miR397-5p and bra-miR171-3p were shown to regulate their target genes in response to Cd stress via RNA-seq and degradome sequencing.

Metal ion transporters and Cd accumulation
Cd is a non-essential element for, and toxic to, plants (Lin and Aarts 2012).
However, Cd may occupy the transport channels for other essential metal elements such as iron, zinc and copper because of their similar chemical properties (Sarwar et al. 2010). When B. parachinensis was treated with a high concentration of Cd, some transport-related genes, including metal ion transport genes and transmembrane transport genes, were up-regulated in the two tested cultivars (Fig.   3). Our findings suggest that the up-regulation of metal ion transport related genes is a common response in the two cultivars and may even be a common response in B. parachinensis. In B. napus, Cd stress similarly up-regulates transport-related genes (Zhang et al. 2018b). Up-regulation of these transporters may lead to the non-selective, competitive absorption and transportation of both Cd 2+ and other metal ions and increase Cd accumulation in plants. This model is consistent with the observation that the application of trace element fertilizers can reduce Cd accumulation in plants, which has been proposed as a strategy for the prevention and control of Cd pollution (Sarwar et al. 2010). Thus, blindly mutating individual transporter genes would be a misguided strategy to reduce Cd accumulation in plants. Instead, Cd-specific transporters or specific transporter domains that confer Cd selectivity should be identified and targeted for the breeding of crops with low Cd accumulation. While many metal transporter genes are commonly up-regulated upon Cd stress in the two cultivars, some show different expression profiles. These may help to understand the differential Cd transport efficiencies between the two cultivars. Six genes in the ABCC family (ABCC3-3, ABCC3-1, ABCC1-2, ABCC15, ABCC4-1 and ABCC7) were up-regulated at 24 HAC in YQ leaves but not in TC leaves (Fig. 3). As some members of the ABCC family (e.g., Arabidopsis ABCC3, ABCC1 and ABCC2) serve as transporters of phytochelatin-Cd complexes (Brunetti et al. 2015), the upregulation of these six ABCC family genes may reflect an increase in the transport of the phytochelatin-Cd complexes in YQ leaves (Fig. 8a). A similar phenomenon was also observed for YSL family genes. It was reported that maize YSL1 transports nicotianamine-metal complexes, and could transport Cd at a low rate (Schaaf et al. 2004). Two YSL family genes, YSL1-3 and YSL1-1, were markedly up-regulated in YQ leaves at 24 HAC, which may increase Cd transport in YQ leaves (Fig. 8a). It was reported that some HMA family genes that are specifically expressed in roots (e.g., rice HMA3 ) are involved in the sequestration of Cd in vacuoles in root cells (Ueno et al. 2010). We found that two HMA family genes, HMA2-1 and HMA3, were more highly up-regulated in TC roots than in YQ roots. These two genes may contribute to the lower root-to-shoot Cd translocation efficiency in TC (Fig. 8a).
Due to the different Cd translocation efficiency, the number of DEGs in leaves was larger in the YQ cultivar. The expression profiles also showed that most of the Cdresponsive DEGs, especially those related to basal metabolism, signaling and stimulus response, maintained their expression changes in YQ for a longer time, while most DEGs in TC recovered their normal expression levels at 48 HAC. In addition, some chloroplast genes and PPR genes, which are in the nuclear genome but regulate gene expression in mitochondria and chloroplasts (Manna 2015;Cheng et al. 2016b), were also affected by Cd to a larger extent in YQ leaves (Fig. 4).
These findings indicate that, as compared to TC, YQ exhibits a more extensive and lasting response to Cd stress in leaves.

PPR genes involved in Cd response in B. parachinensis
To our knowledge, an sRNA peak at 22 nt in leaves has never been observed in 23 sRNA-seq studies of other Brassica species, including B. rapa ssp. pekinensis cv.
napus (Jian et al. 2018), and 22 nt sRNAs were not examined in previous sRNA-seq studies of Cd stress responses in plants (Han et al. 2016;Kang et al. 2017;Jian et al. 2018). In the present study, we found a 22 nt ndhB sRNA showing a Cd-induced decrease in abundance in leaves of both B. parachinensis cultivars (Fig. 4). The ndhB sRNA reduction probably reflects reduced maturation of the ndhB transcript from the polycistronic precursor rps12-rps7-ndhB by Cd stress. The PPR gene CRR2 encodes an RNA-binding protein that specifically associates with the rps12-rps7-ndhB RNA through the 22-nt sequence and promotes ndhB maturation from the polycistronic precursor RNA (Hashimoto et al. 2003;Ruwe and Schmitz-Linneweber 2012). In crr2 mutants, abundance of the mature ndhB transcript is drastically reduced but that of rps12-rps7 is unaffected (Hashimoto et al. 2003). We found that Cd stress had a similar effect on the transcripts from the polycistronic precursor as crr2 mutations, suggesting that CRR2's ability to process the polycistronic transcript is compromised by Cd stress (Fig. 8a). While CRR2 expression was not affected by Cd at the transcript level, the reduction in the CRR2 footprint (the 22 nt sRNA) implies that either CRR2 protein level were reduced or the binding of CRR2 to ndhB RNA was compromised by Cd.
PPR proteins are important regulators of RNA editing, maturation, stabilization or intron splicing in chloroplasts and mitochondria (Manna 2015;Cheng et al. 2016b).
Strikingly, we found that many PPR genes were down-regulated upon Cd stress in both cultivars. This implies that Cd exerts a broad effect on organellar gene expression at the post-transcriptional level. We observed that Cd stress had a stronger effect on the expression of some chloroplast genes and nuclear PPR genes in the cultivar with higher Cd translocation efficiency ( Fig. 4e and 4f). The expression of most of down-regulated PPR genes recovered in the cultivar with lower Cd accumulation, but not in the cultivar with higher Cd accumulation. The lasting and strong suppression of PPR genes in the cultivar with higher Cd accumulation further corroborates the hypothesis that Cd compromises the nuclear control of chloroplast genes.
A number of conserved miRNAs were also identified as DEMs. Among them, bra-miR156-5p, bra-miR171-3p and bra-miR396-5p under Cd stress, were consistent with previous findings in rice from microarray analysis (Ding et al. 2011). However, the expression of these conserved miRNAs was not affected by Cd stress in previous sRNA-seq studies of rapeseed (Zhou et al. 2012a) or B. parachinensis (Zhou et al. 2017). In the present study, the differential expression of these miRNAs was verified by RNA gel blot analysis and their confirmed target genes were correspondingly affected.
It is noteworthy that Cd stress led to significant up-regulation of both bra-miR397-5p and bra-miR398-3p and significant down-regulation of their target genes, LAC4-2/11-1 and a mavicyanin-like protein gene, respectively. These changes are similar to the findings in Arabidopsis thaliana (Gielen et al. 2016), and these two miRNAs have also been implicated in copper homeostasis regulation (Ding and Zhu 2009). It has been found that SPL7 and phytochelatin synthase induce miR397 and miR398 levels after Cd stress (Gielen et al. 2016(Gielen et al. , 2017. A possible mechanism of phytochelatin-mediated induction of miR398 is that Cd increases the levels of phytochelatins, which then chelate both Cd and Cu to induce cellular Cu deficiency; and the increase in miR398 is a well-known Cu deficiency response (Gielen et al. 2016(Gielen et al. , 2017. It is also known that MIR397 and MIR398 expression is regulated by SPL7 through promoter binding (Yamasaki et al. 2009). In this study, we found that bra-miR156-5p and its target SPL genes were down-and up-regulated, respectively, by Cd stress (Fig. 5). We therefore propose a regulatory pathway involving the following components: (Cd 2+ /Cu 2+ ) -miR156 -SPL7 -(miR397 / miR398) -(LAC4-2/11-1 / mavicyanin-like protein gene) (Fig. 8b). Phytochelatin-mediated regulation of miR398 probably interfaces with this pathway to form a network of gene expression responsive to Cd stress.
In summary, in response to Cd stress, the expression of B. parachinensis genes related to metal ion transport and transmembrane transport significantly increased.
The chloroplast gene ndhB and PPR genes were significantly reduced in expression by Cd stress. The differential severity of gene expression changes in the two cultivars was correlated to differences in their Cd translocation efficiency. Many miRNAs, including conserved and none-conserved, are likely involved in the initial response to Cd. Among them, miR156, miR397, miR398 and their target genes form a possible regulatory pathway in Cd stress response. These findings help better understand the mechanisms of Cd stress response in plants.