Differential evolution in 3′UTRs leads to specific gene expression in Staphylococcus

Abstract The evolution of gene expression regulation has contributed to species differentiation. The 3′ untranslated regions (3′UTRs) of mRNAs include regulatory elements that modulate gene expression; however, our knowledge of their implications in the divergence of bacterial species is currently limited. In this study, we performed genome-wide comparative analyses of mRNAs encoding orthologous proteins from the genus Staphylococcus and found that mRNA conservation was lost mostly downstream of the coding sequence (CDS), indicating the presence of high sequence diversity in the 3′UTRs of orthologous genes. Transcriptomic mapping of different staphylococcal species confirmed that 3′UTRs were also variable in length. We constructed chimeric mRNAs carrying the 3′UTR of orthologous genes and demonstrated that 3′UTR sequence variations affect protein production. This suggested that species-specific functional 3′UTRs might be specifically selected during evolution. 3′UTR variations may occur through different processes, including gene rearrangements, local nucleotide changes, and the transposition of insertion sequences. By extending the conservation analyses to specific 3′UTRs, as well as the entire set of Escherichia coli and Bacillus subtilis mRNAs, we showed that 3′UTR variability is widespread in bacteria. In summary, our work unveils an evolutionary bias within 3′UTRs that results in species-specific non-coding sequences that may contribute to bacterial diversity.

About cis-acting 3 UTRs, we previously showed that the long 3 UTR of the icaR mRNA, which encodes the main repressor of the PIA-PNAG exopolysaccharide synthesis in Staphylococcus aureus, interacts with the 5 UTR of the same mRNA to inhibit IcaR translation. This occurs through base-pairing of the anti-SD motif (UCCCC sequence located in the 3 UTR) at the ribosome binding site (RBS), impairing ribosome accessibility by creating a double stranded RNA substrate that is cleaved by the endoribonuclease III (RNase III) (5). Similar 5 -3 UTR interactions were recently described in the hbs mRNA of Bacillus subtilis. The cleavage of the hbs 5 UTR by RNase Y is countered by a 14-nt perfect base-pairing interaction between the 5 and the 3 UTR, which partially overlaps with the cleavage site (8). Other bacterial 3 UTRs are often used as direct entry points for ribonucleases to initiate mRNA degradation. This is the case for the hmsT mRNA, which encodes a protein that modulates c-di-GMP synthesis, required for allosteric activation of polysaccharide production in Yersinia pestis biofilm formation. Polynucleotide phosphorylase (PNPase)-dependent cleavage of the hmsT 3 UTR affects the turnover of the mRNA and, ultimately, HsmT expression (6,9). Additionally, several AUrich 3 UTRs of Y. pestis are targeted by PNPase when the transcriptional terminator (TT) is Rho-dependent (9). In Salmonella, HilD, one of the main transcriptional activator factors of the pathogenicity island 1 (SPI-1) (20), is processed through its 3 UTR by RNase E and PNPase (21). Similarly, the aceA mRNA in Corynebacterium glutamicum, which encodes the isocitrate lyase protein that is part of the glyoxylate cycle (22), is processed by endoribonucleases E/G (RNase E/G) at the 3 UTR (4). On the contrary, there are examples of 3 UTRs that prevent mRNA degradation. When Escherichia coli grows under iron-limiting conditions, the apo-AcnB protein interacts with a stem-loop structure located at the 3 UTR of its own mRNA, and prevents RyhB-induced degradation. This occurs because the apo-AcnB binds close to the RNase E cleavage site, which is essential for acnB mRNA degradation (7).
Regarding trans-acting regulatory RNAs, S. aureus RNAIII was for decades the sole example of a prokaryotic mRNA with a regulatory 3 UTR. RNAIII consists of a 5 UTR of 84-nt, a short CDS encoding the ␦-toxin (hld), and a long structured 3 UTR (352-nt) that controls the temporal expression of several virulence factors and the transcriptional repression of Rot exotoxins by targeting their mRNAs (23,24). More recently, several small RNAs (sR-NAs) derived from 3 UTRs have been described, some of which are expressed from alternative mRNA-internal promoters (Type I 3 UTR-derived sRNAs), while others are the consequences of specific mRNA processing events (Type II 3 UTR-derived sRNAs) (18).
Previous studies have shown the sequence conservation of some 3 UTRs to be variable depending on the bacterial species. For example, the 3 UTR of dapB mRNA, which produces the type I dapZ 3 UTR-derived sRNA, is not conserved among enterobacterial species (except for the R1 seed and the Rho-independent terminator) (11). Similarly, the sdhX sRNA, which is generated by processing the 3 end of the sdh-suc operon by RNase E in E. coli and Salmonella, also shows nucleotide variations between different enterobacterial species (16,17). Interestingly, sequence differences in the sdhX sRNA could establish an mRNA target divergence between Salmonella and E. coli (16).
Although, these and other pioneering studies have indicated the importance of 3 UTRs as a new class of post-transcriptional regulatory elements controlling relevant physiological processes in bacteria, several questions remain unanswered. Are 3 UTR sequences preserved within and between bacterial species? How often do conserved genes from closely related bacteria contain different 3 UTRs and how is 3 UTR variability originated? Do differences in 3 UTRs between orthologous genes have consequences in their expression at the protein level?
In this study, we aimed to evaluate the evolutionary relationship between the CDSs and their corresponding 3 UTR sequences in bacteria. To this end, we performed genomewide comparative analyses of orthologous mRNAs among closely related species of the genus Staphylococcus. As a result, we discovered that most of the 3 UTRs from these mRNAs were rather variable at the sequence level. In most cases, the inter-species nucleotide conservation of mRNAs was lost downstream of the coding sequence, suggesting a large variability within the 3 UTRs. In vivo mapping using RNA-seq of the transcript boundaries in three species of the genus Staphylococcus confirmed that 3 UTRs vary in length and sequence when compared with the S. aureus transcriptomic data. We demonstrated that such differences may have an impact on the expression of orthologous genes and, thus, the functionality of the 3 UTRs could differ depending on the analysed species. Finally, by extending the 3 UTR conservation analysis to other bacteria genera, we found similar sequence variations. In summary, this study proposes that 3 UTRs from orthologous genes may be selectively targeted by evolution to create sequence differences that contribute to diversity in bacteria.

