Bioinformatic Exploration of the Targets of Xylem Sap miRNAs in Maize under Cadmium Stress

Cadmium (Cd) has the potential to be chronically toxic to humans through contaminated crop products. MicroRNAs (miRNAs) can move systemically in plants. To investigate the roles of long-distance moving xylem miRNAs in regulating maize response to Cd stress, three xylem sap small RNA (sRNA) libraries were constructed for high-throughput sequencing to identify potential mobile miRNAs in Cd-stressed maize seedlings and their putative targets in maize transcriptomes. In total, about 199 miRNAs (20–22 nucleotides) were identified in xylem sap from maize seedlings, including 97 newly discovered miRNAs and 102 known miRNAs. Among them, 10 miRNAs showed differential expression in xylem sap after 1 h of Cd treatment. Two miRNAs target prediction tools, psRNAtarget (reporting the inhibition pattern of cleavage) and DPMIND (discovering Plant MiRNA-Target Interaction with degradome evidence), were used in combination to identify, via bioinformatics, the targets of 199 significantly expressed miRNAs in maize xylem sap. The integrative results of these two bioinformatic tools suggested that 27 xylem sap miRNAs inhibit 34 genes through cleavage with degradome evidence. Moreover, nearly 300 other genes were also the potential miRNAs cleavable targets without available degradome data support, and the majority of them were enriched in abiotic stress response, cell signaling, transcription regulation, as well as metal handling. These approaches and results not only enhanced our understanding of the Cd-responsive long-distance transported miRNAs from the view of xylem sap, but also provided novel insights for predicting the molecular genetic mechanisms mediated by miRNAs.


