Gossypium arboreum seeds were kindly provided by Prof. Paterson. Plants were sown under control, drought, salinity and heat conditions in green house or growth chamber. For salinity stress plants at 4–5 leaf stage were watered with 200 mM sodium chloride solution for 1 week and leaves were collected. For drought plants at 4–5 leaves growing in pots were not watered for a week and leaves were collected. For heat stress plants were grown in a growth chamber at 50^o for 1 week and leaves were collected. These different stresses were used to identify treatment specific transcripts, which are not expressed in one condition. The purpose here was to identify maximum number of miRNA but the role of miRNAs induced under a particular condition is the subject of another manuscript. Total RNA was extracted and samples were sent to LC Sciences^® United States¹ for sequencing of sRNA using NGS.

Single-end sequencing of sRNA was conducted on Illumina GA-II NGS platform. The Illumina TruSeq sRNA kit as used to prepare the library and size selection of 145–160nt was performed after ligation. Confirmation of RNA sample concentration was performed on a nanodrop spectrophotometer and RNA quality on an Agilent BioAnalyzer.

The reference genome was obtained from the Cotton Genome Project website², assembled by Beijing Genomic Institute (BGI), release1.0 (Li et al., 2014). We used miRBase release 21 (Griffiths-Jones et al., 2006) for plants miRNA mature and stem-loops sequence reference using linux command line. For filtering of non-coding RNAs other than miRNAs and other RNA structural families, we have used Rfam release 12.0 (Nawrocki et al., 2014). The sequence data was deposited into NCBI public databases under following accession; SRA: SRP095309, BioProject: PRJNA357977 and BioSamples: SAMN06163076, SAMN06163077, SAMN06163078, and SAMN06163079.

The sRNA deep sequencing data retrieved from Next Generation Sequencer was preprocessed to remove adaptor sequences and quality filtering using Cutadapt v1.1 (Martin, 2011) and validated using CLC Genomics Workbench v7.5 as well. Both adaptors were removed whereas the some portion of 3′ adaptor “TGGAATTCTCGGGTGCCAAGG” was still present in 36nt raw sRNA reads as the actual payload was 20–24 bp longer (Martin, 2011). The first 6 bases (TGGAAT) were mapped to raw reads to locate their position and the adaptor is then trimmed from the raw reads, e.g., a raw read of TCCATTGTCGTCCAGTGGTTAGGATATGGAATTCTC is trimmed to a read of TCCATTGTCGTCCAGTGGTTAGGATA. Reads were filtered out if the adaptor cannot map to the raw reads or length of the read is less than 15nt after trimming.

We used local blastn v2.2.26 tool from National Center for Biotechnology Information (NCBI) for searching in Rfam v12.0 database using word size of 11. Reads with zero mismatches, zero gaps, 100% identity and 100% query coverage were removed from miRNAome. Similarly, the NCBI blastx tool was used to search in uniprot protein databases using default word size and the same criteria were used to filter out possible messenger RNA (mRNA) residues.

To predict miRNAs from small RNA-Seq data, using expression patterns and following the criteria for plant miRNA annotation, miR-PREFer (Lei and Sun, 2014) was used due to its low false-positive rate and running time. The mappable reads were aligned to the G. arboreum reference genome using bowtie2 (Langmead and Salzberg, 2012) with zero mismatch and end to end read mapping and multiple mappings allowed. The resultant SAM file was fed into miR-PREFer pipeline along with G. arboreum CDS annotation to exclude all predictions within exons. The precursor’s length was set to 150 and reads depth cutoff was set to 5. The minimum mature sequence length range was set to 18 to 24 with no star sequence expression allowed. The resultant candidate precursors were checked against miRNA primary and ancillary annotation criteria (Meyers et al., 2008). Precursor secondary structures was predicted using RNALfold (Lorenz et al., 2011) and minimum free energy index (MFEI) (Zhang et al., 2006b) was selected as less than or equal to -0.85 and necessary downstream analysis was performed using shell and perl scripting under Centos 6.5 Linux environment.

The predicted candidate precursors were clustered into miRNA families using Cd-hit v4.6.1 (Li and Godzik, 2006) using sequence identity cut-off 0.9 and bandwidth of alignment set to 20. To identify conserved miRNAs, we used blastn to match the 8496 mirBase release 21 from 73 plant families with not more than 3 nucleotide difference, i.e., no more than 3 mismatches with zero gaps in mature sequences. The best match was selected if multiple hits are found to satisfy matching criteria. For partial matches, we’ve analyzed the flanking regions matching the previously known sequences and if the matching criterion is satisfied, then that candidate is assigned that particular family.

