E. coli strains and plasmids used in this study are listed in Table S1. DH5α and BL21(DE3)pLysS were utilized for cloning and overexpression of Hfq proteins, respectively. Transconjugant E. coli J53 harboring pNDM-HK was a laboratory stock from PL Ho's (Ho et al., 2011). Wild-type E. coli strain MG1655 was adopted for assays and phenotypic studies. Bacteria were grown in LB broth at 37°C under shaking at 250 rpm to the phases indicated. Antibiotic concentrations in growth media were applied as below: ampicillin 100 μg/ml, kanamycin 20 μg/ml, or chloramphenicol 25 μg/ml.

Plasmid preparation, DNA purification, restriction endonuclease cleavage, ligation, and transformation followed protocols of kits or standard methods. The in-frame knockout of hfq in MG1655 followed the methods using IPTG-induced recombinase from pKM208 and electroporation (Murphy and Campellone, 2003). The transcription unit (TU) of NDM-sR3 (sR3) was amplified from the plasmid DNA of pNDM-HK employing primers XhoI-sR3-F and XhoI-sR3-R, and inserted into the Xho I site of pTL01, a derivative of pACYC184 carrying an additional Xho I restriction site, generating pTL02.

E. coli cell pellets were re-suspended in extraction buffer (10 mM Tris pH 8.0 and 1 mM EDTA) and incubated with 20 mg/ml lysozyme (Sigma) for 5 min at room temperature. The mixtures were then mixed in three volumes of TRIzol reagent (Invitrogen) and RNA was extracted by adding one volume of chloroform followed by centrifugation. Total RNA was precipitated in isopropanol and its quality and quantity were determined with a NanoDrop ND-1000 spectrophotometer (Thermo) and TAE agarose gel electrophoresis.

Small RNA (sRNA) was separated and enriched from total RNA by utilizing the mirVana™ miRNA Isolation Kit (Life Technologies) and subjected to the MICROBExpress Kit (Ambion) and Ribo-Zero rRNA Removal Kit (Epicenter) to eliminate rRNA according to the manufacturer's instructions. The concentration of sRNAs was also confirmed by ND-1000 (Thermo), and its quality and integrity were then monitored by Bioanalyzer (Agilent) using RNA 6000 Pico Kit (Life Technologies).

The rRNA-depleted RNA was used to construct the library with the Ion Total RNA-Seq Kit v2 (Ambion) according to the manufacturer's protocol. Libraries were next sequenced using the Ion Torrent Sequencing platform on Ion 316 Chips (Life Technologies). Reads were mapped to reference genome E. coli str. K-12 substr. MG1655 (GenBank accession NC_000913.3) and plasmid pNDM-HK (GenBank accession: NC_019063.1) using TMAP (Smith and Waterman, 1981; Ning et al., 2001; Li and Durbin, 2009, 2010; Li, 2012) and only reads with high mapping quality were kept for downstream analysis. Mapping quality was defined as the rate of uniquely mapped reads. The unique mapping rate of E. coli J53 and J53 carrying pNDM-HK were 88 and 98%, respectively. The mapped sequencing reads were visualized by Integrated Genome Viewer (IGV) 2.3.34 (Robinson et al., 2011). Small RNAs were searched in antisense and intergenic regions based on read-mapping patterns. The initial and terminal bases of sRNAs were determined by choosing nucleotides that had more than 10% coverage of its maximum reads. The sequencing data was deposited in the NCBI Sequence Read Archive (SRA) database (accession number: SRR6703077).

5′-Rapid amplification of cDNA ends (5′RACE) was performed with FirstChoice RLM-RACE Kit (Life Technologies) according to the manufacturer's instructions. Briefly, total RNA was reverse transcribed into cDNA by M-MLV Reverse Transcriptase. PCR was first performed with 5′ RACE Outer Primer (Life Technologies) and a gene-specific outer primer. A nested PCR was next performed to enhance the specificity with 5′ RACE Inner Primer (Life Technologies) and gene-specific inner primer. The round two PCR products were excised and ligated with pGEM-T Easy vector (Promega) and transformed into DH5α. Single colonies were selected and sequenced with M13 forward primer.

