Data were obtained by two methods: 1) all publicly available (NCBI) sRNA and degradome libraries from A. thaliana were compiled to provide a diverse representation of sRNA expression and regulation under varied conditions; these were derived from multiple tissues, developmental stages, and biotic and abiotic stress conditions; and 2) paired sets of sRNA-Seq and degradome data from aliquots of individual RNA extracts were produced as part of this study for 14 independent plant samples. There were four plant-virus treatments (described below) with four biological replicates per treatment for both sRNA-Seq and degradome analyses; two sRNA-seq libraries were low quality and removed. All the sRNA and degradome data from (1) and (2) were combined to identify an sRNA-mediated regulatory meta-network (described below).

Two-week-old A. thaliana Col. plants grown at 22 C with a 10 h photoperiod were mechanically inoculated with Cucumber mosaic virus or rubbed without virus as mock controls. (These treatments are from a separate study and the effect of virus is not addressed in this report.) Leaf tissue was collected 10 days post-inoculation, ground in liquid nitrogen, and total RNA extracted using Trizol (Thermo-Fisher) as recommended by the manufacturer. Each resulting RNA preparation was divided into two aliquots to be used as input for sRNA-Seq and degradome libraries. sRNA libraries were prepared from 1 μg of total RNA using methods described previously (Vargas-Asencio et al., 2017). For the degradome libraries, ~40 μg of total RNA was used. Degradome libraries were constructed using the method described by Zhai et al. (2014), but with the following modifications: a) different adapters and primer sequences were used ( Additional file 1: Table S1 ), b) the PCR clean-up step was performed using Axygen™ AxyPrep Mag™ PCR Clean-up (Fisher) instead of Agencourt AMPure XP beads (Beckman Coulter), and c) EcoP151 (NEB) was used for the restriction enzyme digestion step instead of MmeI. Sequencing was performed using an Illumina Hiseq 4000 at the Genomics Resources Core Facility, Weill Cornell, NY, to obtain single-end 51-nt reads for both sRNA (accessions: SRR6234880- SRR6234893) and degradome libraries (accessions: SRR6235006- SRR6235021).

A custom bioinformatics pipeline was implemented to identify sRNA-mediated networks. A detailed description is provided in the following sections. The overall strategy was to gather all available sRNA and degradome data, and to combine it with existing genome annotations and sRNA databases to produce a data-driven, degradome-supported network of interactions between sRNAs and transcripts. There are two types of nodes in the proposed network: sRNAs and transcripts. sRNAs include miRNA and phasiRNAs, and transcripts include miRNA precursors, PHAS loci, and mRNA transcripts targeted by sRNAs. Annotations are available for miRNAs, miRNA precursors, and potential target transcripts, while for PHAS loci, their sRNA triggers, and the resulting phasiRNAs, there are no genome-wide annotation available. The identification of these components and their interactions was therefore part of the tasks included in the pipeline. Newly generated annotations were combined with available genome and known miRNA annotations to perform a genome-wide-level search for sRNA–target interactions. Once all components and their interactions were identified and experimentally validated, they were consolidated into a network for downstream analysis.

The TAIR10 version for A. thaliana provided the reference genome (Swarbreck et al., 2008). Genome annotations were obtained from Araport11 (Cheng et al., 2017). Known miRNA and precursor sequences were obtained from miRBase (Kozomara and Griffiths-Jones, 2014) release 22. Gene ontology terms were obtained from Ensembl Genomes release 37 (Kersey et al., 2017).

Fourteen sRNA and 16 degradome libraries were produced in this study. These data were complemented with all publicly available sRNA datasets representing different tissues, stress conditions (biotic and abiotic), and developmental stages in A. thaliana ( Additional file 2: Table S2 ), as well as all available degradome datasets ( Additional file 3: Table S3 ).