To elucidate the origin of miRNAs and validate the ancillary criteria of miRNA annotations (Meyers et al., 2008), overlapped miRNA positions were checked using transcripts annotations of the genome. The annotations were converted to UCSC bed file format and their overlaps were predicted using bedtools (Quinlan and Hall, 2010). Hence the intronic and intergenic origins were calculated but miRNAs of repetitive origin are still unknown due to unavailability of repeat annotation. The positional clusters of miRNA genes were found using cluster tool from bedtools using strand specific (-s) and window size (-d) of 10000.

RNA was isolated from different tissues of G. arboreum using total plant RNA purification reagent (Invitrogen Cat No. 12322-012). RNA samples (30 μg) were diluted in 2X loading dye (95% formamide, 0.025% xylene cyanol, 0.025% bromophenol blue and 5mM EDTA) incubated at 65°C for 10 min for denaturation and placing on ice. 300pMol primer (complimentary to the probe) was loaded as positive control. Denatured RNA samples were loaded to 15% PAGE (Urea) and run at 100 V in BIO-RAD Mini protean tetrasystem United States. The gel was stained with ethidium bromide and visualized under UV light once the run was finished. RNA was transferred to the nylon membrane (Amersham Hybond^TM-N++, GE Healthcare) using semi dry blotter (Transblot SD BIO-RAD United States) at 15V for 1 h. RNA was cross linked to the membrane using UV Cross Linker (CL100o ultraviolet cross linker, UVP). The DNA Probe was synthesized using DIG Oligonucleotide 3′-End Labeling Kit (2nd generation from Roche Life Science, Catalog No.3353575910). The blot was prehybridised at 42°C for 3h and hybridized with probe at 42°C for 16 h. The blot was washed with 2X SSC shortly then twice with 2xSSC/0.1% SDS for 15mins and 30mins at 60°C. The blot was developed by chromogenic (NBT/BCIP-T) method using supplied protocol (Amin et al., 2011; Zhang et al., 2013).

We used psRNATarget (Dai and Zhao, 2011) for miRNA target predictions against all G. arboreum annotated coding sequences from BGI, with maximum expectation value of 3.0 and hspsize equal to 20bp, allowed maximum energy to unpair the target site (UPE) equal to 25. Flanking length around target sites for target accessibility analysis was set at 17bp in upstream and 13bp in downstream. Range of central mismatch leading to translational inhibition was set from 9 to 11nt. The targets with more stringent miRNA-mRNA pairing were preferred; if multiple targets were found, all were reported as there may be multiple targets of a single miRNA with different specificity. The target genes were functionally annotated using the KEGG pathway database by the BlastKOALA sequence similarity tool³ and by alignment to the PANTHER v11.1 database of protein HMM models⁴ with the PANTHER scoring tool v2.0.

Total read length obtained from NGS sequencer was 36nt after trimming everything above this size. This included both 5′ and 3′ adaptor sequences, but partial 3′ adaptor was still present due to smaller sRNA payload of 20–24nt, which was further removed. The reads in which 3′ adaptors were missing were most probably related to longer payloads and they were not considered as small RNAs –constituted approximately 58% of raw data from all libraries. This is in concordance with our previous understanding that plants have relatively large and complex small RNA populations within which miRNAs are often a minority (Meyers et al., 2008). We have obtained a total of 22.27 million potential sRNA reads out of which 7.64M, 5.095M, 4.635M, and 4.899M belong to four libraries, respectively. The raw reads were grouped into unique clusters, having a maximum of one ambiguous nucleotide, of 1.6M, 0.87M, 0.409M, and 0.724M reads for the four libraries, respectively. Total read length of 15 to 30 nucleotides was obtained after preprocessing including 3′ sequencing adaptor removal and quality filtering of raw reads, which were further mapped to the genome, known as mappable reads. The read quality on phred scale ranged from 24 to 40 and on average 39nt after this step, which increases the reliability of in silico prediction. The summary of raw data statistics is shown in Table 1. The read length distribution, as shown in Figure 1, suggests the distribution of different types of sRNAs in our libraries. A peak read length of ∼21 and 22 suggests that there is some content of sRNAs other than miRNAs (Meyers et al., 2008) present in our data and 22, 23 and 24 expresses the content of miRNAs, which is higher.