Nucleotide conservation analysis
The conservation of UTRs among orthologous monocistronic mRNAs from phylogenetically-related bacterial species was determined by performing blastn comparisons using the Microbial Nucleotide BLAST tool (https: //blast.ncbi.nlm.nih.gov/). Monocistronic mRNA boundaries were manually annotated by visualizing the S. aureus transcriptomic maps from previous studies (5,(25)(26) in the Jbrowse application (27). For the batch comparison of S. aureus 3 UTRs, a file including the query sequences in FASTA format was used. Each sequence included the annotated 3 UTR and the last 200 nt of the CDS of a monocistronic operon (Supplementary Figure S1). This restriction was applied to normalize the starting point for each of the UTRs from the blastn results, facilitating further analysis and plotting. We considered that 200 nt from the CDS were enough to properly identify their corresponding orthologous genes. Similarly, for batch comparison of S. aureus 5 UTRs, a FASTA file including the annotated 5 UTRs and the first 200 nt of the CDSs were used. In order to normalize the starting point of the CDS, we used reversecomplement query sequences for 5 UTR blastn comparisons (Supplementary Figure S1). The UTR conservation length was registered as the last nucleotide from the query sequence showing a positive alignment with the region of the genome under comparison. For this purpose, the information contained in the blastn result TXT files was transformed into tables using a custom R script. Nucleotide alignments that started at values higher than 100, finished at values lower than 100, or were shorter than 80 nucleotides were not considered. The UTR conservation lengths of each Staphylococcus species were plotted against the reference UTR lengths of S. aureus (Supplementary Figure S1). The square correlation coefficient (R 2 ) was calculated in R from the resulting fitted linear model. This analysis was also applied to compare the mapped 3 UTRs of B. subtilis (28) and E. coli (29) with phylogenetically-related species of the Bacillus genus and Enterobacteriaceae family, respectively. Note that in the case of B. subtilis, 3 UTRs longer than 150 nt were not annotated (28).

RNA sequencing and data analysis
In order to perform the transcriptomic maps of the S. simiae CCM 7213T, S. capitis SK14 and S. epidermidis RP62A, strains were grown in TSBg at 37 • C and 200 rpm until exponential phase was reached. Total RNAs were extracted as previously described (30). RNA sequencing and preliminary data analysis was carried out by the Stab Vida company. The RNA-seq reads from the S. simiae CCM 7213T, S. capitis SK14 and S. epidermidis RP62A samples were aligned using the Rockhopper program (31) and the complete genome sequences from S. simiae NCTC13838 (NZ LT906460.1), S. capitis AYP1020 (NZ CP007601.1) and S. epidermidis RP62A (NC 002976.3) as references, respectively. The obtained read coverage files (.wig) from the plus and minus strands were converted to BigWig (.bw) files with the wig2BigWig program. The generated .bw files and the gene annotation files, including the positions of transcriptional terminators (TTs) predicted by TransTerm HP program (32), were loaded into a public web server based on Jbrowse (27) (http://rnamaps.unavarra.es/). The 3 UTR lengths of the orthologous monocistronic mRNAs were manually annotated and plotted against the lengths of previously described S. aureus 3 UTRs (5).

Simultaneous mapping of 5 and 3 ends from RNA molecules
The mapping of 5 and 3 ends of mRNAs was performed using a modified version of the Rapid Amplification of cDNA Ends (mRACE) method (33). Briefly, total RNA samples were treated with the Cap-Clip Acid Pyrophosphatase (Tebu-Bio) following the manufacturer's recommendations. After 1 h of incubation, RNAs were extracted by phenol-chloroform and precipitated with sodium acetate and cold ethanol. Serial dilutions of Cap-Clip and non-Cap-Clip treated RNA were ligated using the T4 RNA Ligase I (New England Biolabs) at 16 • C ON. RT-PCRs for all the ligated RNA dilutions were carried out using the SuperScript™ III One-Step RT-PCR System with the Platinum Taq DNA polymerase (Invitrogen) and the outward primers A and B (Supplementary Table S3). For the mapping of the 3 UTR of rpiRc from S. epidermidis and S. capitis primers RACE-3XFLAG-A and RACE-RpiRc-sau-B were used. PCR products were run in 2.5% agarose gels and bands of the expected size were purified and ligated into the pGEM-T Easy vector (Promega). The resulting reactions were used to transform E. coli XL1-Blue cells (Stratagene). Ten isolated white colonies were then analyzed by Sanger sequencing. Transcript boundaries were determined by blastn analysis and nucleotide frequencies representing the 5 and 3 ends, respectively, were registered as percentages (%).

Chromosomal mutagenesis
The mutants generated in this study (Supplementary Table  S1) were obtained by a two-step homologous recombination that exchanges a specific chromosomic region by the mutant allele present in the pMAD plasmid (34), as previously described (35). The marker-less mutants were verified by PCR using oligonucleotides E and F (Supplementary Table S3) and Sanger sequencing.

Plasmid constructions
Most of the plasmids used in this study were engineered as previously described (36). Briefly, PCR fragments were amplified from chromosomic or plasmidic DNA with the DreamTaq DNA polymerase or Phusion High-Fidelity DNA Polymerase (Thermo Scientific), using the oligonucleotides listed in Supplementary Table S3. The PCR products were run and purified from agarose gels using the NucleoSpin ® Gel and PCR Clean-up Macherey-Nagel kit and ligated into the pJET 1.2 vector (Thermo Scientific). The resulting plasmids were used to transform E. coli XL1-Blue cells. These were then purified from overnight (ON) cultures with the NucleoSpin ® Plasmid Macherey-Nagel kit and verified by Sanger sequencing. When required, DNA fragments were excised using FastDigest restriction enzymes (Thermo Scientific) and ligated into the appropriate vector with the Rapid DNA ligation kit (Thermo Scientific). The final plasmids were introduced into S. aureus strains by electroporation, as previously described (37).
The pMAD plasmids used for chromosomal mutations were constructed by amplifying the flanking sequences (AB and CD) of the target regions using primers A/B and C/D (Supplementary Table S3). After cloning them into pJET they were digested and ligated into pMAD in a two-insert ligation process using BamHI, NheI and EcoRI (Supplementary Table S2 and S3).
The plasmids expressing chimeric icaR mRNAs were constructed taking advantage of conserved natural restriction sites located at the beginning of the icaR 3 UTRs. Specifically, the icaR mRNA regions were amplified by PCR using chromosomic DNA from the corresponding staphylococcal species and the specific oligonucleotides IcaR+1 and IcaR-Term (Supplementary Table S3). The resulting PCR products were inserted into pJET and digested with SpeI and EcoRI. In the case of S. simiae, restriction enzymes HincII and EcoRI were used. The digested fragments were ligated into the p 3xF IcaRm or pIcaRm (5), producing the plasmids carrying the chimeric icaR mRNAs listed in Supplementary Table S2.
The ftnA and rpiRc genes were labelled with the 3xFLAG tag sequence by overlapping PCRs. Specifically, the tagged ftnA gene was generated by linking two partially overlapping PCR fragments that were amplified using oligonucleotides Ftn-A, 3XFLAG-ftn-izqda, 3XFLAG-ftn-drcha and Term-ftn (Supplementary Table S3). A second PCR with oligonucleotides +1-ftn and Term-ftn (Supplementary  Table S3) was performed to amplify the 3xF ftnA mRNA region. The PCR fragment was ligated into pJET, digested with restriction enzymes BamHI and KpnI and inserted into pEW, generating the p 3xF FtnA plasmid (Supplementary Table S2). The pEW plasmid was constructed by transferring the EcoRI-NarI fragment (which includes the transcriptional terminator region) from pCN47 to the pCN40 plasmid (38). Similarly, the rpiRc gene was tagged by introducing the 3xFLAG sequence through overlapping PCRs using oligonucleotides +1-RpiRc, 3XFLAG-RpiRcizqda, 3XFLAG-RpiRc-drcha and Term-RpiRc (Supplementary Table S3). This PCR fragment was digested with BamHI and EcoRI and inserted into pEW, generating the p 3xF RpiRc plasmid (Supplementary Table S2).
To generate the plasmids expressing the chimeric tagged ftnA and rpiRc mRNAs (Supplementary Table S2), the corresponding 3 UTRs from different staphylococcal species were fused to S. aureus FtnA or RpiRc CDSs by overlapping PCRs. First, PCR fragments including the 5 UTR and the 3xFLAG CDS from the ftnA or rpiRc genes were generated by using the p 3xF FtnA and p 3xF RpiRc plasmids as templates and oligonucleotide pairs +1-ftn/CDS-stop-ftn and +1-RpiRc / CDS-stop-RpiRc as primers, respectively. Second, the 3 UTR regions from different staphylococcal species were amplified using their corresponding fw and rvs oligonucleotides (Supplementary Table S3). Finally, the PCR fragments including the 3xFLAG tagged CDS and the 3 UTR were fused by PCR using oligonucleotides +1 and rvs (Supplementary Table S3). Next, they were ligated into pJET and digested with BamHI/ KpnI and BamHI/EcoRI in the case of 3xF FtnA and 3xF RpiRc, respectively, and inserted into pEW.
The GFP reporter plasmids were constructed using the Listeria monocytogenes pAD-cGFP plasmid as a template (39). The hly 5 UTR and GFP sequences were amplified with primers Sal-GFP-fw and BcuI-TT-BamHI-GFP-rvs (Supplementary Table S3). The resulting PCR fragment was cloned into the pEW plasmid using SalI and BamHI. The 3 UTR of rpiRc was amplified using primers BamHI-EcoRI-3UTR-RpiRc-fw and SmaI-3UTR-RpiRc-rvs (Supplementary Table S3) and inserted downstream of the gfp gene using restriction sites BamHI and SmaI. The 3 UTRs carrying the IS were amplified from p 3xF RpiRc+3 UTR IS256 and p 3xF RpiRc+3 UTR IS1181 using BamHI-EcoRI-3UTR-RpiRc-fw and KpnI-term-RpiRc and introduced into the pEW-GFP plasmid using BamHI and KpnI. Lastly, the pGFP-3 UTR-rpiRc was constructed using the pGFP-3 UTR-rpiRc as a template and primers SalI-GFP-fw and KpnI-term-RpiRc. The amplification product was then ligated into the pEW-GFP using SalI and KpnI.

Total protein extraction and Western blotting
Preinocula were grown in 5 ml of TSBg supplemented with Erm (TSBg+Erm) and incubated ON at 37 • C and 200 rpm. Bacterial concentrations were estimated by measuring the optical density (OD 600 ) of the preinocula and then normalized to an OD 600 of 0.02 in Erlenmeyer flasks containing TSBg+Erm. Cultures were grown at 37 • C and 200 rpm until an OD 600 of 0.5 (exponential phase) was reached for strains carrying IcaR and FtnA constructs and around 5-6 (stationary phase) for those carrying the RpiRc plasmids. Culture samples were harvested by centrifugation (10 min at 4,400 g and 4 • C) and bacterial pellets stored at −20 • C until needed. Pellets were thawed, washed and resuspended in 1 ml of phosphate-buffered saline (PBS). Next, bacterial suspensions were transferred to Fast Prep tubes containing acid-washed 100 m glass beads (Sigma) and cells lysed using a FastPrep-24 instrument (MP Biomedicals) at speed 6 for 45 s twice. Total cell extracts were obtained by centrifuging the Fast Prep tubes for 10 min at 21,000 g and 4 • C. The total protein concentration was quantified using the Bio-Rad protein assay kit. Samples were prepared at the desired concentration in Laemmli buffer and stored at −20 • C until needed.
Western blotting was performed as previously described (36). The 3xFLAG tagged IcaR, RpiRc and FtnA samples were incubated with mouse monoclonal anti-FLAG M2-Peroxidase (HRP) antibodies (Sigma) diluted 1:1,000 whereas the GFP samples were incubated with mouse monoclonal anti-GFP antibodies 1:5,000 (Living Colors, Clontech) and peroxidase-conjugated goat antimouse immunoglobulin G and M antibodies 1:2,500 (Pierce-Thermo Scientific). Membranes were developed using the SuperSignal West Pico Chemiluminiscent Substrate kit (Thermo Scientific). Protein bands were quantified by densitometry of Western blot images using ImageJ (http: //rsbweb.nih.gov/ij/). Each of the protein levels was normalized to either the levels of S. aureus (Sau) or the WT sample.

RNA extraction and Northern blotting
Bacteria were grown as described in the previous section. Cultures were centrifuged for 3 min and 4,400 g at 4 • C and pellets stored at −80 • C. Total RNA extraction was performed as previously described (30). Northern blotting was performed as described by Caballero et al (36) with the following modifications. Radiolabelled-RNA probes were synthesized from PCR products carrying the T7 promoter (Supplementary Table S3) using the MAXIscript T7 transcription kit (Ambion) and 32 P-␣-UTP, following the manufacturer's recommendations. Probes were then purified with Illustra MicroSpin G-50 columns (GE Healthcare) and membranes were hybridized with the corresponding RNA probes at 68 • C ON. mRNA levels were quantified by densitometry of Northern blot autoradiographies using ImageJ (http://rsbweb.nih.gov/ij/). Each of the mRNA levels was normalized either to the levels of S. aureus (Sau) or the WT sample.

RNA stability assay
Bacteria were grown as indicated in the previous sections. Once the exponential phase (OD 600 = 0.5) was reached, six aliquots of 20 ml of the culture were transferred to 50 ml Falcon tubes containing 300 g ml −1 of Rifampicin and incubated at 37 • C for 0, 2, 4, 8, 15 and 30 min. Then, 5 ml of the Stop solution (95% ethanol, 5% phenol) were added to the samples and centrifuged for 2 min at 4,400 g. Pellets were frozen in liquid nitrogen and stored at −80 • C. RNA extraction and Northern blot analysis were conducted as described above.

Identification of elements disrupting the rpiRc mRNA sequence
All available S. aureus genome sequences were retrieved from ftp://ftp.ncbi.nlm.nih.gov/genomes/ through the wget command-line utility. The Fasta files were loaded into a local nucleotide database using the Geneious software package. Next, a Blastn algorithm was locally run using the rpiRc mRNA as a query. Genomes presenting pairwise alignment disruption located at the rpiRc 3 UTR were registered. When possible, sequences surrounding the position of alignment disruptions were used to look for insertion sequences in the IS database (https://isfinder.biotoul.fr/) (40). Note that in some cases, due to the assembly protocol, repeated sequences were eliminated from contig ends avoiding the identification of the IS responsible for the alignment disruption.

Haemolysis assay
Bacterial cultures were grown in 15 ml of BHI ON at 37 • C and 200 rpm. Cultures were centrifuged for 10 min and 4,400 g at 4 • C to collect the supernatants containing the haemolysins. Supernatants were equalized using the wet weight of the pellets and the remaining cells removed from by filtration employing 0.2 m pore filters. Supernatants were concentrated 10 times using a SpeedVac Concentrator. 30 l of the concentrated supernatants were transferred into 5-mm holes within Columbia Sheep blood (5%) agar plates (Biomerieux). Plates were incubated at 37 • C ON and then kept at 4 • C for several days.

In vitro transcription
PCR fragments containing the T7 promoter were used as templates for in vitro transcription. The T7 promoter was included in the corresponding Fw oligonucleotides (Supplementary Table S3). In vitro transcription was carried out using the T7 polymerase at 37 • C ON followed by the removal of the PCR templates with a DNase I treatment. Transcripts were run on a denaturing 6% polyacrylamide gel and visualized using UV shadowing. RNAs were excised from the gels and eluted by incubating the gel fractions with elution buffer (0.5 M ammonium acetate, 1 mM EDTA, 0,1% SDS) and phenol (pH 4.5) at 4 • C ON. RNAs were then purified using a phenol-chloroform extraction and precipitated with ethanol. Pellets were washed with 70% ethanol and resuspended in water. The RNA integrity was checked by polyacrylamide gel electrophoresis and their concentration measured with a NanoDrop Instrument (Agilent Technologies).

-End labelling of synthetic RNAs
Prior to labelling, the RNA was dephosphorylated with FastAP (Thermo Scientific) at 37 • C for 15 min followed by a phenol-chloroform extraction. 7 g of dephosphorylated RNA were incubated with 32 P-␥ -ATP and PNK (Thermo Scientific) at 37 • C for 1 h. Labelled nucleic acids were purified on a denaturing 6% polyacrylamide gel as described in the previous section. The efficiency of the labelling was measured using a scintillation counter.

Electrophoretic mobility shift assays
Labelled and non-labelled RNAs were separately prepared at the desired concentrations by adding sterile MilliQ water and 5× renaturing buffer (100 mM Tris-HCl pH 7.5, 300 mM KCl, 200 mM NH 4 Cl, 15 mM DTT). Samples were denatured at 90 • C for 1 min and chilled on ice for an additional minute. Next, 10 mM of MgCl 2 were added and samples were incubated for 10 min at 22 • C. Subsequently, the labelled RNA was mixed with increasing concentrations of non-labelled samples in 1X reaction buffer (20 mM Tris-HCl pH 7.5, 60 mM KCl, 40 mM NH 4 Cl, 3 mM DTT, 10 mM MgCl 2 ) in a total volume of 20 l. The mixtures were then incubated at 37 • C for 15 min before running them in a non-denaturing 6% polyacrylamide gel and 1X Tris-borate buffer containing 1 mM MgCl 2 at 300 V and 4 • C. The gel was developed by autoradiography.

PIA-PNAG quantification
PIA/PNAG exopolysaccharide was extracted and quantified as previously described (5,41). Cell surface extracts (5 l), which were obtained from bacteria grown ON in TSBg, were spotted onto a nitrocellulose membrane using a Bio-Dot microfiltration apparatus (Bio-Rad). The membranes were blocked ON with 5% skimmed milk in phosphate-buffered saline with 0.1% Tween 20 (PBS-Tw20). After washing several times with PBS-Tw20, the membranes were incubated for 2 h with specific anti-PNAG antibodies diluted 1:20,000 (42). Bound antibodies were detected by incubating the membranes for 1 h with peroxidaseconjugated goat anti-rabbit immunoglobulin G antibodies (Jackson ImmunoResearch Laboratories, Inc., Westgrove, PA) diluted 1:10,000. Finally, membranes were washed with PBS-Tw20 and developed using the SuperSignal West Pico Chemiluminescent Substrate (Thermo Scientific) and Amersham Hyperfilm ECL films (GE Healthcare Life Sciences).

UTRs are highly variable among staphylococcal species
Next-generation RNA sequencing (NGS) technologies allow to accurately determine mRNA boundaries and, therefore, determine the lengths of the untranslated regions (UTRs). To evaluate the evolutionary relationship between coding sequences (CDSs) and their corresponding 3 UTRs, we performed genome-wide comparative analyses. Using high-resolution transcriptome maps of S. aureus as a reference, we analysed the 3 UTR sequence conservation of close phylogenetic members of the genus Staphylococcus (Supplementary Figure S1). In order to simplify the study, we focused on monocistronic mRNAs. This ensured that each analysed CDS was flanked by a 5 and a 3 UTR. First, we generated a query database including the 3 UTR sequences plus the last 200 nucleotides (nt) of each monocistronic CDS (∼590 sequences). This stretch of nucleotides was intended to serve as an indicator for the correct identification of orthologous genes among staphylococcal species. Using blastn (43), we compared the selected sequences against 8 representative genomes of the most closely-related staphylococcal species according to the TimeTree knowledge-base (44) (Supplementary Figure S2). The last conserved nucleotide position of each of the orthologous mRNAs was registered (Supplementary Table S4) and plotted against the length of their corresponding S. aureus mRNA (Figure 1A and Supplementary Figure S3). The comparison of S. aureus NCTC 8325 and MW2 strains revealed that 97% of the analysed mRNAs fell on the diagonal axis of the plot, indicating 3 UTRs conservation. Only 17 (3%) out of 589 monocistronic mRNAs showed 3 UTR sequence variations, which could be explained by the presence of unique insertion sequences (ISs) or Staphylococcus aureus repeat (STAR) elements (45).
However, 3 UTRs conservation among the analysed mR-NAs decreased as the phylogenetic distance between S. aureus and the other staphylococcal species increased. This phenomenon is reflected by the dots being positioned along the horizontal line that represents the stop codon in Figure 1A and Supplementary Figure S3A. The lack of conservation was already noticeable for several 3 UTRs of S. argenteus, which is the closest species to S. aureus. This increased considerably in S. simiae and left the remaining species of the genus with just a few conserved 3 UTRs (Figure 1A and Supplementary Figure S3A). As one would expect, among the conserved 3 UTRs (diagonal line) for all the staphylococcal species, we identified RNAIII and seven putative riboswitch-dependent 3 UTRs, meaning that the mRNA transcription ended at the TT of the downstream riboswitch when in the OFF configuration (5,(23)(24)30). We also found two mRNAs, encoding hypothetical proteins SAOUHSC 02781 and SAOUHSC 02702, which produced the 3 UTR-derived sRNAs RsaOT (i.e. srr43 or SAOUH-SCs080) and Teg130 (i.e. SAOUHSCs100), respectively (46)(47)(48). The remaining mRNAs with conserved 3 UTRs encoded enolase (eno), ribosomal protein L32 (rpmF) and a small membrane protein (SAOUHSC 01055a).
It is worth noting that, in most of the mRNAs, nucleotide conservation was lost downstream of the orthologous CDS ( Figure 1A and Supplementary Figure S3A). This is illustrated in Figure 1B  It is expected that non-coding sequences accumulate more changes throughout the course of evolution than their corresponding CDSs, which cannot significantly change without dramatically affecting protein functionality. In agreement with this, we reasoned that 5 UTRs (excluding the ribosome binding region) should suffer a similar degree of nucleotide variation than the 3 UTRs. Therefore, we performed the same blastn comparison, using a query database that included the 5 UTR sequences plus the first 200 nt of the CDS (Supplementary Table S5). As shown by the scatter plots in Supplementary Figure S4, evident variations among the 5 UTRs of orthologous genes existed. However, the number of conserved 5 UTRs was higher than their 3 UTR counterparts. This was significant when comparing the 5 UTR and 3 UTR scatter plots of a specific staphylococcal species, as the number of dots closer to the diagonal line is always greater for the 5 UTRs ( Supplementary Figure S4). Although this analysis may have been biased by the presence of the ribosome binding site, it supported the idea of 3 UTRs being more prone to evolutionary changes than 5 UTRs.

Most of the mRNAs encoding orthologous proteins contain 3 UTRs with different lengths and sequences in Staphylococcaceae
The fact that the conservation of several 3 UTRs was almost fully lost after the CDSs indicated that the mRNAs encoding orthologous proteins had different 3 UTRs. To determine and compare the real length of 3 UTRs among closely-related staphylococcal species, we performed RNA sequencing of whole transcripts from the S. simiae, S. epidermidis and S. capitis strains. The obtained transcriptomic maps were loaded into a Jbrowse-based server (27), which is publicly available at http://rnamaps.unavarra.es. Then, the lengths of the 3 UTRs from each S. aureus orthologous monocistronic gene were manually annotated by combining the RNA-seq data and the Rho-independent transcriptional terminator predictions (using TransTermHP algorithm) (32). Subsequently, we plotted the real 3 UTR lengths from each staphylococcal species against the 3 UTR lengths of their corresponding S. aureus orthologous genes. Figure 2A shows that most of the 3 UTR lengths from S. au-reus did not correlate with those of the three other analysed staphylococcal species (Figure 2A). This lack of correspondence was exemplified when selecting some biologically relevant genes (Figures 2B-D). As reflected in Figures 2B and D, in most cases, the length of the 3 UTR showed variability among all the analysed species, while a similar 3 UTR length for a given gene was only preserved in a few examples. Overall, these data confirmed that most of the mR-NAs encoding orthologous genes in staphylococcal species had 3 UTRs with different lengths in addition to sequence variation.

UTR sequence differences are originated by local nucleotide changes and gene rearrangements
In order to understand how 3 UTR differences originated among mRNAs encoding orthologous proteins, we compared the whole genomic sequence of S. aureus (Sau) against that of S. simiae (Ssim), S. epidermidis (Sepi) and S. capitis (Scap) using Mauve (49). Note that we only focused on the monocistronic mRNAs, as stated above. We ex- pected 3 UTR sequence variations to occur due to different genomic rearrangements, as previously described for some sRNAs located at intergenic regions (IGRs) (50). Therefore, we checked the sequence conservation in genomic regions downstream of the orthologous CDSs, including IGRs and adjacent CDSs, for each genomic pair (Sau versus Ssim, Sau versus Sepi and Sau versus Scap). We assigned a value of 1 or 0 depending on whether the downstream IGR/CDS was conserved or not, respectively (Figure 3A). Confirming previous blastn analysis (Figure 1), we found that only ∼6.8-10.6% of the S. aureus CDSs had a conserved downstream IGR/CDS when compared to the abovementioned species ( Figure 3B). On the other hand, 28-35% of the orthologous CDSs lacked a conserved downstream IGR/CDS, indicating that gene rearrangements in these loci led to 3 UTR variations. The rest of CDSs showed a non-conserved downstream IGR but a preserved downstream CDS (between 58% and 62% depending on the species analysed) ( Figure 3B). When performing this analysis, we realized that such variations substantially altered the IGR lengths ( Figure 3C). Considering that the 3 end conservation was lost around the protein stop codon ( Figure  1), IGR length variations may be attributed to shifts in the position of the stop codon, which may ultimately affect the CDS length. To address such a possibility, we compared the lengths of S. aureus CDSs to those of their corresponding orthologous genes for the three abovementioned staphylococcal species. As Figure 3D shows, most of the orthologous CDSs had similar lengths, indicating that variations in the lengths of the IGRs were not due to differences in the analysed CDS lengths ( Figure 3D). Therefore, we concluded that the differences on the 3 UTR sequences occurred partly due to gene rearrangements and mostly because of sequence variations in the IGR sequences. The causes for local variations were unknown.

Species-specific 3 UTRs variations affect the expression of orthologous genes
3 UTR variations would be physiologically relevant if they led to differences in the expression of the orthologous genes of a particular bacterial species. To study this in detail, we selected the icaR (5), ftnA and rpiRc genes as examples since we had previously found that their expression was affected by the deletion of their corresponding 3 UTR. In addition, according to the single nucleotide polymorphism (SNP) analysis carried out by Joseph and colleagues (51), the 3 UTRs sequences of these mRNAs are highly conserved among different S. aureus strains (Supplementary Figure S5). At the same time, these 3 UTRs present different levels of conservation when compared to their closest phylogenetically-related staphylococcal species ( Figure  4A). In addition, transcript mapping by RNA-seq revealed variations in the 3 UTR lengths ( Figure 2B and D). According to synteny analysis, the icaR 3 UTR showed significant variations in length and sequence due to gene rearrangements (Supplementary Figure S6A), while the differences in the ftnA and rpiRc 3 UTRs were most likely caused by local nucleotide variations in the IGRs, as no gene rearrangements were found by synteny analysis (Supplementary Figure S6B and C). Although a CDS insertion was found downstream of the ftnA in S. argenteus, it did not affect the mRNA sequence (Supplementary Figure S6B).
To confirm the transcript boundaries of the icaR, ftnA and rpiRc mRNAs from S. aureus, S. simiae, S. epidermidis and S. capitis, we performed a simultaneous mapping of their 5 and 3 mRNA ends using a modified version of the rapid amplification of cDNA ends technique (mRACE). Supplementary Figures S7 to S9 show the mRACE results that confirmed the mRNA mappings from the RNA-seq data in combination with the transcriptional terminator predictions. Note that rpiRc in S. epidermidis and S. capitis is the last gene of a polycistronic transcript. Therefore, we were unable to obtain mRACE mapping results from the native gene. We used Mfold to predict the putative transcriptional terminator structures of the transcripts (52). Interestingly, these stem-loops were different for the majority of the cases (Supplementary Figures S7B, S8B and S9B).
Based on the information collected, we constructed chimeric mRNAs that combined the CDSs of the three selected S. aureus genes with the 3 UTR of their corresponding orthologous genes (from the same staphylococcal species analysed in Figure 4A). We fused such 3 UTRs to their corresponding S. aureus CDS, which carried the 3xFLAG sequence at the N-terminal, to monitor their protein expression ( Figure 4B). These plasmids were then electroporated into the mutant strains that lacked the original genes to facilitate the detection of the mRNAs that were expressed from the plasmids. To confirm the mapping of the 3 end of the chimeric rpiRc mRNA carrying the 3 UTRs from S. epidermidis and S. capitis, we performed mRACE. The results showed that the transcripts ended at the expected sites, which were downstream of the predicted stem-loops (Supplementary Figure S9).
Expression analysis of the 3xFLAG-tagged proteins by Western blotting revealed that the deletion of the 3 UTRs produced an increase of the IcaR, Ftn and RpiRc proteins. Interestingly, the chimeric icaR mRNAs harbouring the UCCCC motif (3 UTRs from S. simiae and S. argenteus) ( Figure 4A), which was required for modulation of IcaR translation in S. aureus (5), presented similar IcaR protein levels to those shown by the S. aureus wild-type (WT) mRNA. In contrast, the chimeric icaR mRNAs lacking the UCCCC motif (S. epidermidis and S. capitis icaR 3 UTRs) expressed a comparable amount of the IcaR protein to that of a S. aureus mRNA carrying the 3 UTR deletion ( Figure  4C). Chimeric mRNAs expressing FtnA and RpiRc showed protein yields similar to those of the same mRNAs lacking the 3 UTR. The exception to such behaviour was the chimeric construct carrying the ftnA mRNA in combination with the 3 UTR of S. argenteus, whose sequence was conserved when compared to S. aureus ( Figure 4A). The same samples were then subjected to Northern blot analyses. The results showed that, in most cases, the steady-state mRNA levels correlated with the protein levels ( Figure 4D). This suggested that the increase in protein levels, which occurred when the original 3 UTR sequence was substituted by the one present in its orthologous mRNA, was most likely due to an increase in mRNA stability, as previously described for the icaR 3 UTR mRNA (5). To confirm this, we measured the half-life of the chimeric icaR+3 UTR Sepi mRNA as an example. Supplementary Figure S10 shows that this chimeric mRNA presented a higher half-life in comparison with the one from the icaR WT mRNA (Supplementary Figure S10). Note that the icaR-3 UTR Scap and ftnA-3 UTR Sepi chimeras did not show increased mRNA levels, implying that protein expression could be affected in other ways. Altogether, these results indicated that species-specific variations in the 3 UTR sequences could not reproduce the expression levels found in S. aureus, suggesting putative distinct functional roles for 3 UTRs in the different species.

IS1181 transposition to rpiRc 3 UTR modifies RpiRc expression and haemolysis production
The rpiRc 3 UTR in the S. simiae strain is considerably longer than the 3 UTRs from other staphylococcal species ( Figure 4A). A deep sequence analysis revealed that the S. simiae rpiRc 3 UTR contains three copies of STAR elements (45), which are distributed among staphylococcal species (53). Since the number and genomic locations of the STAR element copies are variable among species and strains (53), we investigated whether insertions of STAR elements may also occur in the rpiRc 3 UTR of S. aureus and, therefore, alter protein expression as observed in the chimeric rpiRc construct carrying the 3 UTR of S. simiae ( Figure 4C). Since SNP analysis were not indicative of the presence of insertion/deletion sequences, we performed rpiRc mRNA sequence pairwise alignments to look for alignment disruptions within the rpiRc 3 UTR. We reasoned that alignment disruptions could help us identify repeated regions. Following this idea, we used the whole rpiRc mRNA sequence as the query to perform blastn against all draft and complete S. aureus genomes available on the NCBI database. The results showed that 33 S. aureus genomes presented alignment disruptions within the rpiRc 3 UTR (Supplementary Table S6). Since many of these genomes were not fully closed, some alignment disruptions occurred due to the query sequence aligning to different contigs. By analysing the extremes of the aligned contigs, we identified that 16 out of the 33 genomes contained IS insertions in different locations of the rpiRc 3 UTR sequence. Specifically, we were able to map nine IS1181 (54), five IS256 (55) and two ISs of the IS30 transposase family ( Figure 5A and Supplementary Table S6). The remaining genomes presented sequence deletions or duplications. Interestingly, we identified duplications of 8-12 nt in some of the strains, which were indicative of the presence/excision of ISs. The fact that no IS sequences were found in the related contigs may be explained by either the use of assembly algorithms that eliminated repeated sequences from the extremes of the contigs or the ISs being excised. However, the latter would not explain why these contigs were not assembled together, suggesting that an unknown IS may be inserted in the rpiRc 3 UTR of these genomes.
In order to analyse the consequences of such insertions, we constructed two plasmids containing the IS256 and IS1181 sequences 42 and 108 nt downstream of the RpiRc stop codon. These mimicked the chromosomic configurations found in the rpiRc region of the S. aureus 2010-60-6511-5 and DAR1183 strains, respectively ( Figure 5B). In both cases, the insertion orientation of the IS placed the transposase gene in the DNA strand opposite to that of rpiRc. Northern blot analysis revealed different chimeric mRNA lengths ( Figure 5C). On the one hand, IS1181 carried a stem-loop downstream of the transposase gene, in the same DNA strand as the rpiRc gene, that acted as a transcriptional terminator ( Figure 5B). Therefore, the IS1181 insertion produced a chimeric rpiRc-IS1181 mRNA of ∼1.4 kb (3 UTR-IS1181 of 222 nt) that included a few nucleotides of the IS1181 sequence. This band migrated in a similar manner as the original rpiRc mRNA. The impact of this insertion was an increase in the mRNA levels (Figure 5B and C), which correlated with an increase in the 3xF RpiRc protein levels ( Figure 5D). On the other hand, since the IS256 lacked stem-loops able to prevent transcription from adjacent genes, a chimeric rpiRc-IS256 mRNA of ∼2.8 kb (3 UTR-IS256 ∼1.5 kb) was produced. In contrast to the IS1181 insertion, the chimeric rpiRc-IS256 mRNA included the whole IS256 sequence ( Figure 5B and C). In this case, there were no considerable changes neither in the mRNA nor 3xF RpiRc protein levels ( Figure 5C and D). Overall, these examples indicated that IS insertions may produce different RpiRc protein expression outcomes de-pending on the type of IS, the insertion site and the orientation of the IS.
To further confirm that the rpiRc 3 UTR had a functional role that could be affected by IS insertions, we analysed the effect of the WT 3 UTR and the chimeric 3 UTR-ISs on the expression of a heterologous gene, such as the green fluorescent protein (gfp). To this end, we fused the rpiRc 3 UTR WT, the 3 UTR-IS1181 and 3 UTR-IS256 sequences downstream of the gfp gene. As a control, we included the 3 UTR plasmid that carried the rpiRc TT downstream of gfp ( Figure 5E). Western blot analyses showed that the rpiRc 3 UTR lead to a lower GFP expression when compared to the 3 UTR version ( Figure 5E). This indicated that the rpiRc 3 UTR could reduce the expression of a heterologous gene and, therefore, act as a functional module. However, the presence of IS1181 in the rpiRc 3 UTR increased GFP expression while IS256 only slightly decreased it. This confirmed that the transposition of ISs in the rpiRc 3 UTR could have different consequences on protein expression depending on the insertion events ( Figure 5E).
Next, we decided to analyse the relevance of the rpiRc 3 UTR on the S. aureus physiology and the implications of the IS insertions in its role. Previous studies showed that RpiRc repressed the expression of RNAIII. RNAIII is the master regulator of S. aureus virulence, which activates many virulence factors, including haemolysins. A mutant lacking the rpiRc gene exhibited an increased haemolysis capacity due to the lack of RNAIII repression (56)(57)(58). To test whether alterations in the rpiRc 3 UTR could modify the haemolytic activity of S. aureus, we constructed chromosomal mutants in S. aureus MW2, a communityacquired agr-positive methicillin-resistant (MRSA) strain. On the one hand, we performed a 115-bp chromosomal deletion that removed most of the 3 UTR of rpiRc while preserving its intrinsic TT (MW2 3 UTR), as we did for the aforementioned construct. On the other hand, we inserted IS1181 in the same rpiRc region as it is naturally found in S. aureus DAR1183 (MW2 3 UTR+IS1181). As a control, we deleted the whole rpiRc gene (MW2 rpiRc). These strains were grown overnight and the production of haemolysins in the culture supernatants was tested in sheepblood agar plates. As expected, the rpiRc mutant showed an increase in haemolysis when compared to the WT strain ( Figure 5F). In contrast, the deletion of the rpiRc 3 UTR and the IS1181 insertion in the 3 UTR, which produced higher protein levels of the RpiRc repressor ( Figure 5D), resulted in lower levels of haemolysis ( Figure 5F). Overall, these results highlighted the relevance of the rpiRc 3 UTR in the physiology of S. aureus while demonstrating that the RpiRc function could be altered by the insertion of IS1181 in the rpiRc 3 UTR.

Species-specific icaR 3 UTRs produce different PNAG levels in S. aureus
The expression analysis of chimeric icaR mRNAs indicated that the icaR 3 UTRs from S. epidermidis and S. capitis were unable to reproduce the 3 UTR-mediated regulation observed in S. aureus ( Figure 4C). We previously showed that the icaR 3 UTR interacted with the 5 UTR of the same mRNA molecule through a UCCCC motif to repress IcaR chimeric mRNAs that are generated when IS1181 and IS256 are inserted in the rpiRc 3 UTR of the S. aureus DAR1183 and 2010-60-6511-5 strains, respectively. Since IS1181 harbours a stem loop that serves as a transcriptional terminator (TT) in the same strand as the rpiRc gene, which is opposite to the transposase gene, the insertion generates a chimeric 3 UTR that is 222 nt-long. In contrast, the IS256 lacks TTs, generating a chimeric mRNA that includes the whole IS256 sequence. (C) Northern blots showing the mRNA levels expressed from the constructs harbouring the WT rpiRc gene and the IS insertions that mimic the configurations found in the S. aureus DAR1183 and 2010-60-6511-5 strains. The plasmid expressing the rpiRc mRNA 3 UTR was included as a control. The mRNA levels of rpiRc were detected with a CDS-specific antisense radiolabelled RNA probe. Ribosomal RNAs stained with Midori Green are shown as loading controls. (D) Western blot showing the 3xF RpiRc protein levels expressed from the constructions used in C. The Western blot was developed using peroxidase conjugated anti-FLAG antibodies. A Coomassie stained gel portion is shown as a loading control. The Western and Northern blot images show the representative results from at least three independent experiments. The protein and mRNA levels were quantified by densitometry of Western blot and Northern blot autoradiographies using ImageJ (http://rsbweb.nih.gov/ij/). Each of the protein or mRNA levels was normalized to the levels of the WT rpiRc mRNA. (E) Western blot showing the GFP protein levels of the constructs harbouring the gfp fused to the 3 UTR, 3 UTR+IS256, 3 UTR+IS1181 of rpiRc and the 3 UTR as a control. Western blot was developed using monoclonal anti-GFP antibodies and quantified as described in D. (F) Haemolytic halos produced by the haemolysins contained in the supernatant of the cultures from S. aureus MW2 wild-type (WT) and their isogenic chromosomic mutants: 3 UTR, rpiRc and 3 UTR+IS1181. The supernatants were concentrated 10 times and loaded into 5-mm holes in Columbia Sheep blood (5%) agar plates. translation and, in return, modulate PIA-PNAG biosynthesis and biofilm formation in S. aureus (5). The absence of the UCCCC motif in the S. epidermidis and S. capitis icaR 3 UTRs may account for the differences in IcaR expression. This may be attributed to the lack of interaction between their 5 and 3 UTRs. To test this hypothesis, we carried out electrophoretic mobility shift assays (EMSAs). The S. aureus icaR 5 UTR was synthesized in vitro, radioactively labelled and then incubated with increasing concentrations of cold icaR 3 UTRs from the five staphylococcal species analysed in Figure 4. Figure 6A shows that the S. aureus icaR 3 UTR produced the expected gel-shifts when interacting with the 5 UTR, as previously described (5). Consistently, a similar behaviour was observed when using the 3 UTR of S. simiae ( Figure 6B). In the case of the S. argenteus icaR 3 UTR, despite sharing a relatively high sequence identity with its S. aureus counterpart (82% of identity), a much lower binding affinity was found ( Figure 6B). This indicated that although S. argenteus preserved the UC-CCC motif, the nucleotide differences might affect the binding affinity in vitro. S. epidermidis and S. capitis 3 UTRs, as one would expect, did not produce gel shifts when incubated with the 5 UTR due to the absence of the UC-CCC motif in their sequences ( Figure 6C). Moreover, recent results showed that the 5 and 3 UTRs from S. epidermidis, unlike what was observed for S. aureus, did not interact between them, excluding the existence of an orthologous mechanism (59). In contrast, it was found that the S. epidermidis-specific sRNA IcaZ, downstream of the icaR 3 UTR, targeted the icaR 5 UTR and prevented IcaR translation (59). Altogether, these results indicated that sequence variations around the 3 UTR determined the presence or absence of specific motifs and, ultimately, the generation of functional distinctions.
In order to demonstrate that the differences observed in terms of IcaR expression by chimeric mRNAs (Figure 4D) were enough to alter PIA-PNAG biosynthesis and, thus, be of biological relevance, we introduced the plasmids expressing these chimeric mRNAs into the S. aureus 15981 strain, which produced high levels of PIA-PNAG exopolysaccharide. The quantification of PIA-PNAG levels by dot-blot using specific antibodies against the exopolysaccharide showed that the strains expressing higher levels of IcaR (the chimeric Scap/Sepi and 3 UTRs strains, which lacked the UCCCC motif) presented lower levels of PIA-PNAG in comparison with the strains producing IcaR from mRNAs containing the UCCCC motif (the WT and chimeric Sarg/Ssim strains) ( Figure 6D). These results confirmed that evolutionary variations in 3 UTRs, which occur among phylogenetically-related bacteria, are enough to create gene expression differences that affect important biological processes, such as exopolysaccharide production. This suggested that the synthesis of orthologous proteins might be differentially affected by 3 UTR sequence variations.

UTR variability among mRNAs encoding orthologous proteins is a common trait in bacteria
Having proved that the inter-species 3 UTR variability had consequences on S. aureus gene expression, we decided to extend the 3 UTR conservation analysis to other bac-terial genera. First, we focused on several previously described functional 3 UTRs (4,(6)(7)(8). We compared the last 200 nt of the CDSs plus the entire 3 UTR sequences from the E. coli acnB, Y. pestis hmsT, C. glutamicum aceA and B. subtilis hbs mRNAs with those of their corresponding phylogenetically-related species ( Supplementary Figures S11A-S14A). Although the CDSs were highly conserved across all closely-related species, there were 3 UTR sequence variations in all four gene examples (Figure 7). The E. coli acnB 3 UTR was conserved in the Citrobacter, Enterobacter and Klebsiella species with only a few nucleotide variations. In contrast, and according to synteny analysis by SyntTax (60), the Salmonella, Cronobacter and Serratia acnB 3 UTRs sequences varied due to gene rearrangements (Supplementary Figure S11B). For example, in S. enterica serovar Typhimurium SL1344, an insertion of the SL1344 0160 gene (which encodes a putative HNH restriction endonuclease) deleted the sequence downstream of the acnB CDS and included a different TT. As a result, the Salmonella acnB 3 UTR was only 34 nt long and it lacked the stem-loop recognized by the apo-AcnB protein (7).
Sequence variations were also present in the C. glutamicum aceA, B. subtilis hbs and Yersinia hmsT 3 UTRs. Although the hbs and hmsT CDSs are widely distributed across the Bacillus and Yersinia genera, their 3 UTRs were only conserved in the closest species: B. amyloliquefaciens, B. licheniformis and B. pumilus in the former and Y. pestis biovars and Y. pseudotuberculosis in the latter (Figure 7). Finally, none of the 10 compared species of the genus Corynebacterium that carried the aceA CDS presented a similar sequence to that of the C. glutamicum aceA 3 UTR (Figure 7). We carried out the prediction of Rho-independent transcriptional terminator structures in the non-conserved 3 UTRs to determine the length variability (Figure 7 and Supplementary Figures S11-S14). Most of these 3 UTRs variations were explained by the gene rearrangements occurring downstream of CDSs, with the exception of the hbs 3 UTR, in which variations may have occurred locally (Supplementary Figures S11-S14).
Finally, to further investigate whether 3 UTR variations could also be distributed among other genes of E. coli and B. subtilis, as observed in S. aureus, we performed similar genome-wide 3 UTR conservation analyses. Since the 3end transcript boundaries were recently released for E. coli and B. subtilis (28,29), we generated a query database including the sequences of all mapped 3 UTRs plus the last 200 nt of their corresponding CDSs. Using blastn, we compared the E. coli and B. subtilis databases against representative genomes of species from the Enterobacteriaceae family and the genus Bacillus, respectively. We found a similar outcome to the one observed for the genus Staphylococcus, meaning that 3 UTR sequences were also variable among mRNAs encoding orthologous proteins (Figure 8, Supplementary Figure S15, Supplementary Tables S7 and  S8). The bacterial species of the Enterobacteriaceae family showed a higher number of conserved 3 UTRs in comparison with the Bacillus species. This was probably due to the phylogenetic distances between the species in the former being shorter than the ones in the latter. However, both bacterial groups presented a very high variation rate of their 3 UTR sequences, as illustrated in Figure 8  mentary Figure S15, where most of the dots fell onto the horizontal line, which represented the protein stop codon. This indicated that nucleotide conservation was also lost downstream of most CDSs (Supplementary Figure S16). On the one hand, these comparative analyses confirmed that the nucleotide conservation analysis might only be useful for detecting a small percentage of inter-species conserved 3 UTRs. On the other hand, the results indicated that the 3 UTRs from orthologous genes underwent different evolutionary events, a phenomenon that seems to be widespread in bacteria.

DISCUSSION
Bacterial IGRs, which have been traditionally defined as the genomic regions between two CDSs, are known for being poorly conserved among phylogenetically-related bacteria (61)(62)(63)(64). For this reason, they have often been regarded as 'junk' genetic material lacking relevant biological functions. However, recent RNA-seq technologies have revealed that IGRs may be more complex than initially anticipated, accommodating a plethora of regulatory elements that are either transcribed independently as sRNAs or fully functional from within the UTRs of mRNAs. In recent years, several studies have shown the biological relevance of 3 UTRs both in bacteria and eukaryotic organisms (for recent reviews see (18)(19)65)). Since the 3 UTR-regulatory elements are physically-and functionally-related to the adjacent CDSs through mRNAs, an evolutionary connection between them might be expected in phylogenetically-related species. This indicates that the 3 UTR sequences would be preserved to the same degree as their corresponding CDSs. However, in this study, we observed this to be true for just a few examples, since most of the 3 UTRs from mRNAs encoding orthologous proteins, in closely-related bacterial species, display different lengths and sequences. The idea of synthetizing hundreds of nucleotides at the 3 ends with-out a relevant function seems contradictory as bacteria tend to efficiently manage energy and clones with a lower fitness are rapidly eliminated from populations. Therefore, an explanation for this relatively high 3 UTR variation rate between related species should be sought. First, it is important to highlight that 3 UTR sequences are more highly conserved in comparison with non-transcribed regions of strains from the same bacterial species (5). This implies that an evolutionary pressure to preserve 3 UTR sequences in a particular species exists. In addition, our data indicates that specific inter-species variations in 3 UTR sequences often occur. Therefore, these variations could differentially affect the expression of orthologous proteins from different species. This could be a strategy used by bacteria to create diversity without changing the CDS, an idea supported by previous studies that show how bacterial genes can shift rapidly between multiple regulatory nodes. Oren et al demonstrated the existence of promoter-regulatory regions that changed among orthologous genes, contributing to expression divergence and conferring distinct fitness advantages (66). Analogously, the finding that 3 UTR sequence variations can affect protein expression depending on the species may be an evolutionarily-selected trait to create species-specific post-transcriptional regulatory processes (16,18,67). Since the 5 UTRs are closely related with the initiation of translation, any significant sequence variations could impair protein synthesis and, thus, explain the bias in variability within 3 UTRs. This hypothesis is in agreement with the role of 3 UTRs in eukaryotes, which have been specifically targeted through the course of evolution to contribute to the divergence of species by accumulating regulatory elements in their sequences (65,(68)(69).
It is worth noting that longer 3 UTRs (including multiple miRNA target sites) are preferred in functionally critical eukaryotic genes that are spatially or temporally expressed. In contrast, housekeeping genes are selected to have shorter 3 UTRs (70) and S. aureus icaR mRNAs (which encode the main transcriptional regulators of the SP-1 pathogenicity island and biofilm formation, respectively) are among the longest bona fide 3 UTRs found in bacteria. Resembling the eukaryotic process in which mRNAs are targeted by miRNAs, it has been recently shown that hilD and icaR 3 UTRs are targeted by sRNAs (71,72). In the former, the Spot 42 sRNA (transcriptionally repressed by CRP-cAMP) targets the hilD 3 UTR, exerting a positive effect on HilD expression (71,73). In the latter, the RsaI sRNA modulates the PIA-PNAG synthesis by interacting with the 3 UTR of the icaR mRNA (72). The RsaI interacting region is located between the UCCCC motif, which is required for the icaR 3 -5 UTR interaction, and the transcriptional terminator of the icaR mRNA (72). Interestingly, while the UCCCC motif is conserved across the S. aureus, S. argenteus and S. simiae icaR mRNAs, the RsaI paring region is only present in the S. aureus icaR mRNA. Note that the RsaI sRNA is conserved in most species of the genus Staphylococcus. The Spot 42 sRNA is present in several bacterial genera, while the hilD gene is exclusive to Salmonella (74). Overall, this shows that the 3 UTRs can accumulate more than one regulatory motif, which may be differentially selected during evolution.
The acquisition of icaR 3 UTR variations may be explained by gene rearrangements that created mRNA chimeras among staphylococcal species. The S. aureus icaR CDS is only present in 5 out of 9 staphylococcal species that encode the icaADBC operon. It is interesting to note that in the four species lacking the icaR CDS either DspB (hexosaminidase), a TetR-like regulator, or proteins of a twocomponent system are encoded in place of the icaR gene. These different genomic organizations indicate that the regulatory locus evolved independently from the icaADBC operon and may also explain the appearance of 3 UTR variations in the icaR mRNA when such a locus was acquired or mobilized. Similarly, gene rearrangements could also be responsible for the disparity of sRNA content in phylogenetically-related bacterial species and may also justify some of the 3 UTR variations (50).
In addition, there may be other phenomena beyond gene rearrangements behind the 3 UTR variations, since only a minor percentage of them (∼28-35%) fall under this category (Figure 3). It is likely that mobile insertion sequences and/or repeated sequences, such as those found in the IGRs of S. aureus, change the functionality of the 3 UTRs of orthologous mRNAs (53,(75)(76). Interestingly, we showed that the presence of STAR elements and ISs that interrupt 3 UTR sequences affect the expression of RpiRc in S. aureus ( Figure 5). This suggests that alternative IS transposition/excision processes may generate local 3 UTR variations that can be evolutionarily selected if they produce an advantage for a specific clone in comparison to the rest of the bacterial population. Alternatively, ISs might be also responsible for phase variation mechanisms (77,78). Here, we showed that the presence of IS1181 in the rpiRc 3 UTR of several S. aureus strains reduced haemolysin production by increasing the levels of the RNAIII-repressor RpiRc ( Figure 5). Further studies would be required to determine whether the transposition of ISs into the rpiRc 3 UTR provide an additional way to regulate S. aureus virulence through the agr locus, as seen for the cases of agrphase-variation (79).
Some of the genomic repeats include stem-loops, which serve as transcriptional terminators ( Figure 5) or provide specific targets for RNase III (80). Moreover, changes in the 3 UTR sequence could also affect mRNA turnover by modifying the 3 UTR accessibility to exoribonucleases. This indicates that the accumulation/elimination of ribonuclease target sites could be an alternative way to specifically modulate the expression of a protein at the post-transcriptional level in a particular species. In agreement with this, the deletion of 3 UTRs (e.g. hilD, aceA, hmsT and icaR), and with them, the putative ribonuclease target sites, increases protein expression (4)(5)(6)9,21). Similarly, the presence of specific target sites for RNA-binding proteins in 3 UTRs could also contribute to regulatory differences among bacterial species (7,11).
In addition to a differential effect on orthologous genes, it should be also considered that variations in 3 UTR sequences may affect the expression of downstream genes. For example, if variations in the 3 UTR sequence alter the strength of the corresponding TT or even eliminate it (e.g. Figure 4 shows that transcriptional read-through occurs in ftnA of some staphylococcal species), transcription would continue to increase the expression of downstream genes or create antisense RNAs if the downstream gene is encoded in the opposite direction (81). Moreover, if a downstream gene is close enough, variations in the 3 UTR sequence could also alter promoter regulatory elements of the next transcription unit, thereby affecting their expression.
Further studies are needed to improve our knowledge of how sequence variations in 3 UTRs may specifically affect gene expression in different bacterial species. The devel-Nucleic Acids Research, 2020, Vol. 48, No. 5 2561 opment of RNA sequencing techniques, mapping the precise 3 -ends of transcripts, as well as ribonuclease cleavage regions and RNA-protein binding sites will considerably contribute to this (28,(82)(83)(84). Combining such information with studies like the one presented here will help improving our understanding of the evolutionary features responsible for the generation of species-specific sequences that may, ultimately, lead to bacterial species differentiation.

DATA AVAILABILITY
The transcriptomic maps of S. simiae CCM 7213T, S. capitis SK14 and S. epidermidis RP62A strains are available at http://rnamaps.unavarra.es/. The transcriptomic data has been also deposited in GEO under accession GSE143600.
Other data or material that support the findings of this study can be made available by the corresponding author upon request.