Reads (51 nt) from sRNA-Seq libraries were filtered using the adaptive adapter trimming function in Trim Galore (Kruger) to account for variability in library construction methodologies. Size was constrained to 20–40 nt after adapter trimming, and non-adapter containing reads were removed. Datasets were collapsed to unique sequences using the Fastx toolkit (Hannon); sequences with fewer than 50 reads were removed. Libraries containing less than 100 unique sequences were considered non-informative and removed. SRA degradome libraries were filtered using the adaptive adapter trimming function in Trim Galore with the minimum size after adapter trimming set to 18 nt. The resulting libraries were evaluated manually, and additional trimming was performed if there was evidence of remaining adapter sequences. For the libraries produced in this study, the first 6 nt derived from the library preparation process were removed. The Fastx toolkit was used to convert reads to fasta format.

PHAS loci detection was performed for each dataset using PhaseTank (Guo et al., 2015). Locus extension was set to zero, and the top 15% of regions with the highest accumulation of mapped reads (described as relative small RNA production regions in Guo et al., 2015) were analyzed for phasiRNA production. Results for all datasets were combined to produce PHAS loci with maximum length from overlapped results. Potential PHAS loci detected in less than 3 of the 902 libraries were discarded. The resulting loci were then extended by 220 nt on each side to perform a search for sRNA triggers associated with phasiRNA production.

PhasiRNA production triggers were searched using the degradome data. Thirty-nine degradome libraries were independently analyzed using CleaveLand4 (Addo-Quaye et al., 2009). Sequences from both strands of the extended PHAS loci were evaluated using known miRNAs as queries. A weighted scoring system (deg_score) to compile the independent degradome analysis results was developed as follows: cleavage events with degradome category zero per CleaveLand4 were given a score of 5, cleavage events with degradome category one were given a score of 4, cleavage events with degradome category two were given a score of 0.5. The scores for each event were added across all 39 degradome libraries. The highest scoring event per PHAS locus was selected as the initial phasiRNA triggering site; a minimum score of 10 was set to assigned triggers. When triggers were found, the polarity of the loci was set to the strand complementary to the triggers.

To identify the phasiRNAs produced by each PHAS locus sRNA reads from each library were mapped to the extended PHAS loci independently. No mismatches were allowed, sRNAs of 21 and 22 nt were accepted, counts for reads mapping to multiple locations were divided between the number of locations, reads with more than 10 mapping locations were removed, and reads mapping outside the original region (before extension) were not considered. Mapped reads were assigned to bins from 1 to 21 (phases) according to their mapping positions from the 5' end. Positions of reverse reads were shifted (+2) due to 3' overhang, to match forward read bin positions. The mapping was performed on each strand of the PHAS loci independently. A scoring system was developed to rank bins by read abundance for each locus across all sRNA libraries. The three most abundant bins per locus per library were used. The most abundant bin was given a score of 5, the second most abundant was given a score of 2, and the third most abundant was given a score of 0.5. The resulting scores from all libraries were added for each bin to produce a ranking of sRNA bins for each PHAS locus.

PhasiRNAs derived from miRNA triggering events were found by matching the phase register set by the degradome-confirmed miRNA triggering events to bin assignments. PhasiRNAs from immediately adjacent bins (-1, +1) were also collected. For the PHAS loci where no trigger was found, sRNAs from the most abundant bin and immediately adjacent bins were collected.

Resulting phasiRNAs were pooled with all known miRNAs to produce a new set of queries to search for phasiRNA production triggers using the degradome-based ranking strategy described above. To identify secondary and/or tertiary triggers, sRNAs whose cleavage events matched the polarity of the primary trigger (highest ranked, with score >10) were kept. The potential secondary/tertiary triggers were evaluated by matching their slicing site coordinates to those corresponding to the three most abundant sRNA bins per PHAS locus. Because 22-nt sRNAs were included in the analysis, which can alter the 21 nt phasing, the bins immediately adjacent (-1,+1) were also considered. In the cases where a match was found, the sRNAs were considered additional phasiRNA triggers. The assignment of secondary/tertiary triggers was further evaluated by determining if the phasiRNAs contained in the matched bins were biologically active (described below). PhasiRNAs derived from secondary and/or tertiary sRNA triggering events were found by matching the phase register set by the degradome-derived sRNA triggering events to bin assignments. The resulting phasiRNAs were pooled with known miRNAs to produce a final set of queries to search for phasiRNA production triggers using the strategy described above in this paragraph.

