S. pyogenes MGAS315 [21] was used for most experiments. MGAS315 is a non-mucoid clinical strain isolated from a patient with streptococcal toxic shock syndrome, and its genome sequence is publically available [21]. S. pyogenes was routinely cultured in Todd-Hewitt medium (BBL) supplemented with 0.2% yeast extract (Difco) at 37°C in sealed tubes without agitation.

In this study, we employed three computational algorithms to increase prediction accuracy: eQRNA, RNAz, and sRNAPredict. The computational approaches are illustrated in Figure 1 and the sources and references of the computational algorithms used in this study are listed in Table 1. The genomic information of the following eight streptococci necessary to run these algorithms, such as genome sequences, loci and names of open reading frames (ORFs), tRNAs, tmRNA, rRNAs and IGRs, was downloaded from the NCBI website (ftp://ftp.ncbi.nih.gov/genomes/Bacteria/): S. pyogenes MGAS315, S. equi subsp. zooepidemicus MGCS10565, S. mutans UA159, S. suis 05ZYH33, S. sanguinis SK36, S. gordonii str. Challis substr. CH1, S. pneumoniae CGSP14, S. agalactiae NEM316.

To run sRNAPredict, the coordinates and the orientations of 1865 ORFs, 67 tRNAs, 18 rRNAs and 1tmRNA in the genome of S. pyogenes MGAS315 were marked and then 1587 IGRs were extracted by Intergenic Sequence Inspector (ISI). To predict the location of Rho-independent terminators of S. pyogenes MGAS315, the RNAmotif program and the TransTerm database were used. RNAMotif searches RNA sequences that match a “motif” describing the interactions of secondary structures, which are defined via a pattern language whose symbols represent helices and single stranded stretches. Matches can be ranked by applying scoring rules that may provide finer distinctions than just matching to a profile. Rho-independent terminators predicted by TransTerm were obtained from the TransTerm website at University of Maryland. Putative terminators whose probabilities (confidence) were greater than 90% were chosen for sRNAPredict analysis. Parameters applied to sRNAPredict analysis are as follows; the minimum distance of predicted terminator from the end of an upstream ORF: 30 nts, the maximum length of gap between a putative terminator and a region of sequence conservation: 20 nts, the values for the minimum and maximum length of putative sRNAs: 30 nts and 550 nts.

For the prediction with eQRNA and RNAz, a file containing intergenic regions (IGRs) of S. pyogenes MGAS315 was generated using ISI. The number of the extracted IGRs from S. pyogenes MGAS315 genome by ISI was 3166 sequences from both DNA strands, totaling about 0.681 MB which represents 17.9% of the full genome. The average sequence length was 138 nts, with the longest one being 1501 nts and the shortest one being 1 nt. Then, the WU-BLAST 2.0 program was used for genome wide sequence homology analysis between IGRs of S. pyogenes MGAS315 and the genomes of the other seven streptococci. The IGRs with a length >12 nts were used as queries. The output data of the BLAST analysis gave pairwise alignments between the query sequences, that are the IGRs of S. pyogenes MGAS315, and the subject sequences, that are the genomic segments of the other seven streptococci. These pairwise alignments of 2118 comparisons were filtered with the parameters of E-values <0.00001 and length >30 nts (Table S1). These pairwise alignments were scanned by CLUSTALW to produce multiple alignments. From the multiple alignments, the long sequences with the length >550 nts were removed. Alignments contained sequences from both DNA strands, and candidates were selected when a signal was identified from either of the two strands. RNAz and eQRNA predictions from the alignments were incorporated into a single predicted RNA locus on the genome. An additional set of alignments was obtained using 69 known RNAs (68 tRNAs and one tmRNA) as queries. These RNA alignments were generated using the same BLAST parameters, and were used to evaluate the sensitivity and specificity of those computational analyses. The multiple sequence alignments that were formatted to CLUSTALW data were used as the input source for RNAz. The CLUSTALW program uses FASTA format as input data. Thus, we transformed the pairwise sequences from BLASTN to FASTA format to execute CLUSTALW. Both RNAz and eQRNA used the window size of 150 nts and the window slide increment of 50 nts. To test sensitivity and specificity of this approach, we used 68 known tRNAs and tmRNA as controls. Through the RNAz and eQRNA analyses, 65 out of 68 S. pyogenes tRNAs and tmRNA were identified (95.6% sensitivity) (Table S2). To test the specificity of the RNAz and eQRNA analyses, we shuffled the sequences of the tRNAs and tmRNA and estimated the false positives if any shuffled sequence was considered as sRNA. For shuffling of sequences, the RNAs were divided into several groups and the sequence locations of the groups (each group has 20 nts and the size of an ordinary tRNA is about 70 nts) were exchanged. This shuffling keeps the sequence conservation but not the conservation of secondary structures. The number of false positives obtained from this process was 1/68 (1.47%), so the specificity was 98.5% (Table S2). RNAz analysis shows similar sensitivity but higher specificity than the eQRNA analysis, and the combination of RNAz and eQRNA analyses increased specificity compared to the individual analysis, as expected (Table S2).

The computer algorithms used in this study and the result of each process can be downloaded from www.cs.siu.edu/~nyu/research.htm.

Total S. pyogenes RNA was extracted using the combination of the miRNeasy kit (Qiagen) and the FastPrep beadbeater (MP biomedicals). An S. pyogenes cell pellet from 10 ml culture was resuspended in 700 μl of the Qiazol lysis reagent provided in the miRNeasy kit and transferred to a Lyse Matrix B blue cap tube (MP biomedicals). Cells were then lysed by the beadbeater, FastPrep 24 (MP biomedicals) at the speed of 6.0 for 40 seconds twice. The remaining procedure for RNA extraction followed the manufacturer's protocol of the miRNeasy kit. The A260/A280 ratio of the extracted RNA was measured with NanoDrop (Thermo Scientific) to determine the RNA concentration and purity (accepted if >1.8). The extracted RNA was mixed with 1 μl of RNasin (Promega Recombinant RNasin Ribonuclease Inhibitor, 40 u/μl), and treated with RNase-Free DNase (Promega DNase I, 1u/µl) according to the manufacturer's protocol.

S. pyogenes MGAS315 was grown in THY medium at 37°C and harvested with centrifugation (at 6,000 × g for 3 min) at the exponential phase (OD₆₀₀, Optical density at 600 nm,  = ∼0.5), the early stationary phase (OD₆₀₀  = ∼1.2) or the late stationary phase (cells were grown for 3 hrs more from the early stationary phase, OD₆₀₀  = ∼1.5). Then, Total RNA was extracted as described above. The extracted RNA (20 µg) was mixed with loading dye containing 50% (v/v) formaldehyde, loaded onto denaturing polyacrylamide gel (6% polyacrylamide (acrylamide: bis-acrylamide  = 29∶1) with 7% urea) pre-run at 400 V for an hour, and electrophoresed at 300 V. For an RNA size marker, 5 ng of Low Range SSRCNA Ladder (New England BioLab, 500 μg/ml) was loaded in a well. The separated RNA was transferred onto a nylon membrane (Zeta Probe Blotting Membrane, Bio-Rad) with a semi-dry electroblotter (Bio-Rad Trans-Blot SD Transfer Cell) at 400 mA for two hours at 4°C. The RNA on the nylon membrane was cross-linked to the membrane with 1-ethyl-3-(3-dimethylaminopropyl)-carbodiimide (EDC) [22]. The membrane was then prehybridized with 5X SSPE/2% SDS hybridization buffer for 30 min and hybridized for 18 hrs with the same buffer containing a 40 nM single stranded DNA oligonucleotide probe (Table S3). The probe had been ³²P-labeled at 5′-end with γ-³²P ATP (10 mCi/ml, PerkinElmer) by T4 polynucleotide kinase (Epicentre Technologies). The probes were designed to bind at the center of putative small RNAs predicted by the computational algorithms. Because the list of sRNA candidates from the computational analysis did not provide information on which DNA strand encodes each sRNA candidate, we designed and tested two probes annealing to each strand. The sequence of each probe is shown in Table S3. After hybridization, the membrane was washed with wash solutions and exposed to X-ray film (autoradiography).

S. pyogenes MGAS315 was grown in THY medium at 37°C until it reached an OD₆₀₀ of 0.5. The cells were harvested with centrifugation and total RNA was extracted as described in “RNA Extraction from S. pyogenes MGAS315”. The extracted RNA was converted to cDNA using reverse transcriptase (ImProm II Reverse Transcriptase, Promega) according to the manufacturer's protocol. Briefly, RNA (2 µg) was mixed with 500 ng of random primers (Promega) and adjusted to 5 μl with RNase free water. The RNA mixture was incubated at 65°C for five minutes and then chilled at 4°C for five minutes. A mix containing 2.4 μl of MgCl₂ (25 nM), 4 μl of ImProm-II 5× reaction buffer, 1 μl of ImProm-II Reverse Transcriptase, 1 μl of dNTPs (10 nM) and 6.6 μl of RNase Free water was added to the RNA mix and incubated at 25°C for five minutes, followed by 42°C for one hour, and then heat inactivated at 70°C for fifteen minutes. Regular PCR was performed with the cDNA as a template. RNA (without reverse transcriptase reaction) was used as a control to confirm that a PCR product was not from chromosomal DNA contamination (RNA was rejected before producing cDNA if PCR amplification performed with the RNA template indicated the presence of contaminating DNA). PCR products and a DNA ladder (1 Kb Plus DNA Ladders, Invitrogen) were electrophoresed on a 2% agarose gel.

Primers for RT-PCR were designed with the parameters of 1 bp GC clamp at 3′ end, 20 nts size, 100–200 nts product size, and 60°C melting temperature (Table S4).

The transcriptional start and stop sites of sRNA candidates were determined using circular RACE (Rapid Amplification of cDNA ends) as described elsewhere [23]. Briefly, RNA was extracted from S. pyogenes cultures in the exponential phase (OD₆₀₀ of 0.5) as described above. The RNA was treated with Tobacco Acid Pyrophosphatase (TAP) (Epicentre) to remove pyrophosphate from the 5′ end. The 5′ end was then ligated to the 3′ end with T4 RNA ligase (Epicentre) to make circular transcripts. The circular transcripts were reverse transcribed using gene specific primers to make first strand cDNA. The first strand cDNAs were amplified with PCR. The PCR products were cloned into pGEM-T (Promega) and sequenced to determine transcriptional start and stop sites. Primers for the circular RACE used in this study are listed in Table S5.

To create a deletion mutant of SSRC21, ΔSSRC21cat, 144 bps of the internal part of SSRC21 was deleted and replaced with a chloramphenicol acetyltransferase gene (cat) [24] through a double cross over homologous recombination. To achieve this, first, two DNA fragments flanking SSRC21 on each side were amplified with primers and joined to delete SSRC21. The primers 5outSSRC21far and 3inSSRC21 were used to amplify 1.23 kbps DNA fragment upstream of SSRC21, and 5inSSRC21 and 3outSSRC21PstI were used to amplify 0.69 kbps DNA fragment downstream of SSRC21 (Table 2). These two fragments were digested with XmaI, ligated, and PCR amplified with the primers 5outSSRC21 and 3outSSRC21PstI (Table 2). The PCR-amplified 1.36 kbps DNA was sequenced to confirm the deletion of SSRC21 and inserted into vector pJRS233 [25]. The deleted SSRC21 was replaced with the cat gene (0.98 kbps) that was PCR-amplified from pABG5 [26] with the primers of 5catXmaI and 3catXmaI (Table 2). The cat gene flanked with ∼0.7 kbp streptococcal DNAs on each side was amplified with the same primers used to amplify the 1.36 kbp PCR fragment (5outSSRC21 and 3outSSRC21PstI), and then used to transform MGAS315 (wild type) by electroporation. The mutant showing chloramphenicol resistance was selected and the chromosomal structure was confirmed by PCR.

Total RNA was extracted from S. pyogenes cultures in the exponential phase (OD₆₀₀ of 0.5) as described above, and submitted to Otogenetics Corporation (Norcross, GA USA) for RNA-Seq assays. Briefly, the integrity and purity of total RNA were assessed using Agilent Bioanalyzer and OD260/280. Up to 5 μg of total RNA was subjected to rRNA depletion using the RiboZero Meta-Bacteria kit (Epicentre Biotechnologies, Madison, WI USA, catalog # MRZMB126) and cDNA was generated from the depleted RNA using the NEBNext mRNA Sample Prep kit (New England Biolabs, Ipswich, MA USA, catalog# E6110). cDNA was profiled using Agilent Bioanalyzer, and subjected to Illumina library preparation using NEBNext reagents (New England Biolabs, Ipswich, MA USA, catalog# E6040). The quality and quantity and the size distribution of the Illumina libraries were determined using an Agilent Bioanalyzer 2100. The libraries were then submitted for Illumina HiSeq2000 sequencing according to the standard operation. Paired-end 90 or 100 nucleotide (nt) reads were generated and subjected to data analysis using the platform provided by Center for Biotechnology and Computational Biology (University of Maryland, College Park, MD USA) as previously described [27]. The data sets generated from RNA-Seq were mapped against GenBank AE014074 (http://www.ncbi.nlm.nih.gov/nuccore/21905618) with Bowtie2 (V2.0.0.5). Hits on regions defined by GenBank were then counted with bedtools. To determine the difference of gene expression between samples, EdgeR (Empirical analysis of digital gene expression data in R) was used. We deposited our RNA seq dataset to NIH Short Read Archive with the accession number SRP020234.

Along with experimental strategies based on shotgun cloning or microarray methods, computational predictions and validation with Northern blot have been a popular method used to identify many sRNAs (for a review, refer to [28]). Most computational algorithms developed for genome-wide screening of small RNAs are based on ‘sequence or structural conservation’ among closely related species. The algorithms that seek sequence conservation such as sRNAPredict [29], [30], and GMMI [31] search first for transcriptional signals such as promoters and terminators, and then examine nucleotide conservation. Since the sequence homology information is based only on the primary structure of RNA, the accuracy of this method may not be sufficiently adequate. Hence, some algorithms seek for phylogenetic conservation of secondary structure and/or thermal stability. This type of algorithm includes Pfold [32], MSARI [33], eQRNA [34], RNAz [35], etc. In this study, we employed a combination of three algorithms that search for different forms of conservation as a way to increase prediction accuracy: sRNAPredict [29], [30], eQRNA [34], and RNAz [35] (Figure 1). The sRNAPredict algorithm uses the information of the location of transcriptional signals and primary sequence conservation of intergenic regions (IGRs). The eQRNA algorithm examines the conservation of secondary structures of RNA for prediction. It identifies base substitution patterns in pairwise alignments likely corresponding to a conserved RNA secondary structure rather than to a conserved coding frame or other genomic features. RNAz even measures thermodynamic stability, which is normalized with respect to both sequence length and base composition in addition to RNA consensus secondary structure. The eQRNA and RNAz algorithms are comparatively strict methods for secondary structure conservation analysis compared to other methods such as Pfold [32] and MSARI [33]. To search for those conservations between closely-related streptococcal species, genome sequences and annotations of the following eight streptococci were used. S. pyogenes MGAS315, S. equi subsp. zooepidemicus MGCS10565, S. mutans UA159, S. suis 05ZYH33, S. sanguinis SK36, S. gordonii str. Challis substr. CH1, S. pneumoniae CGSP14, S. agalactiae NEM316.

Each algorithm of sRNAPredict, eQRNA, and RNAz respectively predicted 191, 312, and 187 intergenic genomic segments as putative sRNAs in S. pyogenes MGAS315. Among these predicted ones, the sequences immediately upstream of annotated open reading frames (ORFs) were removed from the candidates because they most likely correspond to putative riboswitches or other cis-regulatory elements. In addition, predicted candidates in prophages, all of which were located next to the integrase genes, were removed. Then, the intergenic genomic segments predicted by any two algorithms were considered as putative sRNA candidates. The final number of sRNA candidates left from these processes was 45 (Figure 1B, Table 3, and Table S6). Encouragingly, Pel, FasX, and RivX, which have been previously studied sRNAs, were included in these 45 streptococcal small RNA candidates (SSRCs). Among the 45 putative candidates, six (SSRC 8, 12, 15, 30, 32, 33) were predicted by all the three algorithms (Table 3).

Before applying Northern blotting to verify the expression of streptococcal sRNA candidates, we screened the candidates with RT-PCR. RT-PCR can at least detect the presence or absence of their transcripts, even though it cannot distinguish whether the expressed transcripts are cis-elements attached to mRNAs or independently expressed sRNAs. In the RT-PCR analysis, we did not detect the expression of SSRC 1, 3, 6, 7, 22, 24, 28, 39, 40, so we eliminated them from the sRNA candidates.

Then, we examined the remaining candidates with Northern blotting. Since most S. pyogenes sRNAs previously described were highly expressed in the exponential phase of growth [20], RNA was extracted from cells at the exponential phase and used for Northern blotting. The predicted sRNA candidates could be expressed from either strand of DNA, so we designed and tested two probes for each candidate (Table S3).

Among the candidates, the following 14 SSRCs showed consistent signals on Northern blots: SSRC4 (FasX), SSRC8, SSRC10, SSRC12 (Pel/sagA), SSRC13, SSRC21, SSRC27, SSRC29, SSRC30, SSRC31, SSRC32, SSRC34, SSRC38, and SSRC41 (Figure 2). Some candidates (detected with 5′ probes in the Table S3) were expressed from the top DNA strand of the MGAS315 chromosome sequence (GenBank: AE014074.1, http://www.ncbi.nlm.nih.gov/nuccore/21905618?report=fasta) and the others (detected with 3′ probes in the Table S3) were from the complementary DNA strand; The candidates, SSRC21, SSRC27, SSRC29, SSRC30, SSRC32, and SSRC34 were detected with their 5′ Northern blot probes, and SSRC8, SSRC10, SSRC13, SSRC31, SSRC38, and SSRC41 were detected with their 3′ probes. FasX and Pel/sagA were detected as expected. However, RivX was not detected probably because of its extremely low expression in the wild type [15]. Since FasX gave a constant signal, it was used as a control throughout the Northern blot analysis.

The size of each SSRC detected by Northern blot was calculated by measuring traveled distance of each band from a well of a polyacrylamide gel and compared to a standard curve obtained from the RNA ladder that was run with each sRNA side by side. The calculated approximate sizes of the SSRCs are shown in Figure 2. The calculated size of FasX (202 nts) was very close to the actual size (203 nts). The probe designed for Pel/sagA-transcript detected a 133 nts band only, which is the size close to the sagA-transcript, not Pel (459 nts).

Sequence analyses based on the predicted SSRC sequences (Table S6) and next-generation sequencing, RNA-Seq, revealed that SSRC13, SSRC31 and SSRC32 are probably not trans-acting regulatory sRNAs. SSRC13 and SSRC31 appear to be cis-elements in RNA-Seq analysis. The sequence of SSRC13 contains a T-box leader element, which is typically found upstream of aminoacyl-tRNA synthetase genes and some amino acid biosynthesis genes and involved in the regulation of those genes' expression by forming a transcription anti-terminator when uncharged tRNA binds the leader sequence [36]. SSRC13 is located upstream of phenylalanyl-tRNA synthetase subunit alpha. SSRC32 appears to be ribonuclease P (RNase P) RNA, RnpB. RNase P is a ribozyme, which is composed of two components, RnpA (protein) and RnpB (RNA), and cleaves the 5′ leader sequence of precursor tRNAs to produce mature tRNAs [37].

We performed circular RACE (Rapid Amplification of cDNA Ends) to determine the 5′ and 3′ ends of selected small RNA candidates. Based on the calculated sizes, Northern blot probe-binding sites, and/or putative transcriptional signatures (promoters and Rho-independent terminators), we designed primers for each SSRC. From this analysis, we could determine the sequences of SSRC8, SSRC10, SSRC21, SSRC29, SSRC34 and SSRC41 (Figure 3). All the determined sequences contained a Rho-independent transcriptional termination signal (hairpin structure ending with or followed by thymidines). By examining the sequence of promoter regions, we could map putative −10 and −35 promoter sequences of SSRC10, SSRC21, SSRC29, and a putative CovR-binding sequence [38] upstream of SSRC34 (Figure 3).

The three computer algorithms used in this study compared sequences between closely related streptococcal species to search for putative sRNAs, so homologs of the newly discovered small RNA would exist in not only S. pyogenes but also other streptococcal pathogens. Thus, we searched for the homologs of the SSRCs in the other streptococci whose genome sequences are available publically. Each streptococcal genome was blasted against the SSRC sequences determined through circular RACE. As expected, homologs of the SSRCs existed in many streptococcal bacteria, and more SSRC homologs were found in streptococci more closly related to S. pyogenes (Table 4) [39]. For example, S. dysgalactiae subsp. equisimilis, S. equi, S. parauberis, S. salivarius, S. thermophilus had more than 4 homologs. Notably, S. dysgalactiae subsp. equisimilis, which is a beta-hemolytic streptococcus very closely related to S. pyogenes, contained all the six SSRCS. Exceptionally, some closely related streptococci such as S. canis and S. iniae did not contain any homolog of the SSRCs. The most frequent SSRC found in streptococci among the six SSRCs was SSRC10.

We determined the intracellular abundance of the putative trans-acting SSRCs at the exponential (EX), early stationary (ES), and late stationary (LS) growth phases of cells through Northern blotting (Figure 4). Agreeing with a previously reported result, FasX was most abundant at the exponential phase [20]. Most SSRCs exhibited variation of abundance between growth phases and were expressed abundantly at the exponential and early stationary phase and least abundantly at the late stationary phase. The abundance of SSRC 8, 10, 21, 29, 30, and 38 was dramatically reduced at the late stationary phase, showing a similar pattern to that of the FasX transcript over the course of growth. On the other hand, SSRC 34 and 41 exhibited similar abundance throughout all growth phases.

To screen putative target mRNAs of an sRNA detected in this study, we employed differential RNA sequencing analysis. SSRC21 was chosen for this analysis because its expression was changed between growth phases, and initial computational prediction suggested that SSRC21 might influence the expression of several virulence factors. The transcript abundance profile in an SSRC21 deletion mutant was determined through Next-Generation Sequencing, RNA-Seq, and compared to that of the wild type MGAS315. To construct the SSRC21 deletion mutant, ΔSSRC21cat, the SSRC21 gene in the chromosome was replaced with the chloramphenicol acetyl transferase gene, cat. In the RNA-Seq result, 142 transcripts exhibited differential abundance with the criteria of fold change in the mutant over the wild type, FC, greater than 2 or less than −2, p<0.01, and false discovery rate, FDR, <5% (Table S7). The transcripts (with known putative function) showing notable differential abundance between the mutant and the wild type were those of ntp genes encoding V-type ATPase subunits and their putative regulator (SpyM3_0113 through SpyM3_0122; fold changes in the mutant over wild type, FC, 4.2∼7.4). Some transcripts encoding virulence factors also showed differential abundance: streptokinase A, FC 2.6; C5A peptidase, FC −2.1; M protein, FC −2.4; streptococcal phospholipase A2, FC −4.5).