Ten μg total RNA was separated on a 6% polyacrylamide gel containing 8 M urea and transferred onto a Hybond-N nylon membrane (GE Healthcare) at 150 mA for 1 h. RNA was cross-linked with the membrane under UV for 2 min. The blot was then pre-hybridized with 5 ml Ultrasensitive Hybridization buffer (Ambion) at 42°C for 1 h. Oligonucleotide probes 5′-end labeled with [γ-³²P]-ATP (PerkinElmer) using T4 polynucleotide kinase (NEB) were added for blot hybridization at 42°C overnight. For NDM-sR6, the ³²P-labeled RNA probe was synthesized by in vitro transcription using MEGAscript T7 Transcription Kit (Ambion) with addition of [α-³²P]-UTP (PerkinElmer). The blot was washed twice with 20 ml SSC buffer for 10 min, and an image was obtained by phosphor-imager. E. coli 5S rRNA was probed as the loading control.

E. coli MG1655 and MG1655Δhfq harboring pTL02 were grown to exponential phase OD₆₀₀~0.6. Cells were treated with 250 μg/ml rifampicin to terminate RNA biosynthesis. Cultures were, respectively, harvested at time points 0, 1, 2, 5, and 10 min after addition of rifampicin. Total RNA of these samples was extracted and NDM-sR3 level was examined by Northern blot. The RNA levels over time were calculated as the percentages to that at time point zero. The in vivo degradation curve was calculated following the fitted equation N(t) = N(0) ∗ 2^−t/t1/2, where N(t) is the RNA level at specified time point, t, N(0) is the initial RNA level, and t_1/2 is the half-life.

The open-reading frame of E. coli hfq was amplified from MG1655 and cloned into pET28a to generate pET28a-hfq with a C-terminal 6 × His tag, and this plasmid was transformed into E. coli BL21(DE3)pLysS for protein purification. Cells containing pET28a-hfq were grown to OD₆₀₀ ~0.6 at 37°C, and Hfq protein was induced by adding 2 mM IPTG. After shaking for 4 h, cells were harvested by centrifugation and resuspended in lysis buffer (20 mM Tris-HCl, 150 mM NaCl, 0.5% Triton X-100, 0.4 mM PMSF, and 2 mM ß-mercaptoethanol, pH 7.5). Then, bacteria were lysed by sonication and precipitate was removed by centrifugation. Supernatant was passed over TALON Cobalt column (Clontech), followed by washing with five bed volumes of high salt buffer (20 mM Tris-HCl, 500 mM NaCl, 20 mM imidazole, 2 mM ß-mercaptoethanol, pH 7.5) and low salt buffer (20 mM Tris-HCl, 150 mM NaCl, 20 mM imidazole and 2 mM ß-mercaptoethanol, pH 7.5). Hfq protein was eluted by elution buffer (20 mM Tris-HCl, 150 mM NaCl, 250 mM imidazole, and 2 mM ß-mercaptoethanol, pH 7.5). Purified Hfq was dialyzed against storage buffer (20 mM Tris-HCl, 200 mM NaCl and 10% glycerol, pH 7.5) at 4°C overnight. Hfq purity and concentration were determined by A₂₈₀/A₂₆₀ and SDS-PAGE (purity > 95%).

[α-³²P]-UTP labeled NDM-sR3 RNA (0.4 pmol) and indicated amounts of purified His-Hfq₆ were mixed in 20 μl of binding buffer (10 mM Tris-HCl, 0.1 mM EDTA, 10 mM NH₄Cl, 10 mM NaCl, 10 mM KCl, pH 7.5) containing 1 U RNasin (Promega). The mixtures were incubated for 30 min at 37°C and stopped by adding 4 μl loading buffer (25% glycerol, 0.05% bromophenol blue). The samples were analysed on an 8% native polyacrylamide gel containing 5% glycerol. The gel was subjected to autoradiography after electrophoresis.