Corresponding trigger, PHAS locus and phasiRNA sets were evaluated and confirmed manually to produce a miRNA-PHAS loci-phasiRNA annotation. A novel nomenclature is proposed for phasiRNAs in order to provide consistent and detailed information about their biogenesis. To assign a PHAS loci to a gene ID, the PHAS loci with polarity assigned based on confirmed sRNA triggers were compared to the araport11 genome annotation, and if the locus had significant overlap (>70%) and matching polarity to annotated features (genes, transposons), the locus was assigned to the feature. If more than one feature matched a locus. If no trigger was found but the PHAS locus overlapped with an annotation, a tentative assignation notated with lowercase was used; if the PHAS locus did not match any annotation, the forward genomic orientation was kept and the loci were named using their coordinates. For phasiRNAs, they were named using the PHAS locus from which they derived, followed by up to four descriptors: 1) the number of registers (21 nt) from the 5' end of the transcript; 2) in parenthesis, offset to main phased register, if any; 3) polarity, a “+” was used if the phasiRNA derived from the mRNA strand or “-” if derived from the complementary sense strand; and 4) size was indicated in the case of 22 nt long phasiRNAs by adding “_22” to the end.

To determine if the selection of canonical and non-canonical phased registers within phasi loci was adequate, all sRNAs (>50 copies per library) of 21 and 22 nt produced in this study that mapped to the A. thaliana genome were mapped to the regions where PHAS loci were detected. The mapped reads were evaluated according to whether their position corresponded to sites described in the resulting annotation.

To evaluate the improvement in identification of sRNAs from sRNA-Seq datasets, all sRNAs of presumed biological relevance (e.g., with >50 copies per library) were identified using the annotation produced in this study. The relative abundance of the distinct types of sRNAs under consideration was evaluated based on the abundance of unique and total reads.

The miRNA-PHAS loci-phasiRNA annotation was used to identify and quantify miRNAs and phasiRNAs as described above; an arbitrary threshold of 50 combined raw count was established to select candidates for transcript targets. Degradome datasets were analyzed independently using CleaveLand4 (Addo-Quaye et al., 2009) to find target transcripts for selected sRNAs. A custom scoring system (target_deg_score) was developed to evaluate the confidence and repeatability of sRNA–target transcript interactions. The following weighted scores were assigned to the degradome categories described in CleaveLand4 (Addo-Quaye et al., 2009): category 0 hits were given a score of 5, category 1 were given a score of 4, and category 2 were given a score of 0.5. Categories 3 and 4 were not considered informative. The scores were summed across all libraries for each cleavage event detected. An empirical cumulative distribution analysis was performed for the target_deg_scores and an arbitrary threshold of 15 was established to select for the 1% most reliable (high quality, most repeatable) sRNA–target transcript interactions.

For the curated set of sRNA–transcript pairs, the location of the target site in the transcript was converted to genomic coordinates using the GET map/cdna/:id/:region API from Ensembl Genomes (Kersey et al., 2017). Genomic features overlapping the resulting genomic coordinates were extracted and tabulated according to their abundance.

For the sRNA-mediated regulatory meta-network, a bipartite directed meta-network was constructed by obtaining sRNAs, transcripts, and their interactions from a combination of existing miRNA and transcript annotations, the newly developed miRNA-PHAS loci-phasiRNA annotation, and the degradome search results from all datasets. Cytoscape (Shannon et al., 2003) was used for visualization and structural analysis. A functional characterization was performed using Bingo (Maere et al., 2005) to determine the representation of GO Slim categories of the genes included in the network as compared to the genome as a whole.