The unique raw mappable read clusters were mapped to uniprot protein databases, to filter all possible degraded protein coding mRNA products. The total reads mapped with perfect homology to uniprot are 1071 (0.014%), 849 (0.017%), 290 (0.006%), 1299 (0.026%) from libraries 1 to 4, respectively. Moreover, the mappable reads were searched in the RNA family (Rfam release 12.0) (Nawrocki et al., 2014) database to filter out other known non-coding RNA (tRNA, rRNA, snRNA, snoRNA, etc.) products and the total raw read count mapped were 460425 (6.02%), 363352 (7.13%), 90441 (1.95%), 279084 (5.7%) for the four libraries, respectively, as shown in Table 1. The uniprot and RFAM filtered data was mapped to the G. arboreum genome sequence (Li et al., 2014) and 1.19M (73.78%), 0.598M (68.72%), 0.298M (72.9%), and 0.452M (62.49%) of the perfectly matched reads were mapped at 4515822, 2332365, 1364987, and 1793749 different genomic locations from the four libraries, respectively. The alignment was obtained in Sequence Alignment Map (SAM) format and fed into miR-PREFer pipeline (Lei and Sun, 2014) for potential precursor prediction. The 250 base flanking regions were excised by considering both 5′ and 3′ mature sequence possibilities. The pipeline, when subjected to all libraries simultaneously, generated 948, 914, 924, and 914 potential precursors from the four libraries, respectively. These potential candidates were short listed based on total reads mapping depth to a precursor and presence in multiple libraries. The precursors with total mature transcript depth greater than 10 in at least one library were considered to be significant, as shown in Supplementary Table S1. The mappable reads were mapped to mature, star and antisense regions and positional and length variants were observed in different transcript mappings. This was because of varying DICER activity or data preprocessing artifacts. The miRNA/miRNA^∗ duplex with positional variation less than 25% were selected for higher confidence (Meyers et al., 2008). Exact duplicates or overlapping predicted precursors were removed from potential precursors. The flowchart of complete data analysis has been shown in Figure 2. The precursor size ranged from 55 to 296nt with the majority within 101 to 150nt (30.8%), 55–80nt (16.7%), 81–100nt (20.7%), and 151–200nt (20.5%) and greater than 200 (11.3%) as shown in Figure 3. The mature sequences have A-U content greater than G-C throughout all mature transcripts except at position 23. Moreover, the total A-U content is greater than 50% in all positions except at position 9, which is immediately after seed region, as shown in Figure 4. The genomic coordinates of the predicted miRNAs are provided in UCSC, BED specification, see Supplementary File (G.arboreum_miRNAs.bed).

The miRNAs that are present in at least one major ancient clade of land plants are referred as ‘conserved’ (Axtell and Bowman, 2008) and those that are more likely to be evolved much more recently in our species of interest or has very limited phylogenetic distribution, are termed ‘novel’ or ‘species specific’ or ‘non-conserved.’ Altogether, 446 candidate precursors were shortlisted based on above-mentioned criteria and clustered into 224 miRNA families. Among the 224 families, 223 were present in more than one library. Candidate mature sequences mapped to 8496 mature miRNAs from 73 plants with no more than 3 nucleotide differences, i.e., sum of gaps and substitutions was taken to be less than 4nt, indicating 48 conserved and 176 G. arboreum specific or novel miRNAs as shown in Table 2 (Voinnet, 2009). In particular, all 21 conserved miRNA families (Axtell and Bowman, 2008) (156, 159, 160, 162, 164, 166–169, 171, 172, 319, 390, 393–399, and 408) were found.

We have found a total of 1274 targets for 206 miRNA families in G. arboreum, which complements the prediction reliability as the ancillary criteria of miRNA annotation (Meyers et al., 2008) and provided in Supplementary Table S2. The rest of 16 miRNAs may be either target-less elements that are evolutionarily transient or their targets are undiscovered due to lack of precision of target prediction tools. All of the conserved miRNAs do have their associated target genes but this is not true for all novel miRNAs. The highest number of target genes per miRNA was found to be 48 and average number of targets for conserved miRNAs (9) is larger than average number of targets for novels (4). The set of 1274 predicted target proteins were functionally characterized in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database⁵ by the BLASTKOALA (Kanehisa et al., 2015) sequence similarity tool. The BLASTKOALA annotated 441 (34.5%) out of 1274 to the KEGG orthologs database and 441 proteins were enriched in 216 different pathways (TS3). The distribution of these Kegg orthologs in different cellular functional categories has been shown in Figure 5. Altogether, 64 genes are involved in genetic information processing, whereas, other processes include environmental information processing (78), general cellular processes (58), disease related pathways (99), organismal systems (58) and metabolism of carbohydrates (32), energy (17), glycans (7), nucleotides (16), amino acids (25), cofactors and vitamins (4), lipids (16) and terpenoids and polyketides (10) and biosynthesis of other secondary metabolites (8). The rest were categorized in other pathways, as shown in TS3. Moreover, the 1274 target genes were annotated using PANTHER^® gene ontology and classification system (Mi et al., 2013). Altogether, 1184 out of 1274 target proteins were annotated to 614 protein families and 130 unique protein classes and 293 unique GO terms. The highest number of targets were related to DNA and nucleic acid binding proteins (150) and transferases (128) followed by RNA binding proteins (107), transporters (101), hydrolases (82), kinases (133), oxidoreductases (66) and enzyme modulators (46), etc. The miRNAs gene network, regulating different functional pathways, along with their target proteins in G. arboreum has been shown in Figure 6. The top 24 PANTHER^® classification having greater than 1% of target genes in terms of protein classes and GO terms are shown in Figures 7, 8, respectively, whereas the rest of the targets are listed in Supplementary Table S4.

