The IGRs of M. huakuii 7653R with a length ≥50 nt are 5125, and these were compared with four other Mesorhizobial strains by WU-BLASTN to identify conserved candidate sequences in the IGRs. The resulting homologies of IGRs revealed by these genome-wide comparisons were further analyzed. Raw data from ~1,500 sequences were detected after WU-BLAST analysis. Each of the conserved intergenic regions in these candidate sequences was scanned individually using eQRNA and RNAz, and overlapping sequences were determined. Meanwhile, the promoters and Rho-independent terminators were predicted using Promoter 2.0, RNAMotif and Erpin, and the resulting information was used to further assess the sequences mentioned above (Tables S1, S2). Finally, 40 conserved sequences were identified as potential sRNAs (Table 1).

To find the annotated known sRNAs in 7653R, we performed Rfam searches using the above potential sRNAs. Of the 40 predicted sRNAs, 24 were annotated in the Rfam database, thus, the other 16 sRNAs are regarded as putative novel sRNAs (Table 1). Several known sRNAs, which are highly conserved in bacteria, were successfully predicted, including 6S RNA (MH_32), 4.5S (MH_s40), RNase P class A (MH_s20), and tmRNA (MH_s19), suggesting that the methods used in this study worked well. A few glycine, cobalamin and TPP riboswitches were also detected. Notably, most of the sRNAs (35 of 40) were located on the chromosome, with only a few (three in pMhu7653Ra and two in pMhu7653Rb) located on megaplasmids, indicating possible differential conservation of the IGR sequences between the chromosome and megaplasmids among the tested bacteria. Predicting the sRNAs that locate on the symbiotic megaplasmids seems to be difficult.

As noted above, the genomic comparison identified 40 potential sRNAs. Further analyses were conducted to describe their transcription units and transcriptional start regions using Illumina high-throughput sequencing data. These analyses were performed on RNA isolated from 28- and 50-dpi (days post-infection) symbiotic nodules and free-living cells of 7653R (see Materials and Methods) that were treated with TAP. The depths of RNA-seq were shown in Tables S3, S4.

A comparison of the analysis between the predicted sRNAs and RNA-seq products showed that, apart from 5 riboswitches, 19 functional annotated sRNAs and 15 of 16 putative novel sRNAs were detected in the sequencing dataset, and only one was not expressed in plant nodules or in free-living cultured cells (Table S5). Of the riboswitches (cis-regulators of mRNA that control the translation of mRNA), MH_s14, was observed to be highly expressed in all sequenced samples. In total, other than the 5 riboswitches, 34 of the 35 predicted sRNAs had corresponding RNA-seq sequencing products, and only MH_s6 was absent in the sequenced samples. These results suggest that MH_s6 expression might be induced by different stress conditions rather than the logarithmic growth state or symbiosis or that the MH_s6 transcript levels were below the threshold of detection using RNA-seq. Notably, except for a few sRNAs (MH_s1, MH_8, MH_s10, and MH_s20), the lengths of the sequenced sRNAs in the RNA-seq data were slightly longer (10 of 34 sRNAs) or shorter (20 of 34 sRNAs) than the predicted lengths. The predicted transcriptional start and stop regions of the sRNAs were also compared with the RNA-seq data, and the results showed that 12 and 5 predicted sRNAs had approximate matches to their sequenced start and stop positions, respectively. Thus, more than 50% of the predicted sRNAs showed transcript length discrepancies with their RNA-seq products. This result has also been observed by others (Gonzalez et al., 2008; Valverde et al., 2008; Otto et al., 2012; Mentz et al., 2013), and the main reason may be ascribed to RNA processing or degradation in the cell (Madhugiri et al., 2012; Saramago et al., 2014).

There are four sRNAs, including 6S RNA, 4.5S RNA, bacterial RNase P class A, and tmRNA, that are highly conserved and regarded as housekeeping sRNAs in bacteria. In this study, these four housekeeping sRNAs were detected by both bioinformatic analysis and RNA-seq, and their expression profiles and expression levels are respectively shown in Figure 1 and Table 2. In agreement with the earlier reports, RNase P (MH_S20), an omnipresent endoribonuclease, was highly expressed under both free-living and symbiotic conditions and displayed the highest expression levels compared with other housekeeping sRNAs (Regalia et al., 2002; Vercruysse et al., 2010). The 4.5S RNA (MH_S40), which forms the signal recognition particle (SRP) with the Ffh protein, was also relatively highly expressed, especially in symbiotic bacteroids (~1.5-fold higher than in free-living cells). The 6S sRNA (MH_S32) has been reported to be more abundant during the stationary phase of cell growth. Its function is to block δ⁷⁰-dependent transcription in alpha-Proteobacteria. In our sequencing results, 6S was strongly expressed in exponential-phase cells (OD₆₀₀ = 0.8), and abundant transcripts were also detected in mature (28 dpi) and senescent (50 dpi) nodules. The expression of tmRNA (MH_S19), which rescues stalled ribosomes and tags incomplete polypeptides for degradation, was lower in both free-living and symbiotic conditions, and down-regulated in mature and senescent nodules compared with exponential-phase cells (>2-fold).