Currently there are no genomic features or sequence signatures that allow the identification of PHAS loci (regions of phasiRNA production) from a genome sequence; instead their detection depends on a search for phased patterns in sRNA-Seq data. Because of the observed variability in size between sRNA libraries and the assumption that phasiRNA expression may depend on specific environmental queues, the strategy used to identify PHAS loci was to independently evaluate sRNA-Seq datasets and then merge overlapping results to produce PHAS loci with maximum length. The number of consensus PHAS loci detected was variable, ranging from zero to more than 120 per library ( Figure 1A ).

A total of 942 PHAS loci were identified from the combined libraries (n=902; Additional file 4: Table S4 ). The consistency of PHAS loci detection was evaluated by determining the number of recognition events for each locus across all libraries. To remove spurious results, only PHAS loci detected in at least three libraries were included. The number of recognition events varied, with 107 PHAS loci independently detected in at least three libraries ( Figure 1B ). A failure to detect any given locus in a specific library could be due to expression limited to specific experiment conditions, e.g., stress, developmental stage, or tissue type or to a limitation in sensitivity.

Once PHAS loci had been identified, a recursive method was designed to identify triggers by extending a search up and downstream of the detected PHAS loci, followed by searches for secondary or tertiary triggers that would explain the production of non-canonical phasiRNAs (22 nt long or derived from an alternative phased register). sRNAs were considered triggers if their predicted targeted position on the PHAS locus was consistent with the distribution of sRNAs and was supported by degradome data according to the quantitative criteria (deg_score) detailed in Materials and Methods. For the 107 PHAS loci evaluated, triggers were assigned for 57 of them. From the 108 unique sRNAs triggers identified, 16 corresponded to miRNAs and 92 were phasiRNAs; in some cases, sRNAs were assigned to multiple PHAS loci ( Additional file 5: Table S5 ). Among the triggers, there were 16 canonical phasiRNA; 45 were 22 nt long, 64 were from a secondary phased register; 33 were both 22 nt long and derived from a secondary phased register. Consistent with Rajeswaran et al. (2012), multiple triggers per PHAS locus were detected in some cases ( Figure 1C ). The length of the PHAS locus was positively correlated with the number of putative triggers (R² = 0.59, t = 5.54, df = 55, p < 0.01, Figure 1D ); however, the magnitude of degradome support for putative triggers only showed a weak negative correlation with PHAS locus length (R²= -0.199, t = -2.8093, df = 190, p < 0.01, Figure 1E ). Together these results indicate that even though more putative triggers were found for larger PHAS loci, the degradome support for sRNA triggers in larger PHAS loci was essentially retained.

To identify the genes or genome features corresponding to PHAS loci, regions containing PHAS loci were compared to the Araport11 (Cheng et al., 2017) genome annotation. An sRNA trigger acts on the transcribed strand of the PHAS locus, and when available, the polarity of the transcripts was inferred from the trigger assignments described in Additional file 5: Table S5 . Knowledge of the transcript's polarity increased the confidence in the identification of the corresponding gene or feature and helped to resolve situations where annotated features existed for both genomic strands. Ninety-two of the 107 PHAS loci overlapped with annotated genomic locations, of which 49 were specifically assigned to genes based on their location and inferred strand polarity ( Table 1 ). Six of the PHAS regions overlapped with two genomic features (when genes are present at the + and - strands) for a total of 98 genomic features that overlapped with the identified PHAS loci. These included: i) 41 PHAS loci overlapping with gene types previously shown to be involved in phasiRNA production, including six TAS genes, 20 (one pseudogene) Pentatricopeptide repeat/Tetratricopeptide repeat (PPR/TPR) protein genes, six F-box containing /Auxin response genes (ARFs), five (one pseudogene) nucleotide binding site-leucine rich repeat (NBS-LRR) domain containing genes (two PHAS loci were found in AT5G38850) (Wang and Chekanova, 2016), AGO1, a cation exchanger, a basic helix-loop-helix (bHLH) transcription factor and a Curculin-like (mannose-binding) lectin family protein (Chen et al., 2007; Xia et al., 2015); ii) six genes involved in gene expression regulation including DCL1, three GRAS family transcription factors, one LOB domain-containing transcription factor, and a DNA methyl transferase were found to produce phasiRNAs; iii) 24 genes involved in metabolism, structure, and other functions, with no previous connections to phasiRNAs; iv) 14 detected PHAS loci overlapping with 12 transposable elements (three PHAS loci were found in AT1TE51040); and v) a mixture of 13 overlapping sites, including five annotated as “long non-coding RNA,” four “Novel transcribed regions,” one “Other RNA,” and three with “natural antisense transcript”; these regions of phasiRNA production can be re-annotated as PHAS loci based on these results. There were 15 additional PHAS loci located in unannotated regions of the A. thaliana genome.