Introduction
Heavy metal accumulation in soils is of concern in agricultural production due to the adverse effects on food safety. Cadmium (Cd) is a non-essential element for plants, however it can be absorbed by the roots from the soil and transported to the aboveground parts; thus, it can not only affect the growth and the subsequent productivity of crops, but can also pose a great threat to human health because of its accumulation in the consumable parts of food crops [1][2][3][4]. MicroRNAs (miRNAs) are the most studied 20-to 22-nucleotide non-protein coding RNAs and are at the heart of regulating gene expression in multiple developmental and signaling pathways [5][6][7]. MiRNAs are hypersensitive to different heavy metals, such as Cd, aluminum, and lead in some crop plants (including rice, maize, oilseed rape, and radishes), and mounting evidence has revealed that miRNAs and their targets weaved networks play important regulatory roles in plant adaptation to different heavy metal stresses [4,[8][9][10][11].
In plants, one of the most fascinating aspects of RNA silencing is its mobile nature, and the movement of small RNA (sRNA) molecules can "non-cell-autonomously" orchestrate developmental and stress responses [12,13]. Results obtained with grafting techniques and transient expression systems have shown that sequence-specific short interfering RNAs with a size of 21-24 nucleotides travel to distant organs [14,15]. Phloem exudates contain diverse miRNAs and at least two of them, miR395 and miR399, involved in responses to nutrient availability, are transmitted through grafts, indicating long-distance movement [12,16,17]. Similarly, siRNA signals produced in source or sink tissues move from cell-to-cell and travel long distances via the phloem to apical tissues [18]. Though the long-distance transport of sRNAs was intensively investigated in phloem, small RNAs have also been isolated from the developing xylem of Populus stems, and a majority of these miRNAs have been predicted to target developmental-and stress/defense-related genes, including those associated with the biosynthesis of cell wall metabolites [19]. Arabidopsis miR857, specifically expressed in the vascular tissues of seedlings, is involved in regulating lignin content and consequently morphogenesis of the secondary xylem in by regulating the expression of its target gene LACCASE7 [20].
MiRNAs negatively regulate their target gene expression at transcriptional and post-transcriptional levels by regulating both messenger RNA (mRNA) degradation and translational inhibition based on miRNA/target sequence complementarity [6]. High-throughput degradome sequencing has been successfully established and adapted to validate miRNA splicing targets in a variety of plant species, such as hyperaccumulator Sedum alfredii [21], Populus [22], rice [23,24], soybean [25], canola [5], and maize [26]. Recently, an integrated web-based tool, DPMIND (Degradome-based Plant MiRNA-Target Interaction and Network Database), was developed to scan sRNA targets in multiple plant species [27].
To thoroughly predict the roles of long-distance moving Cd-responsive maize miRNAs, three xylem sap sRNA libraries of Cd-stressed maize were constructed for high-throughput sequencing. Then, an integrative bioinformatic approach composed of psRNATarget and DPMIND was employed to predict potential targets of Cd-responsive xylem miRNAs. Intriguingly, 34 high-confidence cleavable targets for 27 xylem sap miRNAs were identified. Moreover, nearly 300 other genes were also the potential miRNAs cleavable targets, and the majority of them were enriched in abiotic stress response, cell signaling, transcription regulation, as well as metal handling, chelation, and storage. This investigation, therefore, would provide aid to elucidate the molecular genetic mechanisms underlying plant responses to Cd stress from the aspect of mobile miRNAs.

High-Throughput Sequencing of sRNAs in Maize Xylem Sap
To investigate the differences in maize xylem sap sRNA profiles after Cd treatment, we collected xylem sap samples from Cd-treated seedlings, and three sRNA libraries (including two control samples) for sequencing were generated from the maize xylem sap.
We obtained about 2.7, 3.8, and 3.8 M total reads, represented by 0.96, 1.25, and 1.39 M unique sRNA reads, respectively, from the untreated 0 h (C0), untreated 1 h (C1), and Cd-treated 1 h (Cd1) libraries of xylem sap collected at the indicated time-point and treatment. The length distribution of reads showed that the majority of the reads were 20-24 nt in size, which was within the typical size range for Dicer-derived products [28].
To confirm the expression of sRNAs identified by deep sequencing, eight sRNAs (lengths of 19-25 nt) were randomly selected for quantitative real-time RT-PCR (qRT-PCR) analysis, and these contained three miRNAs (PC-3p-33282_23, zma-miR169l-5p, and zma-miR395a-5p_R-1) with lengths of 20-22 nt (Supplementary Table S1). The comparison indicated that the expression patterns of Cd-responsive sRNAs from high-throughput sequencing and qRT-qPCR exhibited a good concordance (Figure 1, Supplementary Table S1), implying the reliability of the sRNA-seq profiling data for the following analysis. of 20-22 nt (Supplementary Table S1). The comparison indicated that the expression patterns of Cd-responsive sRNAs from high-throughput sequencing and qRT-qPCR exhibited a good concordance ( Figure 1, Supplementary Table S1), implying the reliability of the sRNA-seq profiling data for the following analysis. The expression levels of sRNAs were compared between 1 h of Cd treatment (Cd1) and the control C1 samples. Data of qRT-qPCR are means ± SD from three independent biological replicates.

Identification of Xylem Sap Cd-Responsive miRNAs
After length filtration of the sequenced sRNAs, about 199 miRNAs (20-22 nt) were identified in xylem sap from Cd-treated or -untreated maize, including 97 newly discovered and 102 known miRNAs, which were homologous to the sequences in miRBase (Supplementary Table S2). Based on the number of reads (>10 at least in one sample) and MFEI (≥0.85) [29,30], about 20 new high-confidence miRNAs with relatively high expression levels were identified in the three samples (Table 1).  The expression levels of sRNAs were compared between 1 h of Cd treatment (Cd1) and the control C1 samples. Data of qRT-qPCR are means ± SD from three independent biological replicates.

Identification of Xylem Sap Cd-Responsive miRNAs
After length filtration of the sequenced sRNAs, about 199 miRNAs (20-22 nt) were identified in xylem sap from Cd-treated or -untreated maize, including 97 newly discovered and 102 known miRNAs, which were homologous to the sequences in miRBase (Supplementary Table S2). Based on the number of reads (>10 at least in one sample) and MFEI (≥0.85) [29,30], about 20 new high-confidence miRNAs with relatively high expression levels were identified in the three samples (Table 1). MiRNAs detected in C1 and Cd1 libraries were used for differential expression analysis using the stringent criteria (|log 2 Ratio| ≥1, p ≤ 0.05). Finally, 10 miRNAs showed differential expression in xylem sap after 1 h of Cd treatment ( Table 2). Among them, the expressions of three newly identified miRNAs (PC-3p-10246_108, PC-3p-33282_23, and PC-3p-65413_10) was significantly regulated by Cd exposure (p ≤ 0.01). Regarding the 10 Cd-modulated miRNAs, two highly expressed miRNAs (zma-miR169l-5p, zma-miR398a-3p), and two new miRNAs (PC-3p-10246_108 and PC-3p-33282_23) were upregulated by Cd exposure; Cd significantly negatively regulated the remaining miRNAs (Table 2).

Target Predictions of Xylem Sap miRNAs
To better understand the biological functions of long-distance transported miRNAs, 199 significantly expressed miRNAs (p ≤ 0.05, in at least one dataset) from maize xylem sap were identified for target scanning (Supplementary File S1). The putative target sites in maize cDNAs were predicted using two plant sRNA target prediction tools (psRNAtarget and PsRobot).
With the application of psRNAtarget using the inhibition pattern of 'Cleavage', we identified a total of 2184 transcripts from 1436 maize genes, to be the targets of 196 xylem sap miRNAs (Supplementary Table S3). Using PsRobot, we obtained a total of 2514 transcripts from 1774 maize genes to be targets of 172 xylem sap miRNAs (Supplementary Table S4). Through the integration, we identified a total of 493 transcripts from 332 genes as the cleavable targets of 115 miRNAs in the intersection of results from psRNAtarget and PsRobot.

The Function Classification of the Predicted miRNAs Targets
To gain insights into the functionality of the miRNA targets, all of these 493 transcripts were functionally grouped by agriGO [31] and visualized in the candidate pathway networks with MapMan software [32].
Among the genes within the 'TF' group, 11 members of MYB family, two WRKYs, and one AP2-EREBP as well as one bHLH-transcription factor were all the targets of miR159 family members, whereas five NACs and six Homeobox-transcription factors were targeted by miR164s and miR166s, respectively (Table 3, Figure 2, Supplementary Table S3).  With regard to the targets mapped to "Secondary metabolism" category, nine laccases were exclusively coupled to MIR397 family members. With regard to 'Abiotic stress' response, five genes (including two DNAJ proteins and one ERD ortholog) were targeted by five miRNAs individually. Similarly, within the 'metal handling, chelation, and storage' group, two major facilitator proteins and one MATE efflux family protein were uniquely targeted by three miRNAs. However, ZM2G058032 (heavy-metal-associated domain protein) and ZM2G407032 (ABC transporter) were co-targeted by miR399s (Table 3, Figure 2).  With regard to the targets mapped to "Secondary metabolism" category, nine laccases were exclusively coupled to MIR397 family members. With regard to 'Abiotic stress' response, five genes (including two DNAJ proteins and one ERD ortholog) were targeted by five miRNAs individually. Similarly, within the 'metal handling, chelation, and storage' group, two major facilitator proteins and one MATE efflux family protein were uniquely targeted by three miRNAs. However, ZM2G058032 (heavy-metal-associated domain protein) and ZM2G407032 (ABC transporter) were co-targeted by miR399s (Table 3, Figure 2).
The predicted cleavable targets inhibited by miRNAs were inputted into MapMan software (3.6.0RC1) for metabolic pathways analysis (using the framework of Arabidopsis seed-Molecular Networks). The colored boxes indicated the expectation score output by psRNAtarget, and the larger value was shown in deep red.
In addition, several genes located in the pathway network were the potential targets of novel miRNAs. Particularly in the 'signaling' category, three genes involved in Ca signaling (ZM2G107575, ZM2G312661, and ZM2G174315) and the protein kinase (ZM2G100454) were uniformly targeted by novel miRNAs (Table 3, Figure 2).

The miRNAs Cleavable Targets
The transcripts of 332 maize genes in the intersection of psRNAtarget and PsRobot output, which were the potential 115 miRNAs cleavable targets, were further evidenced by the maize degradome data. With the aid of the DPMIND webserver [27], which harbors the degradome sequencing data of maize anther and ears, we obtained 34 maize genes harbored target sites of 27 xylem sap miRNAs by combining the outputs of psRNAtarget, PsRobot, and DPMIND (Table 4,  Supplementary Table S5).
Regarding these cleavable candidates, the majority of them were intensively studied, including 11 squamosa promoter binding proteins (SBPs), seven nuclear transcription factor Ys (NFY), and three auxin response transcription factors (ARF) [6]. After excluding these well-known targets of miRNAs, several genes were filtered out as fresh cleavable targets for these xylem sap miRNAs, including the common stress-responsive miRNAs and their targets [6].
Many miRNAs appear to function together via co-targeting to regulate functionally related genes or pathways, and vice versa [33]. For example, the myb74 transcription factor (ZM2G028054_T03) was co-modulated by zma-miR159a-3p_R-1 and zma-miR319a-3p_R+1, while myb138 (ZM2G139688_T01) was specifically targeted by zma-miR159a-3p_R-1 (Table 4). In contrast, two F-box proteins (ZM2G064954_T01 and ZM2G119650_T01) were both the targets of zma-miR394a-5p. Similarly, ZM2G155490_T01 and ZM2G304745_T01 (encoding LRR receptor-like kinase) were targeted by zma-miR390a-5p (Table 4).  The miRNAs in italic and those maize transcripts between parallel lines mean each of the miRNAs can target each transcript successively. psRNAtarget output: Exp for Maximum expectation and the star (*) indicating the largest score of the corresponding miRNA-target combinations, unpaired energy (UPE) for maximum energy to unpaired target site; For DPMIND, miR represents the homolog of the queried miRNA by BLAST, and the dollar label ( $ ) means the least number of the degradome datasets for the miR-target associations.

The Long-Distance Transport of miRNAs
Long-distance transport of signaling molecules, intensively investigated in phloem, is known to be a major component in plant growth regulation, as well as their adaptation to changing environmental conditions [34]. Although some studies have demonstrated the presence of numerous miRNAs in phloem tissues [35][36][37], so far, only three miRNAs (miR399, miR395, and miR172) have been shown to move long distances in plants [16,17,34,38].
However, xylem also plays an important role in the root-to-shoot signaling system [39]. Furthermore, the root-to-shoot Cd translocation process may be more complex than previously thought [40]. Regarding xylem, research on miRNAs is in its infancy. In this study, about 199 miRNAs were identified in Cd-treated or -untreated maize xylem sap, including 97 novel and 102 known maize miRNAs. Moreover, the three famous long-distance moving miRNAs (miR399, miR395, and miR172) were detected in maize xylem sap ( Table 3, Supplementary Table S2). These findings suggested that these miRNAs are potential signal molecules that move systemically.
MATE transporters were involved in the cellular transport and detoxification of Cd [41]. Furthermore, the distribution of allocated transcripts (e.g., MATE) along the root-to-shoot axis was correlated with the siRNA signal spread in hetero-grafted Arabidopsis [42]. Here, we identified zma-miR528a-5p in xylem sap, and predicted ZM2G148937 (MATE family protein) as its potential cleavable target (Table 3).
Six highly conserved amg-miRNA families (amg-miR166, amg-miR172, amg-miR168, amg-miR159, amg-miR394, and amg-miR156) were viewed as potential regulatory sequences of secondary cell wall biosynthesis [43], and Populus Pto-MIR156c might play vital roles in the regulation of wood formation in trees [44]. Moreover, the knockdown of rice MicroRNA166 confers drought resistance by causing leaf rolling and altering stem xylem development [5], and rice miR166 also plays a critical role in Cd accumulation and tolerance [45]. In this study, the miRNAs variants of these six families and nine miR166 isoforms with more than 30 reads in each of the three samples were identified in maize xylem sap (Supplementary Table S2).

The Potential Cleavable Targets of miRNAs in Xylem Sap
In plants, miRNAs and their targets show a pattern of near complementarity, suggesting that plant miRNAs likely act through endonucleolytic cleavage of target mRNAs [46]. In this study, 34 targets were predicted to be inhibited by miRNAs through cleavage with degradome evidence, and most of them were intensively studied targets of common stress-responsive miRNAs, such as NFY, ARFs, SBPs, and GRAS family transcription factors [6]. Concerning the cleavable candidates of xylem sap miRNAs, they were concentrated on the NFY-, ARF-, SBP-, and GRAS-transcription factors, which was targeted by zma-miR169s, zma-miR160f-5p, zma-miR156s, and zma-miR171s, respectively (Table 4).
It is of particular interest to note that three homeobox-transcription factors were the cleavable targets of zma-miR166s, and two of them were annotated as rolled leaf genes (Table 4), which was reminiscent of their rice ortholog OsHB4. Moreover, miR166 plays a critical role in Cd tolerance as well as in drought resistance through regulation of its cleavable HD-Zip target gene OsHB4 in rice [5,45]. Altogether, these results of miR166-mediated regulatory cascade strengthened the pivotal role miR166-HB couple in abiotic stress acclimation. Moreover, among the miRNA targets within the 'TF' group, members of other TF families (e.g., MYBs, WRKYs, NACs) were the specific targets of certain miRNA family members (Table 3, Figure 2, Supplementary Table S3). These TF-type miRNA targets might contribute to elucidate the complicated mechanism of Cd stress from the aspect of transcriptional network through unveiling transcription factor-regulated downstream target genes which were involved in Cd stress response.
In addition to these well-known targets of miRNAs, several genes were filtered out as novel cleavable targets for these xylem sap miRNAs, including the common stress-responsive miRNAs and their targets [6] (Table 4).
Many miRNAs appear to function together via co-targeting to regulate functionally related genes, and vice versa [33]. For example, the myb74 was co-modulated by zma-miR159a-3p_R-1 and zma-miR319a-3p_R+1, whereas two F-box proteins were both the targets of zma-miR394a-5p (Table 4). These co-targeting paradigms indicated that individual miRNA variants have different functions according to their specific targets [33].
Besides the 34 maize genes as the cleavable candidates of miRNAs with degradome evidence from the limited degradome data (Table 4), nearly 300 other genes were also the potential miRNAs cleavable targets ( Table 3, Supplementary Table S3). From a global view, many target genes were prone to be enriched in abiotic stress response, cell signaling, transcription regulation, as well as metal handling, chelation, and storage (Table 3, Figure 2). Using a degradome sequencing approach, leucine-rich repeat (LRR) protein, cation transporting ATPase, and Myb transcription factors, were found to be cleaved by miRNAs under heavy metal stress [25]. Furthermore, a few miRNA cleavable targets, including iron transporter and ABC transporter, were involved in plant responses to Cd stress [2,8]. MATE transporters were involved in the cellular transport and detoxification of Cd [41]. In this study, ZM2G148937 (MATE family protein) was predicted as the cleavable target of zma-miR528a-5p (Table 3). From another perspective, these targets in relation to the Cd stress acclimation highlighted the role of the corresponding miRNAs in regulating Cd stress response.
Intriguingly, these known cleavable targets of miRNAs were also identified in this investigation (Table 3), though we did not retrieve degradome evidence from the available public dataset. miRNAs were regarded as a new target for genetically improving plant tolerance to certain stresses [6]. From the perspective of miRNA-target couple, the characterization of the miRNAs and the associated targets in responses to Cd exposure will provide a framework for understanding the molecular mechanism of heavy metal tolerance in plants. Thus, it would be interesting to determine the role of these long-transported miRNAs, and whether these xylem sap miRNAs are transported to leaves under heavy metal stress or other stresses. Future investigation on the final location of xylem sap miRNAs, which might be achieved through analyzing the difference between the exudates from the node incisions where ears or leaves separated and xylem sap below the node, together with the target identification through degradome and proteome or ribosome profiling of the detached ears or leaves, will aid to illustrate the effect of xylem sap miRNAs on their potential targets in their final destination.

Plant Materials and Cd Treatment
The seedlings of maize (Zea mays L. cv. Nongda 108; China) were cultivated in a hydroponic system in a growth chamber with a temperature of 22 • C (night) to 28 • C (day), photosynthetic active radiation of 200 µmol·m −2 ·s −1 , and a 14/10-h day/night photoperiod. All hydroponic solutions were continuously aerated and renewed every three days. When the third leaves were fully expanded, the seedlings were transferred into fresh growing solutions containing 100 uM CdCl 2 , according to previous reports [2,3].

Sampling of Xylem Sap
The seedlings were separated into three groups-untreated 0 h (C0), untreated 1 h (C1), and Cd-treated 1 h (Cd1). For each group, 30 maize seedlings at the indicated timepoint/treatment were de-topped by cutting the stem with a razor blade just above the first internode, and the remaining parts without stem were used for xylem sap collection, according to previous reports with minor modifications [47][48][49][50][51]. Then, the cut surface was rinsed twice with distilled water and the liquid drawn in the first 5 min was discarded. Finally, the bleeding saps from 30 maize plants were harvested with 10 uL syringe for 1 h after cutting and mixed in a tube containing 1 mL Trizol as one sample replication, and single biological replication for each sample was used for sRNA sequencing.

Small RNA Library Preparation and Sequencing
Total RNA was extracted from C0, C1, and Cd1 xylem sap samples using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's procedure. Approximately 1 µg of total RNA were used to prepare a small RNA library (single biological replication for each sample) according to protocol of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, USA). Then, we performed the single-end sequencing (36 bp) on an Illumina Hiseq2500 at the LC-BIO (Hangzhou, China) following the vendor's recommended protocol. The data was uploaded to NCBI/SRA with accession number SRP073229 (https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP073229), and contained miRNAs reads of C0, C1, and Cd1 xylem sap samples (single biological replication for each sample).

Small RNA Analysis
Data processing followed the procedures as described previously [23,24]. The raw reads were subjected to the Illumina pipeline filter (Solexa 0.3), and then the dataset was further processed with an in-house program, ACGT101-V4.2 (LC Sciences, Houston, TX, USA) to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA), and repeats.

