Computational prediction of microRNAs in marine bacteria of the genus Thalassospira

MicroRNAs (miRNAs) are key players in regulation of gene expression at post-transcription level in eukaryotic cells. MiRNAs have been intensively studied in plants, animals and viruses. The investigations of bacterial miRNAs have gained less attention, except for the recent studies on miRNAs derived from Streptococcus mutans ATCC 25175 and Escherichia coli DH10B. In this study, high-throughput sequencing approach was employed to investigate the miRNA population in bacteria of the genus Thalassospira using both the miRDeep2 and CID-miRNA methods. A total of 984 putative miRNAs were identified in 9 species of the genus Thalassospira using both miRDeep and CID-miRNA analyses. Fifty seven conserved putative miRNAs were found in different species of the genus Thalassospira, and up to 6 miRNAs were found to be present at different locations in the T. alkalitolerans JCM 18968T, T. lucentensis QMT2T and T. xianhensis P-4T. None of the putative miRNAs was found to share sequence to the reported miRNAs in E. coli DH10B and S. mutans ATCC 25175. The findings provide a comprehensive list of computationally identified miRNAs in 9 bacterial species of the genus Thalassospira and addressed the existing knowledge gap on the presence of miRNAs in the Thalassospira genomes.


Introduction
MicroRNA (miRNA) is a class of small, non-coding RNA molecules containing 19-22 nucleotides. Since the discovery of miRNA in Caenorhabditis elegans [1], a large number of predicted miRNA molecules has been reported in animals, plants and viruses as key players in regulation of gene expression networks [2][3][4][5]. In bacteria, the small non-coding RNAs (sRNAs) have been demonstrated to have a similar function to eukaryotic miRNA, in modulating the target mRNAs in various ways at a post-transcriptional level [6]. A number of sRNAs have been identified in bacteria, some of which were identified in marine bacteria such as Vibrio spp. and Synechococcus spp. that are functional analogues to plant miRNAs in response to environmental changes [7,8]. The investigations of bacterial miRNAs, however, have gained little a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 attention, except for the recent studies on miRNAs derived from Streptococcus mutans ATCC 25175 and Escherichia coli DH10B [9,10].
The genus Thalassospira includes Gram-negative, aerobic and halophilic bacteria dwelling in a marine environment [11]. Bacteria of this genus are involved in the biodegradation of a variety of hydrocarbons [12][13][14]. For example, T. tepidiphila 1-1B T and T. povalilytica Zumi 95 T have the ability to degrade polycyclic aromatic hydrocarbons and polyvinyl alcohol [12,13]; and T. australica NP 3b2 T is able to utilise poly (ethylene) terephthalate (PET) plastic as a carbon source [14]. Currently the genus Thalassospira is comprised of 10 validly named species, of which the whole genome sequences of four species (T. australica NP 3b2 T , T. lucentensis QMT2 T , T. profundimaris WP0211 T and T. xiamenensis M-5 T ) have been assembled and deposited in public databases [15][16][17]. The presence of miRNAs in bacteria of this genus has, however, not been demonstrated. Here, the investigation of miRNA populations was carried out using 9 validly named Thalassospira species, i.e., T. alkalitolerans JCM 18968 T , T. lucentensis QMT2 T , T. mesophila JCM 18969 T , T. povalilytica Zumi 95 T , T. profundimaris WP0211 T , T. tepidiphila 1-1B T , T. xiamenensis M-5 T , T. xianhensis P-4 T and T. australica NP 3b2 T . Since Thalassospira lohafexi 139Z-12 T was not validly published until 2015 [11], this species was not included in this study. A deep sequencing approach has been used to identify miRNAs in bacteria in other studies [10]. Thus, the aim of this work was to identify the potential miRNAs in bacteria of the genus Thalassospira using computational approaches from small RNA sequence dataset generated by high-throughput sequencing technology.

Bacterial strains and growth conditions
Nine type strains of Thalassospira species were obtained from various culture collections and used in this study.  [18,19]. Bacterial cells were collected by centrifugation and used for RNA extraction or stored at -80˚C in marine broth 2216 (BD, USA) supplemented with 20% (v/v) glycerol.