The total RNA for qRT-PCR was treated with 2U TURBO DNase (Ambion) at 37°C for 30 min twice, and incubated with 1/10 (v/v) inactivation reagent (Ambion) to inactivate DNase. The recovered RNA was subjected to DNA contamination testing by PCR using gapA primers. One μg qualified RNA was next reverse transcribed into cDNA by the Superscript III First-Strand Synthesis System (Life Technologies) following the instructions of the manufacturer. Primers for qRT-PCR were designed with Primer3 software and are listed in Table S2. The gene, gapA, served as the endogenous control for normalization of target genes. qRT-PCR was set up with 5 μl Power SYBR Green PCR Master Mix (Life Technologies), 1 μl cDNA, 2 μM forward and reverse primers, and appropriate nuclease-free water to a total of 10 μl per reaction. The PCR was run on a 7500 Fast Real-time PCR System (ABI) with a program for 95°C 5 min for 1 cycle, and 95°C 15 s, 60°C 1 min for 35 cycles. The Ct-values from all qRT-PCR reactions in triplicate were analyzed to detect target gene expressions.

Small RNA or pNDM-HK sequences were taken as input queries and submitted to the Blast (Altschul et al., 1997; Tatusova and Madden, 1999; Brilli et al., 2008) website to search for plasmids that contained similar input sequences under stringent criteria (megablast). Multiple sRNAs were concatenated into one long query sequence. ClustalW was applied to generate multiple alignments and phylogenetic trees (Thompson et al., 1994; Chenna, 2003; Larkin et al., 2007). Only high coverage and identity (at least 70%) alignments were kept for phylogenetic tree construction employing the neighbor-joining method. Max sequence difference was set at 0.1. The resulting phylogenetic tree was plotted with MEGA6 (Tamura et al., 2013).

We previously sequenced the pNDM-HK plasmid isolated from an MDR E. coli strain in Hong Kong. The plasmid belonged to the broad host range IncL/M incompatibility group. In order to identify novel sRNAs encoded in this plasmid, we performed sRNA sequencing of the pNDM-HK transconjugant. The pNDM-HK plasmid was conjugated into E. coli K-12 J53 (Yi et al., 2012) and the transconjugant was selected and isolated as per our previous work (Ho et al., 2011). Over 150,000 reads from each sample were obtained. A majority of the reads (90%) were mapped to the coding regions (CDS) and known ncRNAs of the E. coli genome whereas 10% of the total reads were mapped onto the plasmid sequence (Table S3). Based on the read-mapping patterns, we identified six plasmid-encoded sRNAs that were all located at intergenic regions. Their details are summarized in Table 1. The sRNA sequences are presented in Table 2.

As depicted in Figure 1, these six plasmid-encoded sRNAs were distributed across different regions of the pNDM-HK plasmids, of which four (NDM-sR1, sR2, sR4, and sR5a-d) were located at the backbone region and two (NDM-sR3 and sR6) fell within the variable region. NDM-sR1 and sR2, which are located between repC and repB in the replication region, are two overlapping sRNAs with opposite orientations. They are 76 nt (NDM-sR1) and 64 nt (NDM-sR2) in length with 15 nt overlapping at the 3′-ends. Intriguingly, NDM-sR2 possessed the characteristic of the known counter-transcribed RNAs (ctRNA) because of the orientation (antisense RNA between repC and repB), sequence homology (Figure S1), and secondary structure (Athanasopoulos et al., 1995). As a negative post-transcriptional regulator of repA, ctRNA determines plasmid incompatibility and controls replication (Athanasopoulos et al., 1999; Izquierdo et al., 2005; Cervantes-Rivera et al., 2010). Unexpectedly, we identified a novel sRNA (NDM-sR1) that is antisense to this ctRNA, suggesting the possibility of a more intricate and complex replication control system in the pNDM-HK plasmid.