PhasiRNAs have the capacity to act as triggers to induce phasiRNA production in cis and/or trans by targeting their progenitor transcripts or unrelated transcripts, respectively, leading to regulatory modules that can involve multiple PHAS loci with the potential to integrate the regulation of phasiRNA production. Though a number of interactions of this kind (referred to herein as high-level interactions) have been described (Rajeswaran et al., 2012; Xia et al., 2017), a genome-level evaluation is not available.

To better understand these high-level interactions between PHAS loci, regulatory modules were constructed from PHAS loci connected by the phasiRNA triggering events tabulated in Additional file 5: Table S5 . Seventeen regulatory modules of varying size and complexity were identified and are shown in Figure 2 . Specific details on any component of the network can be found by cross-referencing Figure 2 with the tabulated information in Additional file 5: Table S5 . Consistent with the validity of the methodological approaches taken, the results included previously described interactions such as miR173-TAS1-TAS2-PPR/TPR (Allen et al., 2005; Howell et al., 2007), miR161-PPR (Howell et al., 2007), mir390-TAS3-ARF (Ozerova et al., 2013; Xia et al., 2017), miR393-AFB (Si-Ammour et al., 2011), miR828-TAS4 (Rajagopalan et al., 2006), miR472-NBS-LRR (Boccara et al., 2014), and miR168-AGO1 (Mallory and Vaucheret, 2009). The genome-wide approach used in this study allowed an extension of some of these modules. The miR173-TAS1/2-PPR/TPR was found to include seven additional PPR/TPR genes, a connection to miR161.2 and to interact with phasiRNAs-producing non-coding RNAs. A Basic helix-loop-helix (bHLH) type transcription factor and one non-coding RNA were linked to the miR393-AFB-F-box module, suggesting the existence of complex hormone-sRNA-TF regulatory interactions. The miR168-AGO1 (Mallory and Vaucheret, 2009) was extended to include interactions with non-coding RNAs, chlorophyll A/B binding proteins, and Proline-rich proteins 4. The miR2939-F-box module makes the third example of F-box containing regulatory modules, indicating an important connection between hormone and sRNA-mediated regulation. Two previously undescribed modules involving important transcriptional level regulators were identified. miR170-miR171-GRAS TF-F-box-LOB and miR773-DNA-MET provide further evidence for a direct connection between transcriptional and post-transcriptional regulation. Eight additional modules included combinations of non-coding RNAs, unannotated regions, and other genes not previously reported to produce phasiRNAs.

The identification of phased registers in PHAS loci is complicated by the possibility of multiple phasiRNA producing triggers and the presence of shift-inducing 22 nt long phasiRNAs. To evaluate the accuracy in phased register selection within PHAS loci, the relative abundance of reads mapping to the PHAS loci regions that match the selected phased registers was evaluated; assuming the abundance of sRNAs is indicative of biological activity, a high proportion of recalled sRNAs would indicate the correct registers were selected. Over 75% of all reads mapping to regions where PHAS loci were detected belong to the registers included in the annotation ( Figure 3A ). These results indicate that the annotation produced from a composite analysis of all available sRNAs libraries in A. thaliana resulted in an adequate representation of PHAS regions that can be used to successfully identify phasiRNAs in individual libraries.