ctRNA (counter-transcribed RNA) is a highly conserved small RNA located in repB–repC intergenic regions. It regulates the replication and incompatibility of repABC-type plasmids (MacLellan et al., 2005; Vercruysse et al., 2010; Rivera-Urbalejo et al., 2015). M. huakuii 7653R possesses two symbiotic megaplasmids required for symbiosis and nitrogen fixation. In our results, the levels of ctRNA (MH_s1 in pMhu7653Ra and MH_s4 in pMhu7653Rb) were very high in free-living cells, approximately 3- to 4-fold higher in 28 dpi and 50 dpi nodules, respectively (Figure 2). Moreover, a cis-encoded highly abundant sequencing product was observed in the repB–repC intergenic regions of pMhu7653Ra (located at 2601–2727) and pMhu7653Rb (located at 2344–2474), respectively, so we speculated that this cis-encoded product could be formed by ctRNA terminating the transcription of repC. The results further showed that ctRNA is trans-transcribed and highly expressed in rhizobia (Yip et al., 2015).

To confirm the existence of the predicted novel sRNAs and to understand their basic functions, 9 of the 15 potential novel sRNAs were selected randomly and detected by conducting Northern blotting on the total RNA extracted from cultured cells grown under non-stressed and stressed conditions (Figure 3), and 5S RNA was used as a positive control (Figure S1). As shown, the sizes of the nine sRNAs that were tested roughly consistent with those of the RNA sequencing dataset. Two candidate sRNAs (MH_s10 and MH_s22) exhibited two bands, and MH_s25 showed a complex banding pattern. Six of the nine tested sRNAs (MH_3, MH_s10, MH_s15, MH_s22, MH_s25, and MH_s36) displayed strong signals in response to a variety of external stimuli. MH_s7, MH_s11 and MH_s39 were expressed under only rare growth conditions. In addition, after analyzing the sequencing products again, another highly expressed sequence (MH_s25) was detected that located between chromosomes 4531972 and 4532239. Because this transcript product lacked typical conserved structures and 4532144–4532330 is close to the predicted region (4532142–4532343), we regard MH_s25 (4532144–4532330) as a potential sRNAs. Apart from 4532144–4532330 and 4532142–4532343, the further information on other Northern blotting signals of MH_s25 could not be detected in the RNA-seq data. The main reason may be that the bio-informational analysis and compilation were performed to deal with an enormous number of non-coding RNAs in the RNA-seq data, which led to some shorter or less abundant transcripts usually be merged into the longer ones in a same transcription unit, so the sequence information for complex banding pattern appearing in Northern blots was difficult to obtain.

The expression and abundance of sRNAs were affected by various external stimuli, and these condition-specific expression patterns help to infer the functions of different sRNAs (Vercruysse et al., 2010). In our study, most of the tested sRNAs (MH_3, MH_s10, MH_s15, MH_s22, MH_s25, and MH_s36) were highly abundant in symbiotic, micro-oxygen, oxidative and salt-stress conditions, but few sRNAs accumulated under heat, cold, acid, or alkali stress conditions. MH_s15 and MH_s25 were widely expressed sRNAs and accumulated heavily under most of the tested induction conditions, suggesting that they function in multiple cellar processes. The expression patterns of MH_s7, MH_s11 and MH_s39 were unique, and these sRNAs accumulated under only a few of the tested conditions. Nearly all the detected sRNAs (8 of 9, MH_s11 being the exception) were expressed under symbiotic conditions or symbiosis-mimicking (micro-oxygen and H₂O₂) stimuli, indicating that a large number of sRNAs are required during rhizobial symbiotic nitrogen fixation. However, more studies are required to characterize the functions and mechanisms of action of these sRNAs.

To further examine the characteristics of the detected sRNAs, a BLASTN search (NCBI) was performed to analyze the conservation of all the detected sRNAs in M. huakuii 7653R. The results showed that most of the detected sRNAs are highly conserved and only distributed in Mesorhizobium spp. and Bradyrhizobium japonicum, except for MH_s11, which has homologs in Bartonella spp. The detected sRNAs contain highly conserved secondary structures and more than two loops (Figure 4). MH_s15 has highly conserved structures, with a single, clear band shown on the Northern blots, and abundant polyA/U in the loop and tail of its secondary structure, which is the putative Hfq-binding site. In addition, the transcription levels of the genes upstream and downstream of MH_s15 were upregulated in 28 and 50 dpi nodules, respectively, compared with free-living cells (from 1.5- to 7-fold), according to the RNA-seq results. Therefore, the target genes of MH_s15 were analyzed by using the IntaRNA software and the potential conserved interactive regions were shown in Table S6, but further experiments are needed to characterize the interactions between MH_s15 and its target genes and determine the influence of MH_s15 on the symbiosis between M. huakuii 7653R and host.