NDM-sR3 is a 76-nt sRNA located between ΔtnpATn1 and insL, suggesting the association of transposable-element shuffling or dissemination between plasmids. NDM-sR4 is a 101-nt sRNA in the stability region and downstream to a transposase gene tnpA of IS26. For NDM-sR5, it is noteworthy that four highly similar copies were found and deemed NDM-sR5a-d. The sequence of NDM-sR5b and -sR5c has up to 96% identity. The occurrence of NDM-sR5 at several genomic locations with high sequence similarity indicates that duplication events might have taken place during the course of evolution. As these isoforms are all located at the stability region of the plasmid, they may play roles in plasmid partitioning or post-segregational stability. NDM-sR6 was detected at the 3′-UTR of aacC2, a gene conferring resistance to aminoglycoside. Intriguingly, 3′-UTRs was reported as a source of regulatory RNAs in bacteria (Gößringer and Hartmann, 2012), and we found three transcription start sites (TSSs) of NDM-sR6 (Figure S2). These results suggest that the aacC2 gene may initiate other downstream consequences to the bacterial host through NDM-sR6. Whether NDM-sR6 is generated through internal processing of aacC2 or is transcribed independently requires further investigation.

In order to validate and characterize the newly found intergenic sRNAs encoded from the pNDM-HK plasmid, we performed Northern blot analysis to identify their expression patterns. Total RNAs were extracted from transconjugants grown at stationary phase and separated by denaturing gel electrophoresis. The RNAs were then transferred to the membrane, hybridized with the sRNA-specific radioactive labeled RNA probes and analyzed by autoradiography. As portrayed in Figure 2, all plasmid-encoded sRNAs were detected in the transconjugant but not in the J53 control, indicating that these sRNAs were produced from the pNDM-HK plasmid. The sRNAs ranged between 50 and 150 nt, and all of them showed a single major band as indicated by the Northern blot. The original images of the Northern blots are found in Figure S3.

To identify the TSS of these six novel sRNAs, 5′-RACEs of each newly transcribed sRNA was performed. As shown in Figure 3, the TSS of each sRNA (indicated by arrows) was identified on the basis of a stronger signal in Tobacco Acid Pyrophosphatase (TAP)-treated samples. All of these sRNAs exhibited single TSS except NDM-sR6 that possessed three different TSSs (Figure S2). The calculated sizes of NDM-sR6 were 147, 250, and 321 nucleotides. The secondary structures of the sRNAs were also predicted with the RNAfold program (Gruber et al., 2008) and are depicted in Figure 4. Most contain at least one stem loop and a stretch of uridine sequences, which indicates the typical structure of bacterial sRNA.

Phylogenetic profiling has been instrumental in inferring the evolutionary relationships between plasmids, especially the acquisition of genes by integration, transposition, and recombination. Phylogenetic trees are normally constructed from either the backbone sequence or the MDR regions of plasmids (Ma et al., 2007; Brilli et al., 2008; Norberg et al., 2011). In order to understand the contribution of plasmid-encoded sRNAs to the phylogenetic relationship of plasmids, we constructed and compared phylogenetic trees based on either the whole plasmid sequence or six novel sRNAs (refers to section Materials and Methods). As seen in Figure 5, the pNDM-HK and pNDM-OM plasmids are closely related phylogenetically, just as previously reported. Remarkably, we found that the plasmids in our sRNA-phylogenetic tree all belonged to the IncL/M group. Most of the resistance genes in pNDM-HK were acquired from these plasmids through transposition and recombination (Bonnin et al., 2013; Adamczuk et al., 2015). Plasmids such as pCTX-M3, pEl1573, and pNDM-OM were reported to be involved in the evolutionary pathway of pNDM-HK. Based on a recently published paper that re-designated IncL and IncM plasmids from the IncL/M group, most of the plasmids in our sRNA-phylogenetic tree were assigned to IncM2 (Carattoli et al., 2015) (pNDM-HK, pNDM-OM, pCTXM360, pEI1573, and pCTX-M3), except for plasmids without re-designation information (pIMP-HB623 and pKPC-CAV1741). We also constructed a phylogenetic tree based on either the backbone region of pNDM-HK or sRNAs at the backbone region, including NDM-sR1, 2, 4, and 5 (Figure S4). Although the resolution of the sRNA-phylogenetic tree is not as good as the tree constructed from the backbone region, we can still identify the plasmids closely related to pNDM-HK, such as pCTX-M360, pCTX-M3, pEl1573, and pNDM-OM. These results strongly suggest that plasmid-encoded sRNAs can be used for phylogenetic analysis to reveal evolutionary relationships in other MDR plasmids.