Multiple classes of sRNAs are known in A. thaliana, but an adequate annotation is only available for miRNAs (n=428); phasiRNAs have been reported sparsely, with Rajeswaran et al. (2012) providing the only case where phasiRNAs derived from the TAS1c locus are reported with details of their origin and phased register allowing their identification. Using the resulting miRNA-PHAS loci-phasiRNA annotation from our bioinformatics analysis, ~30 to 50% of all unique 21 and 22 nt long sRNAs from the sRNA-Seq datasets produced in this study were identified ( Figure 3B ); this nearly doubles in most cases the number of unique sequences identified, though a significant number remained unidentified. Evaluation of the total reads ( Figure 3C ) indicates that miRNAs make up the large majority of sRNA reads. There are ~10% of reads corresponding to phasiRNAs, both canonical and non-canonical, and ~15% of reads that remained unassigned.

In plants, sRNAs act mainly through cleavage of their transcripts, yet there are examples of other mechanisms such as translational repression (Borges and Martienssen, 2015; Wang and Chekanova, 2016). Also, it has been shown that not all phasiRNAs produced from a PHAS locus are active; instead only some of them appear to be competent for loading into argonaute (AGO) containing complexes where they exert their activities (Fei et al., 2013). Therefore, in this study known miRNAs and phasiRNAs derived from the detected PHAS loci (including non-canonical phasiRNAs) were evaluated for biological activity using degradome data.

For A. thaliana, a limited number of degradome libraries are publicly available ( Additional file 3: Table S3 ), including 11 datasets corresponding to inflorescence tissue, 6 to leaf tissue, 5 to seedling tissue, and 1 whole plant (Addo-Quaye et al., 2008; Creasey et al., 2014; Thatcher et al., 2015; Hou et al., 2016; Lin et al., 2017). Sixteen new degradome libraries from Cucumber mosaic virus-infected leaf tissue were produced as part of this study (accessions: SRR6235006-SRR6235021) and all available libraries were evaluated based on their yield ( Figure 4 ). The data produced in this study represented a significant increase (~20%) in the total amount of degradome data available for A. thaliana in the NCBI SRA database.

To select for high confidence degradome supported sRNA–transcript interactions, quality and repeatability was evaluated using a custom scoring system (deg_score). The distribution of deg_scores for cleavage events was evaluated across all degradome libraries. Only the top 1% set of interactions (deg_score > 15) was considered for downstream analysis ( Figure 5 ).

sRNAs were annotated as active if at least one of their predicted targets was confirmed ( Additional file 6: Table S6 ). Experimental support was found for the targeting and cleavage activity of 201 out of 428 of annotated miRNAs, and the number of targets per active miRNA ranged from 1 to 46. In the case of phasiRNAs, 5047 out of 28,464 (~18%) were found to be active, and the number of targets per active phasiRNA ranged from 1 to 39.

Previous studies on smaller datasets had reported the presence of sRNA targets in the 5'UTR (untranslated region), CDS (coding sequence), and 3'UTR regions of the genes (Kamthan et al., 2015). Given the significantly larger number of sRNA target sites identified here, the distribution of sites per region in the genes of target transcripts was evaluated to determine if there is a bias towards specific regions. CDS regions had the highest number of sRNA target sites, followed by 3' UTR and 5' UTRs ( Figure 6 ). It should be noted that these results were not weighted to match the length of these gene features, and the observed distribution could be proportional to the total size of these features in the genome. A very small fraction (~0.03%) of target sites resided in non-coding RNAs.

With the additional degradome data from this study ( Figure 4 ), this is the most comprehensive evaluation of sRNA cleavage activity to date for A. thaliana.

phasiRNAs whose activities were experimentally validated were evaluated based on their sizes and phased registers. Based on their size, a majority of active phasiRNAs were 21 nt long (~71%), with the remainder being 22 nt long phasiRNAs ( Figure 7 ). phasiRNAs derived from non-canonical registers represent a significant majority in both size categories indicating the relevance of these sRNAs ( Figure 7 ).