In order to better understand the potential function of the 16 new sRNA candidates identified in this work, the target genes of 16 putative novel sRNAs were predicted by using the webserver IntaRNA. These candidate target genes were subjective to GO annotation and enrichment analysis. The WEGO output of the target genes was shown in Figures S2, S3. The enriched GO terms and functional categorization of candidate small RNA target genes were summarized and demonstrated in a pie chart (Figure 5, Table S8). As it showed in Figure 5A, classified in terms of cellular component, a majority of the predicted target genes is enriched in the function of membrane systems (40%), cell part (34%) and intracellular part (15%). When classified as molecular functions, a total of 38% of the enriched GO terms are related to catalytic activity, and 21% of the target genes are involved in the binding of nucleic acid, nucleoside and nucleotide (Figure 5B). Regarding the enriched GO term of biological process, as shown in Figure 5C, the target genes were distributed widely in various biological functions, from cellular process, to biological regulation, to primary metabolic process, to nitrogen compound metabolic process, and to response to stimulus and stress. Importantly to point out, 48% of the target genes are enriched in all kinds of metabolic processes, such as primary metabolism, second metabolism, nitrogen compound metabolism, alcohol metabolism etc. Remarkably, 9% of the target genes are involved in the nitrogen compound metabolic process, indicating the important roles of the sRNAs playing in the nitrogen fixing and utilization. Moreover, 16 and 17% of the target genes are respectively involved in biological regulation and cellular process, and 2% of the target genes are connected with response to stimulus and stress (Figure 5C).