RNA isolation and construction of small RNA libraries
Total RNAs were extracted from bacterial cells of Thalassospira species using TRIsure reagent (Bioline, Australia) according to the manufacturer's instructions. The RNA concentration and purity were determined by measuring the absorbance at 260 nm and A260/A280 ratio using Nanodrop (ThermoFisher, Australia). The integrity of extracted RNAs was examined by a Bioanalyzer (Agilent Technology, USA). A RIN value � 5 was applied for selection of good RNA quality with RIN (RNA integrity number) generated for each sample based on the ratio of ribosomal bands and the presence or absence of degradation products on the electrophoretic image. Small RNA-Seq libraries were prepared using the TruSeq Small RNA library preparation kit for Illumina (New England Biolabs, USA). Briefly, rRNA was depleted with Ribo Zero. RNAs were fragmented by heat and divalent cations and the 1 st strand cDNA was synthesised with SuperScript II Reverse Transcriptase (Invitrogen). The DNA fragments were 3'adenylated, ligated with 3' and 5' adapters and amplified via PCR (13 cycles). The amplified products were loaded on 6% polyacrylamide for small RNA (18-35 nucleotides) selection. The amplicons of small RNA were purified and sequenced using Illumina Hiseq2500 sequencing platform with 50 bp single end chemistry. The raw sequencing data generated from this study was submitted to the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA505357.