Identification of Known and Novel miRNAs
Unique sequences with lengths of 20-22 nucleotide [52] were mapped to monocot plant precursors in miRBase Release 22.1 (October 2018, mirbase.org) by BLAST search to identify known miRNAs and novel 3p-and 5p-derived miRNAs.
The unique sequences mapping to maize mature miRNAs in hairpin arms were identified as known miRNAs. The unique sequences mapping to the other arm of known maize precursor hairpin opposite to the annotated mature miRNA-containing arm were considered to be novel 5p-or 3p-derived miRNA candidates. The remaining sequences were mapped to other monocot precursors (with the exclusion of maize) in miRBase 20.0 by BLAST search, and the mapped pre-miRNAs were further BLASTed against the maize genomes (ftp://ftp.maizesequence.org/pub/maize/release-5b/assembly/ ZmB73_RefGen_v2.tar.gz with genome annotation file ZmB73_5b_FGS.gff.gz) to determine their genomic locations [53].
The sequences unmapped to maize mature miRNAs or other monocot miRNAs precursors were further BLASTed against the maize genomes, and the hairpin RNA structures containing mappable sequences were predicted from the flank 120 nt sequences using RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi).The criteria used to annotate potential miRNAs was as described previously [54][55][56][57]. Sequences that met the criteria for secondary structure prediction were then considered to be novel miRNA precursors [53,58,59]. Minimum free energy index (MFEI) was also taken into account for evaluating the confidence of novel miRNAs, with the expected value (≥0.85) [7,29,30]. Moreover, all of the aforementioned criteria must be fulfilled in at least two distinct sRNA-seq libraries [52].