The comparative analysis of IGRs was carried out among M. austrqalicum WSM2073, M. ciceri WSM1271, M. Loti MAFF303099, M. opportunistum WSM2075 and M. huakuii 7653R. Their genomic sequences and protein annotations were downloaded from the NCBI bacterial databases (ftp://ftp.ncbi.nih.gov/genomes/Bacteria). The predictions of sRNAs were performed as previously described (del Val et al., 2007). IGRs with a length ≥50 nt were scanned and homologous regions were clustered using WU-BLASTN (http://blast.wustl.edu) in the tested genomes with E-values ≤ 1e⁻⁵. The alignments were interrogated using eQRNA version 2.0.4 (ftp://ftp.genetics.wustl.edu/pub/eddy/software/qrna.tar.Z) as described previously (del Val et al., 2007). RNAz version 0.1.1 (http://www.tbi.univie.ac.at/~wash/RNAz) was used to screen the conserved intergenic regions in these candidate genomic sequences, and overlapping sequences between eQRNA and RNAz were selected. Subsequently, pairwise alignments in both eQRNA and RNAz were calculated using RNAfold (ViennaRNA-2.1.1) to determine the stable secondary structures, with the stem-length options set to a minimum of 20 bp. The promoters and Rho-independent terminators were analyzed using Promoter 2.0, RNAMotif and Erpin to further assess the resulting sequences. Finally, the nucleotide sequences of potential sRNAs were further inspected according to the results of RNA-seq deep sequencing. Each predicted sRNA was used to query the Rfam database to capture any annotated RNA regulatory elements.

Rhizobium strains were cultured in tryptone yeast medium (TY) or AMS medium (FAM2 minimal salts medium supplemented with 0.1% vol/vol aniline) at 28°C, and Escherichia coli strains were grown in Luria-Bertani (LB) at 37°C. When required, antibiotics were added at the following concentrations: ampicillin (Ap) 100 μg/ml; gentamicin (Gm) 20 μg/ml; streptomycin (Str) 20 μg/ml and kanamycin (Kan) 50 μg/ml.

In the Northern blotting analyses, free-living rhizobial cells were cultivated in TY at 28°C, and cells were harvested in the log phase (OD₆₀₀ = 0.4–0.5). When required, the growing log cultures were subjected to the following stress conditions as described previously (Madhugiri et al., 2012; Klahn et al., 2015): oxidative stress (2 mM H₂O₂ for 10 min), salt stress (4 M NaCl for 30 min), acid stress (pH5 for 90 min), alkali stress (pH9 for 90 min), cold stress (20°C for 90 min) or heat stress (37°C for 90 min). For micro-oxygen stress, log phase cells growing in TY (OD₆₀₀ = 0.4–0.5) were shifted to a low O₂ environment by continuous sparging of the culture with purified N₂ for 24 h (N₂ ≥ 99%). Controls were non-stressed exponentially growing cultures. For symbiotic conditions, surface-sterilized seeds of A. sinicus L. cultivar XY202 were cultivated in pots filled with sterile sand and irrigated with nitrogen-free Fahraeus solution, and cultivation was performed in a green house with a 16 h light and 8 h darkness cycle at 22°C and 20°C, respectively. Root nodules were harvested 28 or 50 days after M. huakuii 7653R inoculation.

Total RNA was extracted from M. huakuii 7653R to detect the selected sRNAs by Northern blotting. Following exposure to stress conditions, cells were quickly centrifuged, resuspended in TRIzol Reagent (Roche, USA), and immediately frozen in liquid nitrogen and stored at –80°C. Plant nodules were resuspended in TRIzol Reagent and frozen in liquid nitrogen, then immediately subjected to grinding using a baked mortar and pestle to crush the nodule tissues. Isolation of total RNA and Northern blotting were conducted according to the manufacturer's protocols. Crude RNA extracts were purified with DNaseI to remove potential genomic DNA, and the resulting RNA was stored at –80°C. Northern hybridizations were performed with end-labeled [γ^−32P]ATP DNA probes (Table S7). Hybridization signals were monitored using a radioisotope imaging system (Fujifilm, FLA 5100). The 5S rRNA was used as a control, and hybridization band intensities were relatively quantified with the Biorad Quantity One software.

Total RNA was extracted using an RNeasy Mini Kit (Qiagen) following the manufacturer's instructions to carry out the RNA deep sequencing. The RNA integrity number (RIN) of the total RNA samples was checked using an Agilent Bioanalyzer 2100 (Agilent). Qualified total RNA was further purified using an RNeasy micro kit (QIAGEN) and RNase-Free DNase Set (QIAGEN) according to the manufacturer's instructions.

Total RNA, including small-sized RNAs from 28 dpi (days post-infection) nodules, 50 dpi nodules and free-living cells at OD₆₀₀ = 0.8, was isolated and separated into small (< 100 nt) and long (> 100 nt) fractions using a miRNeasy Mini Kit (Qiagen) or a mirVana miRNA Isolation Kit (Ambion), respectively, according to the manufacturer's instructions. Total purified RNA was then treated with tobacco acid pyrophosphatase (TAP) (2.5 U/5 μg RNA in 50 μl for 2 h) to maintain the start regions of the transcription products. The quantity and quality of the processed RNA samples were assessed using of a NanoDrop ND-1000 spectrophotometer and Agilent 2100 Bioanalyzer. Qualified RNA was ligated to 5′- and 3′-RNA adapters by PCR, then the free adapters were removed and the remaining cDNA was purified. A Qubit® 2.0 Fluorometer and Agilent2100 instrument were used to monitor the abundance and size of cDNA, respectively. Solexa sequencing was performed using an Illumina HiSeq 2500. After that, all sequencing data were mapped by using the genomic sequence of M. huakuii 7653R to obtain the effective reads.

In order to figure out the potential function of the 16 new sRNA candidates identified in this work, we predicted the target genes for each of these sRNAs by using the webserver IntaRNA (http://rna.informatik.uni-freiburg.de/IntaRNA/Input.jsp). All of the encoded protein sequences of predicted target genes for each candidate sRNA were submitted to InterPro (http://www.ebi.ac.uk/interpro/sequence-search) for protein sequence analysis & classification, and the GO annotations (GO_ID) were extracted from the output. Further, the WEGO (Web Gene Ontology Annotation Plot) webserver (http://wego.genomics.org.cn/cgi-bin/wego/index.pl) was used to display the GO annotation results in histograms. The enriched GO terms of candidate small RNA target genes were summarized and demonstrated in a pie chart.

Searches for promoters in the M. huakuii 7653R genome were performed using Promoter 2.0 (http://www.cbs.dtu.dk/services/Promoter/) and BPROM 2.4.3.1111 (http://www.molquest.com/). Predictions of Rho-independent terminators were implemented with RNAMotif 3.0.7 (http://casegroup.rutgers.edu/casegr-sh-2.5.html), TransTermHP 2.09 (transterm.cbcb.umd.edu/), FindTerm 2.4.3.1111 (http://www.molquest.com/) and Erpin 5.4 (http://tagc.univ-mrs.fr/erpin/). Searches for known functional RNA sequences among the predicted sRNAs were conducted within the Rfam database using Rfam, version 12.0 (http://rfam.xfam.org/search) (Nawrocki et al., 2015). RNA secondary structures were analyzed using RNAfold 2.1.9 (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi). The target genes for the identified sRNAs were predicted using IntaRNA soft (http://rna.informatik.uni-freiburg.de/).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.