Bioinformatics analysis of bacterial miRNAs
The raw reads (18-35 nucleotides) generated from the Illumina sequencing were first filtered by removing any contaminations including rDNA, Illumina small RNA adapter sequences or low-quality reads with < 18 nt length using Cutadapt [20]. The remaining reads after filtering were then aligned against the Thalassospira genome sequences available in the NCBI database using Bowtie (version 1.1.2) with parameters setting (-n 0 -e 80 -l 18 -a -m 5-best-strata) [21]. The whole genome sequences used as references for these 9 Thalassospira species were listed in S1 Table. The mapped sequences to genomes were further processed to remove redundancy and normalised expression was computed as CPM (Reads counts per Million) [22]. After normalisation, these mapped sequences were used to predict miRNAs. Two different miRNA prediction algorithms, miRDeep2 and CID-miRNA were used to detect miRNA in this study. The miRDeep2 software (version 2.0.0.7, http://www.mdc-berlin.de/rajewsky/ miRDeep) was utilised to predict miRNAs based on an investigation of the secondary structure of the miRNA precursor sequences and integration of miRNA precursors with Dicer, providing the precursors and mature miRNA sequences [23]. CID-miRNA [24] is a method to identify miRNA precursors based on secondary structure filter and an algorithm of stochastic context free grammar (SCFG) [24]. This method only predicts the miRNA precursors. Therefore, MatureBayes (http://mirna.imbb.forth.gr/MatureBayes.html) was employed to identify mature miRNAs based on the sequence and structure of the miRNA precursors provided by CID-miRNA [25]. The small sequences mapped to the genome were processed through miR-Deep2 ver2.0.0.7 as well as CID-miRNA and MatureBayes for identifying mature miRNAs with default parameters. The detected mature miRNAs were further investigated by searching any shared miRNA sequences generated in both methods using BLAST algorithm (version 4) with default setting. The potential precursor sequences predicted by both methods were also BLAST-searched to find any shared precursor sequences. Any shared miRNA sequences were then confirmed using ClustalW2 (https://www.ebi.ac.uk/Tools/msa/clustalw2/) with default setting. In order to find any conserved miRNA across the genus, the predicted mature miR-NAs were aligned to those reported in Escherichia coli DH10B and Streptococcus mutans ATCC 25175 [9,10]. The candidates were also aligned to miRNAs predicted within the genus Thalassospira to identify any conserved miRNAs across all bacteria within the genus.

Phylogenetic analysis and genome comparison
The 16S rRNA gene sequences of validly described Thalassospira species were collected from NCBI (https://www.ncbi.nlm.nih.gov/) and compared using the ClustalW2 (https://www.ebi. ac.uk/Tools/msa/clustalw2/) with default setting. Evolutionary phylogenetic trees were then generated using the neighbour-joining (NJ) algorithm [26]. Genetic distances for the NJ tree were calculated using Kimura's two-parameter model [27], with the robustness of 1,000 replications, using MEGA 7 software [28].
T. profundimaris WP0211 T , T. tepidiphila 1-1B T , T. xiamenensis M-5 T , T. xianhensis P-4 T and T. australica NP 3b2 T . Next-generation small RNA-seq was used for detection of the presence of small RNAs. The results yielded 242,204,737 sequence reads of 18-33 nucleotides (nt) length, ranging from 22.99 to 31.43 million reads for 9 libraries ( Table 1). The reads were curated for any contamination or artefacts such as rDNA or Illumina small RNA adaptor sequences using Cutadapt [20]. The cleaned reads were then aligned to the Thalassospira genome sequences available in the NCBI database using Bowtie software [21].
Since the whole genome sequences of T. australica NP 3b2 T , T. lucentensis QMT2 T , T. profundimaris WP0211 T and T. xiamenensis M-5 T were available in the NCBI database, the small RNA sequence reads were aligned to these genome sequences. The remaining 5 Thalassospira species that did not have the whole genome sequence were aligned to these references based on their phylogenetically closest species [14] (S1 Table, S1 Fig). By aligning 34.23 million collapsed reads from nine species of the genus Thalassospira to the genome sequences available in NCBI, it was found that a proportion of cleaned reads was mappable to the reference genomes, ranging from 38.42% to 71.94% for 4 species with available genome sequences and from 2.55% to 28.16% for the remaining 5 species without whole genome sequenced ( Table 1). The aligned reads were further analysed for miRNA identification.

Identification of miRNAs using the mirDeep2 method
Identification of bacterial miRNAs is in early stages in comparison to that in animals and plants. Therefore, no specific method has been developed for investigation of miRNAs in bacteria. miRDeep2 software (http://www.mdc-berlin.de/rajewsky/miRDeep) was first presented as an algorithm to predict miRNA in human and animals [23]. This software was later used to successfully identify miRNAs in plant [29]. Based on the applicability in different organisms, miRDeep2 software was, therefore, employed to detect putative bacterial miRNAs from the millions of short sequences generated from the dataset. The miRDeep2 software was designed to detect putative miRNAs based on an investigation of the secondary structure of the miRNA precursor sequences. The potential miRNA precursors were then integrated with a model for miRNA precursor processing using Dicer, releasing the mature miRNA sequences, star sequences (also called miRNA5p and miRNA3p sequences) and the loop [23]. By using miR-Deep2, 86 potential miRNA precursor sequences could be identified from the dataset. The mature miRNA, star sequence and the loop were also identified within each precursor (Fig 1).
The identified 86 putative mature miRNAs were located in either the 5' arm or the 3' arm of the precursors, with 11 sequences from T. australica NP 3b2 T , 11 from T. lucentensis QMT2 T , 14 from T. profundimaris WP0211 T , 6 from T. xiamenensis M-5 T , 11 from T. alkalitolerans JCM 18968 T , 5 from T. mesophila JCM 18969 T , 4 from T. povalilytica Zumi 95 T , 9 from T. tepidiphila 1-1B T and 15 from T. xianhensis P-4 T (Table 2) being obtained.  The identified putative miRNAs are 18-25 nucleotides in length (Table 2), which is in the range previously reported for animal and plant miRNAs [30]. These results are also consistent with the length distribution of miRNAs recently detected in E. coli DH10B and Streptococcus mutans ATCC 25175 [8,9]. Based on the number of read counts, these putative miRNAs appear to have various expression levels among nine libraries, ranging from one to thousands of reads in each library, under the normal growth condition. Among the 86 detected putative miRNAs, the highest expression was found in T.prof_5p_16873, with 50145 reads, while the lowest expression was shared between T.alka_5p_4136 and T. alka_5p_1889 with only one read each being found in the libraries. Of the 86 detected putative miRNAs, thirteen were in high abundance, with over a thousand reads [10] with three miRNAs from T. australica NP 3b2 T , four from T. lucentensis QMT2 T , four from T. profundimaris WP0211 T and one each from T. tepidiphila 1-1B T and T. xiamenensis M-5 T being identified. The high degree of expression of these putative miRNAs suggested that they may play specific roles involving in the growth and development processes of these bacteria (Fig 2).

Identification of miRNAs using CID-miRNA method
The 86 putative miRNAs detected in nine Thalassospira species using miRDeep method appeared to be rather low in number compared to 400 miRNAs identified in E. coli DH10B [10] and 900 miRNAs in Streptococcus mutans ATCC 25175 [9]. Therefore, an alternative method, CID-miRNA [24] was employed in the expectation of identifying further potential miRNAs from the dataset. CID-miRNA is a web-server developed for the identification of the miRNA precursors based on the secondary structure filter and an algorithm of stochastic context free grammar (SCFG) [24]. The server was first used to predict the potential miRNA precursors in human genome [24] which was later applied in animals such as mouse, zebrafish and sea squirt [31]. Using this method, 449 potential miRNA precursors were identified from over 242 million reads of T. australica NP 3b2 T  In order to identify the mature miRNAs, these precursor sequences were further analysed by employing another web tool, MatureBayes (http://mirna.imbb.forth.gr/MatureBayes.html).
MatureBayes is designed to incorporate to the Naïve Bayes classifier for identifying the putative mature miRNA molecules based on the sequence and structure of the miRNA precursors [25].  Table). These identified putative miRNAs were used for further analysis of miRNA conservation.

Comparison of bacterial miRNAs identified using miRDeep2 and CID-miRNA
In order to determine whether any miRNAs were detected by both methods, miRNA data generated from miRDeep were BLAST-searched to those produced by CID-miRNA. The results showed that there were no shared putative mature miRNA sequences. miRDeep2 predicts putative miRNAs based on distribution of minimum free energy (MFE) and stability of secondary structure established for nematode, C. elegans, while CID-miRNA utilised secondary structure filter and stochastic context free grammar trained on human miRNAs [23,24]. The differences of parameter setting may cause no overlap of putative mature miRNA sequences detected in both methods. Thakur et al [32] also pointed out that application of new parameters improved the accuracy of plant miRNA prediction using miRDeep in compared to the default setting for animals.
Study on miRNA biogenesis showed that one or more mature miRNAs can be produced from one pre-miRNA molecule [33]. These identified putative miRNAs, therefore, were aligned to the potential precursor miRNA sequences generated by both methods. BLASTsearch between putative mature miRNA sequences identified from miRDeep2 and the potential precursor miRNAs detected by CID-miRNA was carried out and vice versa. 5 putative mature miRNA sequences identified by CID-miRNA (T.luce_5p_228121, T.luce_3p_228121, T.xian_5p_2738, T.xian_5p_2740 and T.xian_3p_2740) were found to locate in 3 potential precursor miRNAs producing T.luce_5p_11956, T.xian_5p_4477 and T.xian_3p_4479 from miRDeep2. The precursor miRNA sequences producing 5 putative mature miRNAs from CID-miRNA were then aligned to 3 potential precursor miRNAs from miRDeep2 using ClustalW2. The result showed the overlap of 3 precursor miRNAs predicted in both methods in which the putative mature miRNAs identified by each method located at different positions within the precursor sequences (Fig 3).
The shared sequences from the same Thalassospira species proved the probability of putative miRNAs finding in both methods. Thus, a total of 984 putative miRNA candidates were preliminarily obtained from both methods with default setting together for 9 species (Table 3). Further work, such as real time PCR or Northern blot analysis, is needed to verify these putative miRNAs.

Conservation of miRNAs in the genus Thalassospira and miRNAs previously described in Streptococcus mutans ATCC 25175 and E. coli DH10B
Previous studies showed the conservation of some miRNAs across animals or the plant kingdom [34,35]; however, novel sequences can only be found in a particular species. The conserved and novel miRNAs in bacteria still remain largely unknown, in comparison to intensive studies of miRNAs in eukaryotic organisms and viruses. In this study, 984 putative miRNA candidates of the genus Thalassospira were BLAST-searched against those reported for E. coli DH10B (400 miRNAs) and S. mutans ATCC 25175 (900 miRNAs) [9,10] in order to identify any conserved miRNA. However, no conserved sequences could be found without any mismatch or with three and fewer nucleotide substitution. E. coli is an enteric bacterium commonly found in low intestine of humans and animals [36], S. mutans is an oral pathogen that causes human dental caries [37], while bacteria of the Thalassospira genus are environmental bacteria [14]. The differences in their characteristics may influence the low degree of conservation of miRNAs among these bacteria. The results obtained in this study is in agreement with previous studies that reported the lack of significant sequence similarity in non-coding RNA homologues of different bacterial species [38,39].
In order to identify any conserved miRNAs among species of the genus Thalassospira, the putative miRNA sequences obtained were compared. It appears that the nine species shared 57 common putative miRNAs. Among these, T. profundimaris WP0211 T and T. tepidiphila 1-1B T had the highest number of conserved miRNAs, with 45 sequences presenting in both species. Five putative miRNAs in T. alkalitolerans JCM 18968 T were also found in the T. lucentensis QMT2 T , while T. xianhensis P-4 T and T. xiamenensis M-5 T shared 4 conserved sequences. One sequence was shared between T. australica NP 3b2 T and T. tepidiphila 1-1B T , T. lucentensis QMT2 T and T. tepidiphila 1-1B T and T. lucentensis QMT2 T and T. mesophila JCM 18969 T ( Table 4). As seen from the data, T. lucentensis QMT2 T has identical miRNA sequences as T. alkalitolerans JCM 18968 T (5 sequences), T. tepidiphila 1-1B T (1) and T. mesophila JCM 18969 T (1), while T. tepidiphila 1-1B T also shared sequences with T. profundimaris WP0211 T (45) and T. australica NP 3b2 T (1).
A comparative analysis of 16S rRNA sequence similarities revealed that the highest number of common miRNAs shared by the phylogenetically close species; e.g., T. tepidiphila 1-1B T were found to be phylogenetically closely related to T. profundimaris WP0211 T (99.3% sequence similarity); T. xianhensis P-4 T and T. xiamenensis M-5 T also shared 99.3% 16S rRNA similarity, while T. mesophila JCM 18969 T and T. alkalitolerans JCM 18968 T shared 95.3% and 94.9% of 16S rRNA sequence similarity, respectively, with T. lucentensis QMT2 T [14]. Conserved miRNA sequences found in these species may indicate a close genetic relationship and that these miRNAs have a similar role in the regulation of the growth and development of bacteria. It is interesting to note that two species, T. australica NP 3b2 T and T. tepidiphila 1-1B T , which are capable of hydrocarbon degradation [12,14], shared one conserved miRNAs and the potential targets of this putative miRNA needs to be investigated for any role related to this process. In addition, T. alkalitolerans JCM 18968 T , T. lucentensis QMT2 T and T. xianhensis P-4 T were also found to have the same sequence presented in different genomic locations ( Table 5). These miRNAs can have an influence on their expression and function at different genomic locations [40]. It will be of great interest to identify the target mRNA of these miR-NAs and investigate their roles and the mechanisms of gene regulation in the physiology of these unique environmental bacteria.

Conclusions
Over 242 million reads of 18 to 33 nucleotides length were generated from nine bacterial species using high-throughput sequencing technology. Using miRDeep2 and CID-miRNA analyses, a total of 984 putative miRNAs were eventually identified, with typical miRNA length of 19-25 nucleotides. Compared to other species, these detected putative miRNAs were not conserved to those reported in E. coli DH10B and S. mutans ATCC 25175. This study presents the first comprehensive list of computationally identified miRNAs in 9 bacterial species of the genus Thalassospira without experimental verification. The further work, however, is needed to validate these candidates experimentally. In addition, further identification of miRNA targets will provide insights into the fundamental functions of miRNAs in the physiology of these bacteria.