To integrate datasets and describe the A. thaliana sRNA-mediated regulatory network, a bipartite, directed network was constructed. To differentiate between the theoretical and functional (biologically relevant) components of the network, the network was restricted to the components identified in the miRNA-PHAS loci-phasiRNA annotation and interactions validated by degradome data. sRNA nodes (miRNAs and phasiRNAs) were restricted to those with known precursors and validated targets identified in degradome analysis. Transcript nodes included only those transcripts annotated as precursors (pre-miRNAs, PHAS loci) or validated in the degradome analysis to be direct targets of active sRNAs. The resulting network contained a total of 11,156 nodes, composed of 5,475 sRNA nodes (428 miRNAs and 5047 phasiRNAs) and 5,680 transcript nodes (325 miRNA precursors, 107 PHAS loci, and 5248 target transcripts). These nodes were connected by 13,160 edges; 5,602 of these were involved in the biogenesis of sRNAs and 7,558 edges were involved in sRNA cleavage of transcripts ( Figure 8 , Additional file 7: Table S7 , Additional file 8: Table S8 ).

Three metrics were used to assess the regulatory contribution of the resulting sRNA-mediated network. The proportion, function, and regulatory roles of the genes included in the network were evaluated. There are 33,341 A. thaliana genes defined in the araport11 genome annotation (Cheng et al., 2017). The proportion of annotated genes with degradome data to support interactions with sRNAs was about ~17% (n=5,624) of the total annotated genes. The networks regulatory role was assessed at a functional level using the GO annotations of the genes under sRNA control (i.e., with degradome supported interactions with sRNAs) to determine the biological processes in which sRNAs exert their control. Go slim annotations give a broad overview of the ontology content; using these terms, 22 of 45 functional categories under the biological processes domain were found to be disproportionally enriched in the network (versus a random distribution of categories among annotated genes) ( Figure 9 ). Translation and RNA binding were the only terms out of 45 under the biological processes domain that were found to be underrepresented in the sRNA-mediated regulatory network (corrected p-value= 2.6977E-9 and 5.6680E-12, respectively).

Lastly, to assess the regulatory roles of the genes included in the network, the proportion of genes annotated to be involved in gene expression regulation (GO term: GO: GO:0003700) were evaluated for regulatory interactions with sRNAs. Based on GO annotation from Ensembl genomes release 37 (Kersey et al., 2017), 14% of genes (n=1039) involved in gene expression regulation are under sRNA control, including 21% (n=371) of annotated transcription factors (GO term: GO:0003700). Based on the proportion of total genes (~17%) interacting with sRNAs, the enrichment in the network of most categories under the biological processes domain (22/45), and the proportion of regulatory genes under sRNA control (14%), our results are consistent with the notion that sRNAs play a key regulatory role in plants.

Networks structure analysis allows a multiscale study of complex biological systems such as for the sRNA-mediated regulatory network presented here; global and local features can be identified and compared to related systems. To better understand properties of the resulting sRNA-mediated regulatory network and to determine the interconnectivity between miRNA and phasiRNA, structural features were analyzed. The sRNA-mediated regulatory network consisted of 192 disconnected components, or groups of connected nodes. Given the directed nature of the networks, each one of these corresponded to weakly connected components. The distribution of the number of nodes per component was uneven, with the largest component containing ~96% of the nodes, and the remainder consisting of 15 or less nodes. Network density, a measure of the ratio of the observed number of edges to the maximum number of edges, in this case was very low (<0.001), consistent with results from studies on other biological networks (Leclerc, 2008). By contrast, the clustering coefficient, a measure showing the tendency of a graph to be divided into clusters, was very low (<0.001), whereas in biological systems, higher values are usually observed (Pavlopoulos et al., 2011). In totality, the node degree showed a heavy-tailed distribution, both for indegree and outdegree ( Figure 10A ), describing a limited number of nodes with high degree while the majority have low degree. However, the bipartite nature of the network warrants a separate evaluation of the different types of interactions. Degree distributions for sRNAs and transcripts (lRNAs) were evaluated separately ( Figure 10B ). Negative correlations were found between node degree and abundance, except for the case of sRNA indegree as this indegree distribution is restricted by the nature of sRNA biogenesis and does not allow for testing.