Differential Expression Analysis of sRNAs Under Cd Stress
Data normalization followed the procedures as described in a previous study [60]. sRNA differential expression based on normalized deep-sequencing counts was analyzed by selectively using Fisher's exact test and Chi-squared 2 × 2 test with the selection threshold of 0.05 [23,28,60].
To investigate the differentially expressed miRNAs between libraries, we compared the gene expression patterns of miRNAs in Cd1 and C1 library. Towards this purpose, we considered the following criteria: (1) p-value should be less than 0.05 (p ≤ 0.05) in at least one dataset; and (2) Log 2 ratio of fold change between normalized counts of C1 and Cd1 libraries was greater than 1 or less than −1 [28,56,57,59,61].

The Prediction of miRNA Targets
After filtering, we get 199 significantly expressed miRNAs (p ≤ 0.05, in at least one of the three samples). Then, these xylem sap 199 miRNAs were used to interrogate maize annotated cDNAs sequences (Zea_mays.AGPv3.22, ftp://ftp.ensemblgenomes.org) preloaded in the psRNAtarget web server (plantgrn.noble.org/psRNATarget) [62] for predicting target sites as described previously [61], and the default criteria for target prediction in psRNAtarget website were used.
Prediction of miRNA target genes was also performed by the local version of psRobot_tar scripts in psRobot (omicslab.genetics.ac.cn/psRobot) [63] using the 199 miRNAs for scanning maize cDNAs sequences (http://ftp.maizesequence.org/release-5a/filtered-set/ZmB73_5b_FGS_cdna.fasta.gz) with default settings (moderate model of pre-set the target penalty score 2.5, thus alignments that meet the penalty score cutoff will be reported in the result).
To scan the cleavable targets of the identified known miRNAs and novel miRNAs, these miRNAs were uploaded to DPMIND (http://cbi.njau.edu.cn/DPMIND) [27] to locate the homologous miRNAs for the following degradome data query.
For determining the expression of sRNAs, about 2 µg RNAs were reverse-transcribed by miRcute miRNA First-Strand cDNA Synthesis kit (TIANGEN, Beijing, China). Transcript levels of mature sRNAs were measured by qRT-PCR using a DNA Engine Opticon 2 real-time PCR detection system (Bio-Rad, Hercules, CA, USA) with miRcute miRNA qPCR Detection kit (TIANGEN) according to the manufacturer's instructions. Details of the primers used are listed in Supplementary Table S1. The maize 5S RNA was used as the internal control for RNA template normalization [64]. All reactions were run in triplicate. The relative expression levels of sRNAs were calculated by the comparative threshold cycle (Ct) method. At least three independent biological replicates were used for each small RNA.