As illustrated in Figure 4, NDM-sR3 contains a polyU stretch at the 3′-end stem-loop and this stretch was commonly considered an Hfq-binding sequence (Otaka et al., 2011). As Hfq is the RNA chaperone that plays a significant role in sRNA biogenesis and function, we went on to assess if it interacts with NDM-sR3. We performed an Electrophoretic Mobility Shift Assay (EMSA) using in-vitro transcribed and radioactive-labeled NDM-sR3 RNA with purified recombinant Hfq protein. In Figure 6A, the NDM-sR3-Hfq complex was gradually formed and shifted in the gel as the concentration of Hfq increased. To quantitate the dissociation constant (K_d), the binding curve was plotted against an increasing concentration of Hfq based on the intensity of the bands (Figure 6B). The K_d of NDM-sR3 to Hfq was calculated as 750 nM. This suggests that NDM-sR3 interacts with Hfq at medium affinity compared with other sRNAs, of which the K_d are normally <100 nM (Henderson et al., 2013). We also tested the affinity between Hfq and another plasmid-encoded sRNA, NDM-sR2, as a control. Nearly no sR2-protein complex formed in the presence of 1,600 nM Hfq protein, indicating there was no interaction between Hfq and NDM-sR2 (Figure S5).

Hfq could protect the sRNA from ribonuclease degradation and increase its stability in bacteria. After showing that Hfq interacts with NDM-sR3, we further examined whether Hfq can protect NDM-sR3 in the bacterial host. The abundance of NDM-sR3 was monitored by Northern blot analysis in the presence (MG1655) and absence (MG1655Δhfq) of Hfq after the addition of rifampicin, a bacterial DNA-dependent RNA synthesis inhibitor. The results showed that NDM-sR3 was rapidly degraded in both strains, and it completely disappeared in 10 min after rifampicin treatment (Figure 7). The half-life of NDM-sR3 in MG1655 was determined as 2.19 ± 0.02 min, whereas that in MG1655Δhfq was 1.31 ± 0.33 min, indicative of NDM-sR3 being stabilized and protected by Hfq. The short half-life of NDM-sR3 is consistent with previously reported plasmid-encoded antisense RNA (RNAI) (Wagner and Brantl, 1998).

We have identified six pNDM-HK-encoded sRNAs, and they are constitutively expressed within the bacterial host. In order to further characterize the functional roles of sRNAs in the bacterial host, we selected NDM-sR3 that is located at the variable region and associated with IS26-mediated insertions in pCTX-M3 for further investigation. To determine the regulatory role of NDM-sR3 in the gene expression of bacteria, we cloned the TU of NDM-sR3 into the pTL01 plasmid and transformed it into the model strain of E. coli MG1655. We utilized a target gene prediction program, IntaRNA (Busch et al., 2008; Wright et al., 2014), to predict target genes on the bacterial host chromosome and validated their expression level by qRT-PCR in the NDM-sR3 overexpressed bacteria. As shown in Figure 8A, NDM-sR3 was successfully expressed and detected according to Northern blot analysis. Moreover, we found that overexpression of NDM-sR3 down-regulated the expression level of dinG, osmC, ptsI, and ybhF by approximately 50% (Figure 8B). dinG is an ATP-dependent helicase that confers DNA helicase activity in terms of DNA repair and replication. osmC is an osmotically inducible, stress-inducible membrane protein as well as peroxiredoxin (Cussiol et al., 2010; Cheng et al., 2012). ptsI is a cytoplasmic protein component of the phosphoenolpyruvate:sugar phosphotransferase system (Ginsburg and Peterkofsky, 2002). Finally, ybhF is predicted as part of the ATP-dependent efflux pump (Zhang et al., 2007).

As we have shown that NDM-sR3 is stabilized by Hfq protein in E. coli, we continued to evaluate whether the regulatory effects of NDM-sR3 on its target genes are also Hfq-dependent. We transformed NDM-sR3 into E. coli MG1655Δhfq and measured the expression levels of the target genes. Of note, none of the genes exhibited any influence in a regulatory sense in the absence of Hfq, indicating that the regulatory role of NDM-sR3 is Hfq-dependent.

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.