Varying miRNA copy numbers may affect target regulation capability of plants. These copy number variants (CNVs) precursors may or may not have exactly same nucleotide sequences, but their mature segments are more conserved. The gain or loss of miRNA genes is due to tandem and segmental duplication events (Li and Mao, 2007). Plants have relatively higher numbers of miRNA gene family members as compared to animals but smaller total numbers of miRNA genes. The average number of precursors in 224 genes was 2.0, which is in concordance with previous reports for plants (Li and Mao, 2007). The maximum copy number was observed for gar-miR7504 with 17 members and gar-miR166, gar-miR8771, gar-miR7484, gar-miR171, gar-miR156, gar-miR396, gar-miR7502, gar-miR164, gar-miR395, gar-miR8672 and gar-miR167 have 14, 14, 11, 9, 8, 7, 7, 6, 6, 6 and 5 precursor members, respectively. The other conserved members have precursors greater than or equal to 2 are shown in Figure 9. The relative expression of miRNA genes is associated to their copy numbers or methylation or both (Aure et al., 2013); hence individually, a strong correlation with copy numbers cannot be predicted in complete miRNome. We have used SPSS for determining the bivariate Pearson correlation coefficient (PCC) between samples and the miRNA copy number with sig. (2-tailed) and bootstrap equal to 1000 at 95% confidence interval. We have observed significant but not highest positive correlation, i.e., PCC equal to 0.638 (p-value: 0.001) for highly conserved miRNA gene families (Axtell and Bowman, 2008) and for all 48 conserved miRNAs, the PCC is equal to 0.382 (p-value: 0.007). The correlation of all 176 predicted novel miRNAs is not significant with PCC equal to -0.52 (p-value: 0.494), suggesting their transient nature.

Altogether, 224 miRNA families originate from 446 different genomic origins due to varying number of copies and isoforms. All the chromosomes have nearly symmetrical number of miRNAs distributed over them, except chromosome 5; with only 17 miRNA genes were found. Chromosome 10 has the greatest number of miRNAs (49) predicted, followed by chr8 (48), chr12 (41), chr1 (40), chr4 (36), chr7 (36), chr11 (33), chr3 (32), chr9 (31), chr6 (30), chr13 (28), chr2 (23), and chr5 (17). Altogether, 19 miRNAs originated from 32 different intronic origins of different annotated protein coding genes including their isoforms and copy numbers, listed in Supplementary Table S5. These include conserved miRNA families gar-miR172, gar-miR395, gar-miR7484, gar-miR7502, gar-miR7504, gar-miR8634, gar-miR8674, gar-miR8779 and 11 novel miRNAs. Six out of seven copies of gar-miR7502 were originated from introns of different protein coding genes. Plants have less common phenomenon of clustered miRNA genes than animals (Zhang et al., 2006a; Voinnet, 2009). These clusters are mostly homologous family members of conserved miRNAs (Axtell and Bowman, 2008) and within a few kilobases in distance (Li and Mao, 2007). We have found 9 such clusters with a window size up to 8kb, including gar-miR171 and gar-miR7504 on chr1 and chr8, respectively, having intergenic distance of 2kb. The gar-miR2949 and gar-nov47 made clusters within 4kb distance on chr4 and chr8, respectively. The other clusters were made by gar-miR8677 and some species specific miRNAs within 8kb intergenic distances. Increasing the window size yields a few more clusters but they are less likely to be polycistronic. All of the clusters have two precursor members except miR8677, which contains three such genes in 8kb vicinity, as shown in Supplementary Table S6. The genomic origins of all of the predicted genes are shown in Figure 10, elaborating the distribution of miRNA genes over 13 chromosomes.

In order to validate the effectiveness of our bioinformatics pipeline, we chose two G. arboreum specific miRNAs (gar-novel-165: UUUGUACUUUAGAUGUCUCUC and gar-novel-167: AUCGAUAUUGUUCUCUUUGGG). These miRNAs originate from single genomic locus where gar-novel-165 is transcribed from positive strand and gar-novel-167 from negative strand. Their secondary structure is having MFEI value less than -1 whereas target of one of them (gar-novel-165) is not found. Both of these pre-miRNAs have good base pairing in the stem region MFE -76.7 and -91, respectively. The above selection criteria make these miRNAs ideal for validation and hence extrapolating the results for our in silico mapping based pipeline. We have extracted RNA from leaf, shoot and root tissues and leaves show relatively more expression level. Moreover, wheat leaves RNA was used as negative control and primers (complimentary sequence) as positive controls, as shown in Figure 11. This analysis showed that these miRNAs specifically expressed in G. arboreum and have different expression levels in different tissues.