**Figure 1.** The *glmYqseEGF* gene locus in *P. mirabilis* N2. (**A**) The *P. mirabilis glmYqseEGF* gene locus corresponds to the similar locus in *E. coli* O157:H7 EDL933 and *Salmonella* Typhimurium LT2 with corresponding genes in shadows. An amino acid sequence analysis of the locus in these strains was performed using position-specific iterative BLAST. The percent amino acid identities of PurL, QseE, QseG, QseF and GlnB between *P*. *mirabilis* N2 and *E. coli* O157:H7 EDL933 or *Salmonella* Typhimurium LT2 are shown. The number above each arrow represents the gene length (bp). The intergenic space or overlap in terms of base pairs is also shown. (**B**) The promoter activities of *glmY* (−381 to −3 from putative *glmY* transcriptional start site) and *qseE* (−553 to −6 from *qseE* start codon) upstream region. (**C**) The *glmY*-*qseE-qseG-qseF* constitute an operon by gene-specific reverse transcription PCR. cDNA was synthesized by the *qseF*-specific primer and then the *glmY* DNA fragment was amplified by PCR. no RT, negative control; gDNA (genomic DNA), positive control. (**D**) The sequence of *glmYqseEGF* promoter region in *P. mirabilis* N2. The putative QseF and IHF binding sites are bold and underlined. The conserved σ70 and σ54 binding sites are indicated in boxes and shadows, respectively. **Figure 1.** The *glmYqseEGF* gene locus in *P. mirabilis* N2. (**A**) The *P. mirabilis glmYqseEGF* gene locus corresponds to the similar locus in *E. coli* O157:H7 EDL933 and *Salmonella* Typhimurium LT2 with corresponding genes in shadows. An amino acid sequence analysis of the locus in these strains was performed using position-specific iterative BLAST. The percent amino acid identities of PurL, QseE, QseG, QseF and GlnB between *P*. *mirabilis* N2 and *E. coli* O157:H7 EDL933 or *Salmonella* Typhimurium LT2 are shown. The number above each arrow represents the gene length (bp). The intergenic space or overlap in terms of base pairs is also shown. (**B**) The promoter activities of *glmY* (−381 to −3 from putative *glmY* transcriptional start site) and *qseE* (−553 to −6 from *qseE* start codon) upstream region. (**C**) The *glmY*-*qseE-qseG-qseF* constitute an operon by gene-specific reverse transcription PCR. cDNA was synthesized by the *qseF*-specific primer and then the *glmY* DNA fragment was amplified by PCR. no RT, negative control; gDNA (genomic DNA), positive control. (**D**) The sequence of *glmYqseEGF* promoter region in *P. mirabilis* N2. The putative QseF and IHF binding sites are bold and underlined. The conserved σ 70 and σ 54 binding sites are indicated in boxes and shadows, respectively.

#### *2.2. Identification of P. mirabilis glmY and Co-Transcription of glmY, qseE, qseG and qseF*

In view of conservation of the gene cluster *glmY-glrK* (*qseE*)*-yfhG* (*qseG*)*-glrR* (*qseF*)*-glnB* in *Enterobacteriaceae* [27] and modulation of motility [12] and pathogenesis [29,30] by GlmY, we searched the bacterial small RNA database (BSRD) for *P. mirabilis* GlmY counterpart and found it located upstream of *qseE* and downstream of *purL* (Figure 1A). We first examined the promoter activity in the upstream 548-bp region of *qseE* (nucleotide −553 to −6 of QseE start codon) by the *xylE* reporter assay. No promoter activity was found (Figure 1B). Then, we tested the *glmY* promoter activity using the upstream 379-bp DNA fragment (nucleotide −381 to −3 from the putative transcription start site) of *glmY* and found the promoter activity at 3, 5 and 7 h after incubation (Figure 1B). Furthermore, we demonstrated that *glmY*, *qseE*, *qseG* and *qseF* belong to a transcript by RT (reverse transcription)-PCR using *qseF* specific primer to obtain the cDNA for amplifying the *glmY* fragment (Figure 1C). The *glmY* PCR product of 159 bp was observed, but no product was produced in the no RT control (Figure 1C). These data indicate that *glmY*, *qseE*, *qseG* and *qseF* can share the same promoter upstream of *glmY*. The lack of promoter activity in the 548-bp DNA fragment upstream *qseE* results from the fragment containing only partial *glmY* promoter DNA sequences (201/379 bp). In summary, *glmY*, *qseE*, *qseG* and *qseF* are co-transcribed from the *glmY* promoter in *P. mirabilis.* In addition, we identified one overlapped σ 70/σ 54, two IHF and three QseF putative binding sites in the upstream promoter region of *glmY* (Figure 1D).

#### *2.3. Phenotypic Traits of glmY and qse Mutants*

To investigate the roles of GlmY and QseEF TCS in *P. mirabilis*, we first generated isogenic mutants of *glmY* and *qseF*. Among the phenotypes assayed, both *glmY* and *qseF* mutant strains exhibited reduced swarming and swimming abilities relative to the wildtype and respective complemented strains (Figure 2A,D). Comparable growth of the wildtype and respective complemented strains was observed, and no growth defect was found in the mutant strains compared to the wild-type (data not shown). *glmY* and *qseF* mutants migrated much slower than the wild-type and respective complemented strains during the 8-h period after inoculation on the swarming plate. Mutants lacking *qseE*, *qseG* or *qseEGF* were also constructed with no growth defect compared to the wild-type bacteria. We found similar swarming and swimming patterns among *qseF*, *qseE*, *qseG* and *qseEGF* mutants (data not shown). Therefore, we investigated the QseEGF-mediated motility using *qseF* mutant to represent QseEF malfunction in the following experiments. Significantly reduced swarming-related phenotypes, including cell differentiation, hemolysin activity and flagellin level, were observed in the *glmY* and *qseF* mutants compared to the wild-type and respective complemented strains (Figure 2B,C,E). Transmission electron microscopy (TEM) also showed *glmY* and *qseF* mutant cells were shorter and had fewer flagella than the wild-type strain (Figure 2F).

#### *2.4. Downregulation of glmYqseEGF Operon by qseF Deletion*

Based on the presence of the putative QseF binding site in the *glmYqseEGF* promoter region (Figure 1D), we examined whether *glmY* expression is under the control of QseF by the reporter assay. Just as expected, *qseF* mutant had significantly lower *glmY* promoter activity compared to the wild-type and the *qseF*-complemented strain at 3, 5, 7 h after inoculation and incubation (Figure 3A).

#### *2.5. Altered Expression of flhDC, rcsB and cheA in qseF and glmY Mutants*

With the notion that disruption of *qseF*, encoding a TCS transcriptional regulator QseF, impaired the swarming motility and related phenotypes in *P. mirabilis*, we thus examined the effect of *qseF* loss on *flhDC* expression by the reporter assay. The loss of *qseF* incurred significant reduction in *flhDC* promoter activity compared to the wild-type and the complemented strains at 3, 5, 7 h after inoculation (Figure 3B). Since QseF is a transcriptional regulator, we then tried to find the consensus QseF binding site (TGTCN10GACA) [31] in the promoter region of *flhDC* but failed. Since QseF could regulate *glmY* expression (Figure 3A) and both *qseF* and *glmY* mutants displayed decreased swarming and swimming motility (Figure 2A,D), we examined the promoter activity of *flhDC* in the *glmY* mutant. As is the case for *qseF* mutant, *glmY* mutant had a lower *flhDC* promoter activity than the wild-type and the complemented strains (Figure 3B). The mRNA level of *flhDC* in *qseF* and *glmY* mutants was consistent with the results of the reporter assay (Figure 3C).

**Figure 2.** GlmY and QseF regulated swarming and related phenotypes. Swarming motility (**A**), cell length (**B**), hemolysin activity (**C**) and swimming motility (**D**) were determined in the wild-type, mutants of *qseF* and *glmY* and complemented strains. Swarming migration was monitored at 1-h intervals, and swimming motility was determined after incubation for 16 h. Cell length and hemolysin activity were determined at 5 h after incubation on the LB agar plate. The data are the averages and standard deviations of three independent experiments. Significant difference of *qseF* or *glmY* mutant from the wild-type at 4, 5, 6, 7 and 8 h is indicated (\* *p* < 0.05; \*\* *p* < 0.01; \*\*\* *p* < 0.001 by the Student's *t* test) in (**A)**. Significant difference from the wild-type is indicated (\* *p* < 0.05; \*\*\* *p* < 0.001 by the Student's *t* test) in (**B**–**D**). (**E**) Analysis of flagellin expression by SDS-PAGE (upper panel) and Western blotting (lower panel). The flagellin level of the wild-type, mutants of *qseF* and *glmY* and complemented strains were examined at 5 h after seeding on the swarming plates by the SDS-PAGE and Western blotting as described in Materials and Methods. The representative picture of three independent experiments is shown. M, molecular weight marker. (**F**) TEM pictures of wildtype and mutants of *qseF* and *glmY*. Bacterial cultures were applied onto a carbon-coated grid, cells **Figure 2.** GlmY and QseF regulated swarming and related phenotypes. Swarming motility (**A**), cell length (**B**), hemolysin activity (**C**) and swimming motility (**D**) were determined in the wild-type, mutants of *qseF* and *glmY* and complemented strains. Swarming migration was monitored at 1-h intervals, and swimming motility was determined after incubation for 16 h. Cell length and hemolysin activity were determined at 5 h after incubation on the LB agar plate. The data are the averages and standard deviations of three independent experiments. Significant difference of *qseF* or *glmY* mutant from the wild-type at 4, 5, 6, 7 and 8 h is indicated (\* *p* < 0.05; \*\* *p* < 0.01; \*\*\* *p* < 0.001 by the Student's *t* test) in (**A**). Significant difference from the wild-type is indicated (\* *p* < 0.05; \*\*\* *p* < 0.001 by the Student's *t* test) in (**B**–**D**). (**E**) Analysis of flagellin expression by SDS-PAGE (upper panel) and Western blotting (lower panel). The flagellin level of the wild-type, mutants of *qseF* and *glmY* and complemented strains were examined at 5 h after seeding on the swarming plates by the SDS-PAGE and Western blotting as described in Materials and Methods. The representative picture of three independent experiments is shown. M, molecular weight marker. (**F**) TEM pictures of wild-type and mutants of *qseF* and *glmY*. Bacterial cultures were applied onto a carbon-coated grid, cells were stained with 1% PTA and TEM pictures were taken. Flagella are indicated by arrows. wt, wild-type; qseF, *qseF* mutant; glmY, *glmY* mutant; qseFc, *qseF*-complemented strain; glmYc, *glmY*-complemented strain.

were stained with 1% PTA and TEM pictures were taken. Flagella are indicated by arrows. wt, wildtype; qseF, *qseF* mutant; glmY, *glmY* mutant; qseFc, *qseF*-complemented strain; glmYc, *glmY*-com-

Based on the presence of the putative QseF binding site in the *glmYqseEGF* promoter region (Figure 1D), we examined whether *glmY* expression is under the control of QseF by the reporter assay. Just as expected, *qseF* mutant had significantly lower *glmY* promoter activity compared to the wild-type and the *qseF*-complemented strain at 3, 5, 7 h after

*2.4. Downregulation of glmYqseEGF Operon by qseF Deletion*

plemented strain.

**Figure 3.** QseF-controlled GlmY regulated expression of *flhDC*, *rcsB* and *cheA*. (**A**) The *glmY* promoter activity in the wild-type, *qseF* mutant and *qseF*-complemented strain. (**B**) The *flhDC* promoter activity in the wild-type, mutants of *qseF* and *glmY* and the complemented strains. (**C**) The mRNA level of *flhDC*, *rcsB* and *cheA* in the wild-type, mutants of *qseF* and *glmY* and the complemented strains. In (**A**,**B**), bacteria cultures were spread on LB agar plate and XylE promoter activities were determined as described in Materials and Methods after incubation for 3, 5 and 7 h at 37 °C. In (**C**), the mRNA amount was measured at 5 h after inoculation on the LB agar plate. The value for the wild-type was set at 1, and other data are presented relative to this value. The data are the averages and standard deviations of three independent experiments. Significant difference from the wildtype is indicated (\* *p* < 0.05; \*\* *p* < 0.01; \*\*\* *p* < 0.001 by the Student's *t* test). wt, wild-type; qseF, *qseF* mutant; glmY, *glmY* mutant; qseFc, *qseF-* complemented strain; glmYc, *glmY*-complemented strain. **Figure 3.** QseF-controlled GlmY regulated expression of *flhDC*, *rcsB* and *cheA*. (**A**) The *glmY* promoter activity in the wild-type, *qseF* mutant and *qseF*-complemented strain. (**B**) The *flhDC* promoter activity in the wild-type, mutants of *qseF* and *glmY* and the complemented strains. (**C**) The mRNA level of *flhDC*, *rcsB* and *cheA* in the wild-type, mutants of *qseF* and *glmY* and the complemented strains. In (**A**,**B**), bacteria cultures were spread on LB agar plate and XylE promoter activities were determined as described in Materials and Methods after incubation for 3, 5 and 7 h at 37 ◦C. In (**C**), the mRNA amount was measured at 5 h after inoculation on the LB agar plate. The value for the wild-type was set at 1, and other data are presented relative to this value. The data are the averages and standard deviations of three independent experiments. Significant difference from the wild-type is indicated (\* *p* < 0.05; \*\* *p* < 0.01; \*\*\* *p* < 0.001 by the Student's *t* test). wt, wild-type; qseF, *qseF* mutant; glmY, *glmY* mutant; qseFc, *qseF-* complemented strain; glmYc, *glmY*-complemented strain.

*2.5. Altered Expression of flhDC, rcsB and cheA in qseF and glmY Mutants*  With the notion that disruption of *qseF*, encoding a TCS transcriptional regulator QseF, impaired the swarming motility and related phenotypes in *P. mirabilis*, we thus examined the effect of *qseF* loss on *flhDC* expression by the reporter assay. The loss of *qseF* incurred significant reduction in *flhDC* promoter activity compared to the wild-type and the complemented strains at 3, 5, 7 h after inoculation (Figure 3B). Since QseF is a transcriptional regulator, we then tried to find the consensus QseF binding site (TGTCN10GACA) [31] in the promoter region of *flhDC* but failed. Since QseF could regulate *glmY* expression (Figure 3A) and both *qseF* and *glmY* mutants displayed decreased swarming and swimming motility (Figure 2A,D), we examined the promoter activity of *flhDC* in the *glmY* mutant. As is the case for *qseF* mutant, *glmY* mutant had a lower *flhDC* promoter activity than the wild-type and the complemented strains (Figure 3B). The mRNA level of *flhDC* in *qseF* and *glmY* mutants was consistent with the results of the reporter assay (Figure 3C). Given the positive regulation of swarming and *flhDC* expression by both QseF and GlmY, the regulatory hierarchy from QseF to GlmY and the absence of QseF binding site in *flhDC* promoter region, we surmised that QseF may modulate swarming through GlmY to affect *f lhDC* expression. How could GlmY, an sRNA, regulate both promoter activity and mRNA level of *flhDC* (Figure 3B,C)? In this respect, GlmY may modulate expression of a regulator which affects *flhDC* transcription directly and thus *flhDC* mRNA level. The direct regulation of *flhDC* mRNA level by GlmY may also contribute to the altered mRNA amount. As for *flhDC* regulator, *P. mirabilis* RcsB could modulate swarming by direct regulation of *flhDC* expression in a negative way [8,32] and *rcsB* overexpression led to a reduction of *flhDC* expression [33]. Therefore, we tested the effect of *glmY* mutation on *rcsB* mRNA level and found *glmY* mutant had increased *rcsB* mRNA level compared to the wild-type and *glmY*-complemented strain (Figure 3C). In view of regulation of *glmY* by QseF, *qseF* mutant exhibited a similar increase in *rcsB* mRNA level as the *glmY* mutant relative to the wild-type and the complemented strain (Figure 3C). On the other hand, we

tried to identify whether GlmY could target *f lhDC* or other mRNAs to modulate motility of *P. mirabilis*. We used the full sequence of GlmY as input to search for GlmY targets on the TargetRNA2 website. The tool did report a region in the 50UTR of the *cheA* (but not *f lhDC* or *rcsB*) as a candidate target with lower energy score and the *p*-value less than 0.05. The prediction shows interaction of GlmY (positions 17 to 28) with a 50UTR region of *cheA* mRNA spanning −61 to −50 positions from the translation start site of *cheA* mRNA. CheA kinase, encoded by *cheA*, belongs to a family of two-component sensors responsible for bacterial chemotaxis [34]. As *P. mirabilis* GlmY positively regulated swarming (Figure 2A) and small RNAs exert positive regulation of virulence primarily at the level of mRNA stabilization [35], we assessed whether GlmY affects *cheA* mRNA level. As expected, the amount of *cheA* mRNA was reduced significantly in *glmY* and *qseF* mutants compared to the wild-type and respective complemented strains (Figure 3C).

#### *2.6. GlmY Activates cheA Expression at the Post-Transcriptional Level*

To investigate the role of interaction between GlmY and *cheA* 5 0UTR in *cheA* expression and thus swarming, we generated plasmids carrying the individual *lac* promoter-driven *xylE* translational fusion containing different lengths (520, 253, 138 and 61 bp) of *cheA* 5 0UTR and the first 27 bp of *cheA* ORF (Figure 4A). The individual fusion plasmid was transformed into the wild-type strain and *glmY* mutant to probe the essential region of *cheA* 5 0UTR requiring GlmY for activating *cheA* translation by monitoring the *xylE* activity. While the *glmY* mutant carrying the plasmid of 520, 253 or 138 bp-50UTR translational fusion exhibited significantly lower *xylE* activity than the wild-type strain carrying the respective fusion plasmid, albeit to a lesser extent for the 138 bp-fusion plasmid (Figure 4B), the *xylE* activity of the wild-type strain and the *glmY* mutant harboring the 61 bp-fusion plasmid was comparable (Figure 4B). The results indicate that there is no need for GlmY to activate *cheA* translation when only the 61 bp-50UTR of *cheA* is present. The data indicate all the 520, 253 and 138-bp fragments of *cheA* 5 0UTR had constrained structures requiring GlmY for releasing to activate *cheA* translation.

To further demonstrate whether GlmY affected *cheA* expression via direct base-pairing, we performed site-directed mutagenesis to inactivate the interaction of GlmY and *cheA* 5 0UTR. We introduced an 8-bp mutation into the predicted pairing regions of the *cheA* 5 0UTR (cheAm) and *glmY* (GlmYm) (Figure 4C) using the *cheA* 5 0UTR (253 bp)-*xylE* translational fusion plasmid (pcheA) and *glmY*-containing pBAD33 (pBglmY) as templates to generate pcheAm and pBglmYm, respectively. The wild-type *P. mirabilis* harboring the pcheAm exhibited significantly lower XylE activity compared to that harboring the wild-type *cheA* 5 0UTR-*xylE* plasmid (pcheA) (Figure 4D). The *glmY* mutant carrying pBglmY and pcheAm showed significantly lower XylE activity than the mutant carrying pBglmY and pcheA (0.4 vs. 1 in Figure 4D). GlmY with the 8-bp substitution (GlmYm) could not effectively enhance the expression of *cheA*-*xylE* fusion to the extent of wild-type GlmY (0.6 vs. 1) in the *glmY* mutant (Figure 4D). In addition, introduction of compensating mutations into the *cheA* 5 0UTR could restore the ability of GlmYm to activate the *cheA-xylE* fusion to 80% of the wild-type GlmY and *cheA* interaction (pBglmYm-pcheAm vs. pBglmY-pcheA in Figure 4D).

**Figure 4.** GlmY activated *cheA* expression at the post-transcriptional level. (**A**) Schematic representation of the different constructs used in the *cheA-xylE* translational reporter assay. The various lengths of *cheA* 5′UTR (520, 253, 138 and 61 bp before AUG plus 27 bp after AUG were in frame ligated to *xylE* gene to obtain the translational fusion under the control of *lac* promoter in the pGEM-T easy vector. (**B**) The effect of various lengths of *cheA* 5′UTR (as shown in (**A**)) on XylE activity of the *cheA* 5′UTR-*xylE* translational fusion in the wild-type and *glmY* mutant. The activity of XylE in the translational *cheA* 5′UTR-*xylE* reporter plasmid-transformed *P. mirabilis* strains was determined using the reporter assay at 5 h after incubation on an LB agar plate. The value obtained for the wildtype strain at a length of 5′UTR was set at 1. The data are the averages and standard deviations of three independent experiments. Significant difference from the wild-type is indicated (\* *p* < 0.05; \*\* *p* < 0.01 by the Student's *t* test). (**C**) The predicted region of base-pairing between GlmY and the 5′ UTR of the *cheA* mRNA. The interaction site of *cheA* and GlmY are underlined and the substitutions present in GlmY (GlmYm) and *cheA* 5′ UTR (*cheA*m) are shown. (**i**). Wild-type GlmY and *cheA* 5′ UTR; (**ii**) wild-type GlmY and mutated *cheA* 5′ UTR; (**iii**) wild-type *cheA* 5′ UTR and mutated GlmY; (**iv**) mutated GlmY and the compensatory mutations in *cheA* 5′ UTR. (**D**) Analysis of GlmY interaction with *cheA* 5′ UTR for *cheA* expression. The *cheA* 5′UTR (253 bp)-*xylE* translational fusion plasmid (pcheA) and *glmY*-containing pBAD33 (pBglmY) were used as templates for introducing an 8-bp substitution (shown in (**C**)) into GlmY and *cheA* 5′ UTR by site-directed mutagenesis to produce pcheAm and pBglmYm, respectively. The pcheA and pcheAm were introduced into wild-type *P. mirabilis* separately*,* while *glmY* mutant was transformed with combinations of pBglmY-pcheA, pBglmY-pcheAm, pBglmYm-pcheA or pBglmYm-pcheAm. Then, the activities of XylE in the various *glmY* mutants (in the presence of 0.2% arabinose) and wild-types were determined at 5 h after incubation on an LB agar plate. The relative XylE activity was shown with the value of the pcheAharbored wild-type or *glmY* mutant carrying pBglmY and pcheA set at 1, respectively. Significant difference from the wild-type or *glmY* mutant set at 1 is indicated (\* *p* < 0.05; \*\* *p* < 0.01 by the Student's *t* test). wt, wild-type; glmY, *glmY* mutant. To further demonstrate whether GlmY affected *cheA* expression via direct base-pairing, we performed site-directed mutagenesis to inactivate the interaction of GlmY and **Figure 4.** GlmY activated *cheA* expression at the post-transcriptional level. (**A**) Schematic representation of the different constructs used in the *cheA-xylE* translational reporter assay. The various lengths of *cheA* 5 0UTR (520, 253, 138 and 61 bp before AUG plus 27 bp after AUG were in frame ligated to *xylE* gene to obtain the translational fusion under the control of *lac* promoter in the pGEM-T easy vector. (**B**) The effect of various lengths of *cheA* 5 0UTR (as shown in (**A**)) on XylE activity of the *cheA* 5 0UTR-*xylE* translational fusion in the wild-type and *glmY* mutant. The activity of XylE in the translational *cheA* 5 0UTR-*xylE* reporter plasmid-transformed *P. mirabilis* strains was determined using the reporter assay at 5 h after incubation on an LB agar plate. The value obtained for the wild-type strain at a length of 50UTR was set at 1. The data are the averages and standard deviations of three independent experiments. Significant difference from the wild-type is indicated (\* *p* < 0.05; \*\* *p* < 0.01 by the Student's *t* test). (**C**) The predicted region of base-pairing between GlmY and the 50 UTR of the *cheA* mRNA. The interaction site of *cheA* and GlmY are underlined and the substitutions present in GlmY (GlmYm) and *cheA* 5 0 UTR (*cheA*m) are shown. (**i**). Wild-type GlmY and *cheA* 5 0 UTR; (**ii**) wild-type GlmY and mutated *cheA* 5 0 UTR; (**iii**) wild-type *cheA* 5 0 UTR and mutated GlmY; (**iv**) mutated GlmY and the compensatory mutations in *cheA* 5 0 UTR. (**D**) Analysis of GlmY interaction with *cheA* 5 0 UTR for *cheA* expression. The *cheA* 5 0UTR (253 bp)-*xylE* translational fusion plasmid (pcheA) and *glmY*-containing pBAD33 (pBglmY) were used as templates for introducing an 8-bp substitution (shown in (**C**)) into GlmY and *cheA* 5 0 UTR by site-directed mutagenesis to produce pcheAm and pBglmYm, respectively. The pcheA and pcheAm were introduced into wild-type *P. mirabilis* separately, while *glmY* mutant was transformed with combinations of pBglmY-pcheA, pBglmY-pcheAm, pBglmYm-pcheA or pBglmYm-pcheAm. Then, the activities of XylE in the various *glmY* mutants (in the presence of 0.2% arabinose) and wild-types were determined at 5 h after incubation on an LB agar plate. The relative XylE activity was shown with the value of the pcheA-harbored wild-type or *glmY* mutant carrying pBglmY and pcheA set at 1, respectively. Significant difference from the wild-type or *glmY* mutant set at 1 is indicated (\* *p* < 0.05; \*\* *p* < 0.01 by the Student's *t* test). wt, wild-type; glmY, *glmY* mutant.

These results prove the role of *cheA* 5 0 UTR in interaction with GlmY and the 8-bp direct pairing between GlmY and *cheA* 5 0 UTR critical for facilitating expression of *cheA*. For the first time, this work identifies *cheA* as a novel target of GlmY for modulating swarming motility of *P. mirabilis*. *2.7. Urea Inhibited Promoter Activity of the glmYqseEGF Operon in P. mirabilis*  Based on the fact that *P. mirabilis* swarming was subject to positive regulation of the

#### *2.7. Urea Inhibited Promoter Activity of the glmYqseEGF Operon in P. mirabilis* QseF-GlmY pathway, it was tempting to determine whether urea affect the expression of *glmYqseEGF* operon in the urinary tract, an environment containing a lot of urea. There-

pBglmY-pcheA in Figure 4D).

swarming motility of *P. mirabilis*.

*Int. J. Mol. Sci.* **2022**, *23*, x FOR PEER REVIEW 9 of 17

*cheA* 5′UTR. We introduced an 8-bp mutation into the predicted pairing regions of the *cheA* 5′UTR (cheAm) and *glmY* (GlmYm) (Figure 4C) using the *cheA* 5′UTR (253 bp)-*xylE* translational fusion plasmid (pcheA) and *glmY*-containing pBAD33 (pBglmY) as templates to generate pcheAm and pBglmYm, respectively. The wild-type *P. mirabilis* harboring the pcheAm exhibited significantly lower XylE activity compared to that harboring the wild-type *cheA* 5′UTR-*xylE* plasmid (pcheA) (Figure 4D). The *glmY* mutant carrying pBglmY and pcheAm showed significantly lower XylE activity than the mutant carrying pBglmY and pcheA (0.4 vs. 1 in Figure 4D). GlmY with the 8-bp substitution (GlmYm) could not effectively enhance the expression of *cheA*-*xylE* fusion to the extent of wild-type GlmY (0.6 vs. 1) in the *glmY* mutant (Figure 4D). In addition, introduction of compensating mutations into the *cheA* 5′UTR could restore the ability of GlmYm to activate the *cheAxylE* fusion to 80% of the wild-type GlmY and *cheA* interaction (pBglmYm-pcheAm vs.

These results prove the role of *cheA* 5′ UTR in interaction with GlmY and the 8-bp

direct pairing between GlmY and *cheA* 5′ UTR critical for facilitating expression of *cheA*. For the first time, this work identifies *cheA* as a novel target of GlmY for modulating

Based on the fact that *P. mirabilis* swarming was subject to positive regulation of the QseF-GlmY pathway, it was tempting to determine whether urea affect the expression of *glmYqseEGF* operon in the urinary tract, an environment containing a lot of urea. Therefore, the XylE activities of the *glmY-xylE* reporter plasmid-transformed wild-type and the *qseF* mutant were monitored in the presence of urea. The promoter activity of the *glmY* operon was reduced by urea at 50 mM in the wild-type but not the *qseF* mutant after incubation for 5 h on the LB agar plate (Figure 5A). We then examined the swarming motility of the wild-type and *qseF* mutant in the presence of urea or not. The swarming motility of the wild-type but not the *qseF* mutant was decreased by urea at 50 mM (Figure 5B). The results suggest that urea could be a negative signal for expression of *P. mirabilis glmYqseEGF* operon in the urinary tract. fore, the XylE activities of the *glmY-xylE* reporter plasmid-transformed wild-type and the *qseF* mutant were monitored in the presence of urea. The promoter activity of the *glmY*  operon was reduced by urea at 50 mM in the wild-type but not the *qseF* mutant after incubation for 5 h on the LB agar plate (Figure 5A). We then examined the swarming motility of the wild-type and *qseF* mutant in the presence of urea or not. The swarming motility of the wild-type but not the *qseF* mutant was decreased by urea at 50 mM (Figure 5B). The results suggest that urea could be a negative signal for expression of *P. mirabilis glmYqseEGF* operon in the urinary tract.

**Figure 5.** Urea inhibited *glmY* expression through QseF. (**A**) *glmY* promoter activity in the wild-type and *qseF* mutant in the presence and absence of urea. Bacterial cultures were spread on an LB agar plate with or without 50 mM of urea, and XylE promoter activities were determined at 5 h after incubation at 37 °C. (**B**) Swarming motility of the wild-type and *qseF* mutant in the presence and absence of 50 mM of urea was monitored at 1-h intervals. The data are the averages and standard deviations of three independent experiments. Significant difference between presence and absence of urea is indicated (\* *p* < 0.05; \*\* *p* < 0.01 by the Student's *t* test). wt, wild-type; qseF, *qseF* mutant; **Figure 5.** Urea inhibited *glmY* expression through QseF. (**A**) *glmY* promoter activity in the wild-type and *qseF* mutant in the presence and absence of urea. Bacterial cultures were spread on an LB agar plate with or without 50 mM of urea, and XylE promoter activities were determined at 5 h after incubation at 37 ◦C. (**B**) Swarming motility of the wild-type and *qseF* mutant in the presence and absence of 50 mM of urea was monitored at 1-h intervals. The data are the averages and standard deviations of three independent experiments. Significant difference between presence and absence of urea is indicated (\* *p* < 0.05; \*\* *p* < 0.01 by the Student's *t* test). wt, wild-type; qseF, *qseF* mutant; nil, no urea.

#### nil, no urea. **3. Discussion**

For the first time, in this study, a TCS regulator QseF participating in modulation of swarming motility through GlmY (an sRNA) and the underlying mechanisms were revealed in uropathogenic *P. mirabilis*. GlmY regulated swarming through direct GlmY-*cheA* interaction together with GlmY-*rcsB-* and/or GlmY-mediated *flhDC* expression (Figure 6), whereby affecting swarming-related chemotaxis system and flagellum production (Figure 6). The chemotaxis system plays an essential role in flagellar function and swarm cell differentiation, thus important for *P. mirabilis* to display a vigorous swarming pattern [34,36]. Given decreased not abolished GlmY expression upon loss of *qseF* (Figure 3A), this indicates GlmY is subject to QseF-independent control. Moreover, the finding that loss of *glmY* or *qseF* resulting in a similar motility ability, expression of *flhDC* and *rcsB* and *flhDC*-associated phenotypes (cell length and flagellin level) (Figures 2 and 3) indicates QseF could have

alternative regulation of *flhDC* expression bypassing GlmY (Figure 6). In this regard, our preliminary data showed introduction of GlmY-expressing plasmid into *qseF* mutant could not restore *flhDC* mRNA to the wild-type level. have alternative regulation of *flhDC* expression bypassing GlmY (Figure 6). In this regard, our preliminary data showed introduction of GlmY-expressing plasmid into *qseF* mutant could not restore *flhDC* mRNA to the wild-type level.

For the first time, in this study, a TCS regulator QseF participating in modulation of

swarming motility through GlmY (an sRNA) and the underlying mechanisms were revealed in uropathogenic *P. mirabilis*. GlmY regulated swarming through direct GlmY-*cheA* interaction together with GlmY-*rcsB-* and/or GlmY-mediated *flhDC* expression (Figure 6), whereby affecting swarming-related chemotaxis system and flagellum production (Figure 6). The chemotaxis system plays an essential role in flagellar function and swarm cell differentiation, thus important for *P. mirabilis* to display a vigorous swarming pattern [34,36]. Given decreased not abolished GlmY expression upon loss of *qseF* (Figure 3A), this indicates GlmY is subject to QseF-independent control. Moreover, the finding that loss of *glmY* or *qseF* resulting in a similar motility ability, expression of *flhDC* and *rcsB* and *flhDC*-associated phenotypes (cell length and flagellin level) (Figures 2 and 3) indicates QseF could

*Int. J. Mol. Sci.* **2022**, *23*, x FOR PEER REVIEW 10 of 17

**3. Discussion**

**Figure 6.** Summary of the swarming regulation by the QseEF-GlmY pathway involving expression of *rcsB*, *flhDC* and *cheA*. The sensor kinase QseE senses specific signals and transfers phosphoryl group to response regulator QseF. Then, QseF activates expression of GlmY, followed by *cheA* activation and *rcsB* inhibition. RcsB downregulation and CheA upregulation facilitate flagellum production and swarmer cell differentiation, respectively, thereby enhancing swarming motility. QseG could help activation of QseEF and unknown factor X facilitates GlmY expression. GlmY may interact with *flhDC* and modulate its expression. In addition, QseF could affect *flhDC* expression in a GlmY-independent way. Arrow, positive effect; line with a vertical bar, negative effect; dotted line, direct effect needed to be investigated. **Figure 6.** Summary of the swarming regulation by the QseEF-GlmY pathway involving expression of *rcsB*, *flhDC* and *cheA*. The sensor kinase QseE senses specific signals and transfers phosphoryl group to response regulator QseF. Then, QseF activates expression of GlmY, followed by *cheA* activation and *rcsB* inhibition. RcsB downregulation and CheA upregulation facilitate flagellum production and swarmer cell differentiation, respectively, thereby enhancing swarming motility. QseG could help activation of QseEF and unknown factor X facilitates GlmY expression. GlmY may interact with *flhDC* and modulate its expression. In addition, QseF could affect *flhDC* expression in a GlmY-independent way. Arrow, positive effect; line with a vertical bar, negative effect; dotted line, direct effect needed to be investigated.

An ever-increasing number and variety of sRNAs are being identified to serve regulatory functions for bacteria to respond to environmental cues and thrive in diverse habitats [9–11,37]. *E. coli* GlmY has been coopted to modulate the expression of virulence and be involved in cellular metabolism and architecture, including for biosynthesis of LPS [38– 40], a permeability barrier and a major virulence determinant in pathogenic bacteria [41]. GlmY fine-tunes expression of type III secretion system and its effectors to promote bacterial attachment and subsequent actin rearrangement on host cells through post-transcriptional control of EspFu and the locus-of-enterocyte-effacement (LEE) [29]. In addition, GlmY and GlmZ participate in gene expression of curli adhesion, acid resistance and also tryptophan metabolism [42]. Despite similar genomic gene arrangement of *glmYqseEGF*, there are discrepancies of the QseF and GlmY-related regulation between *P. mirabilis* and *E. coli*. First, a Rho-independent terminator exists in the end of *E. coli* GlmY An ever-increasing number and variety of sRNAs are being identified to serve regulatory functions for bacteria to respond to environmental cues and thrive in diverse habitats [9–11,37]. *E. coli* GlmY has been coopted to modulate the expression of virulence and be involved in cellular metabolism and architecture, including for biosynthesis of LPS [38–40], a permeability barrier and a major virulence determinant in pathogenic bacteria [41]. GlmY fine-tunes expression of type III secretion system and its effectors to promote bacterial attachment and subsequent actin rearrangement on host cells through post-transcriptional control of EspFu and the locus-of-enterocyte-effacement (LEE) [29]. In addition, GlmY and GlmZ participate in gene expression of curli adhesion, acid resistance and also tryptophan metabolism [42]. Despite similar genomic gene arrangement of *glmYqseEGF*, there are discrepancies of the QseF and GlmY-related regulation between *P. mirabilis* and *E. coli*. First, a Rho-independent terminator exists in the end of *E. coli* GlmY [43] but not in that of *P. mirabilis*. Second, neither promoter prediction nor reporter assay showed a promoter present in the intergenic region between *glmY* and *qseE* of *P. mirabilis*, which is not the case for *E. coli*, initiating *qseE* transcription in the *glmY-qseE* intergenic region [31]. Third, *E. coli* QseEF is involved in regulating genes required for pedestal formation but not motility [21]. There was no difference in motility and flagellar expression as seen by Western blotting between the wild-type *E. coli* strain and the *qseE* mutant [21], contrary to the similar defect of both *P. mirabilis qseF* and *qseE* mutants in swarming and swimming abilities. *E. coli glmY* and *qseE* are independently transcribed from different promoters [31], while we

demonstrated that *glmY* and *qseEGF* constitute an operon by RT-PCR assay (Figure 1C). In this regard, the presence of the conserved RNase E-cleavage motif GCCUUAU in GlmY of *P. mirabilis* [38] indicates the *glmYqseEGF* transcript could be processed to produce effective GlmY for its function.

We found swarming, swimming and swarming-related phenotypes (cell length, haemolysin activity and flagellin level) were all reduced in both *qseF* and *glmY* mutants (Figure 2). Both *flhDC* promoter activity and mRNA level were downregulated in *qseF* and *glmY* mutants (Figure 3). Because no QseF binding site in the promoter region of *flhDC*, we speculated QseF may exert its effect on *flhDC* expression through GlmY. TargetRNA2 tool, focusing its search for an sRNA-mRNA interaction in a neighborhood around the rbs of the mRNA, revealed no GlmY binding site in *flhDC* and *rcsB* mRNAs. We then extended mRNA searching by IntaRNA for GlmY-interacting sites and binding sites of −101 to −92 and −284 to −271 from AUG (both outside rbs) for *flhDC* and *rcsB*, respectively, were identified. Based on the GlmY target prediction, there are two possibilities for upregulation of *flhDC* mRNA level by GlmY. One is by direct interaction with *flhDC* mRNA; the other is indirectly through *rcsB*. It also can not be ruled out that both modes of action coexist. Further studies are needed to investigate whether GlmY directly interacts with mRNAs of *flhDC* and *rcsB*. For example, the *flhDC* (or *rcsB*) mRNA amount of the *glmY* mutant harboring wild-type or interacting site-mutated *glmY* on a plasmid will be determined to assess the GlmY-*flhDC* (or *rcsB*) mRNA interaction. It is noteworthy that overexpression of *glmY* in the same *E. coli* strain either had no effect or resulted in motility repression [12,44]. It could be due to different plasmid vectors and assay conditions used.

Translational fusion assay indicated that regions of −520, −253 and −138 to +27 from AUG of *cheA* mRNA contain constrained secondary structures needed to be resolved by GlmY for *cheA* translation, albeit to a lesser extent for the region of −138 to +27 (Figure 4B). To elucidate how *cheA* 5 0UTR interacts with GlmY, we uploaded sequences from −253 to −138 and −138 to +27 of *cheA* mRNA to the IntaRNA webserver and interaction between regions of −233 to −211 and −111 to −94 from AUG was revealed (Figure 7). Likewise, using sequences from −138 to −61 and −61 to +27 of *cheA* mRNA as input fragments showed interaction between regions of −88 to −82 and −15 to −9 by IntaRNA (Figure 7). Inspection of the sequence of −15 to −9 disclosed the existence of two overlapped putative ribosome binding sites (rbs), aagguga (gagguga in *E. coli*) and gaaugag (gaagga in *E. coli*) for translation of long *cheA* and short *cheA* of *E. coli* [45], respectively (Figure 7). Furthermore, IntaRNA revealed interaction of *cheA* −61 to −50 from AUG with GlmY +17 to +28 (Figure 7) by using full-length GlmY and −253 to +27 of *cheA* as the input fragments. These data suggest interaction between regions of −88 to −82 and −15 to −9 from AUG could inhibit *cheA* translation by hiding rbs. Hence GlmY is required for releasing the secondary structure to assist *cheA* translation. The reason for the less extent of translation affected in the absence of *glmY* using the translational fusion comprising −138 to +27 of *cheA* compared to the −253 to +27 of *cheA* fusion could be ascribed to the interaction between regions of −233 to −211 and −111 to −94 from AUG (Figure 7), thereby affecting the ease for GlmY to uncover the rbs for translation. It is interesting to know that GlmY not only activated *cheA* translation but also maintained *cheA* mRNA level. This is in line with that sRNA-mediated stability control is the crucial element of activation of *trans*-encoded mRNAs [35].

EHEC regulated pathogenesis and motility by sensing epinephrine or norepinephrine through QseBC and QseEF two-component signaling systems [21,22]. The membrane kinase QseC autophosphorylates and phosphorylates the QseB response regulator initiating a signaling cascade that activates QseEF to trigger expression of LEE genes, leading to AE lesions on intestinal epithelial cells. In *P. mirabilis* genome, only QseEF homologue of high similarity was found, with no homologue of QseBC existing. We found QseF participated in swarming regulation and deletion of *qseE* or *qseG* also had decreased swarming motility. Previous study has revealed the phosphorylation state of QseE and QseF is governed by interaction with QseG in response to epinephrine for post-transcriptional regulation of

virulence genes through GlmY [46]. In addition, QseEGF has been shown to modulate transcription of *phoPQ*, linking to the virulence regulation [47]. The similarity of the transcriptome profiles of *qseE*, *qseF* and *qseG* mutants also indicates that these proteins work together [47]. Therefore, *P. mirabilis* QseEF likely should exert functions other than swarming to affect virulence. Accordingly, our preliminary data showing *qseF* mutant had a significantly impaired ability to colonize mouse bladders and kidneys compared to the wild-type. The investigation of virulence traits such as cytotoxicity, urothelial cell invasion and survival in macrophages is underway. *Int. J. Mol. Sci.* **2022**, *23*, x FOR PEER REVIEW 12 of 17

**Figure 7.** A model for regulation of *cheA* by GlmY sRNA. (**A**) The model structure for the 5′ UTR of *cheA* mRNA in the absence of GlmY. The expression of *cheA* is silenced at the post-transcriptional level by pairing −111 to −94 with −233 to −210 from AUG of *cheA* 5′ UTR and −88 to −82 from AUG with the putative ribosome binding site of *cheA*. (**B**) The putative structure for the 5′ UTR of *cheA* mRNA in the presence of GlmY. GlmY pairs with −61 to −50 from AUG of *cheA* 5′ UTR freeing the ribosome binding site and causing the *cheA* mRNA to be translated. Two overlapped putative ribosome binding sequences (AAGGUGA and GAAUGAG) are indicated in boxes, and the RNA-RNA base-pairing sequences predicted by intaRNA are underlined. **Figure 7.** A model for regulation of *cheA* by GlmY sRNA. (**A**) The model structure for the 50 UTR of *cheA* mRNA in the absence of GlmY. The expression of *cheA* is silenced at the post-transcriptional level by pairing −111 to −94 with −233 to −210 from AUG of *cheA* 5 0 UTR and −88 to −82 from AUG with the putative ribosome binding site of *cheA*. (**B**) The putative structure for the 50 UTR of *cheA* mRNA in the presence of GlmY. GlmY pairs with −61 to −50 from AUG of *cheA* 5 0 UTR freeing the ribosome binding site and causing the *cheA* mRNA to be translated. Two overlapped putative ribosome binding sequences (AAGGUGA and GAAUGAG) are indicated in boxes, and the RNA-RNA base-pairing sequences predicted by intaRNA are underlined.

EHEC regulated pathogenesis and motility by sensing epinephrine or norepinephrine through QseBC and QseEF two-component signaling systems [21,22]. The membrane kinase QseC autophosphorylates and phosphorylates the QseB response regulator initiating a signaling cascade that activates QseEF to trigger expression of LEE genes, leading to AE lesions on intestinal epithelial cells. In *P. mirabilis* genome, only QseEF homologue of high similarity was found, with no homologue of QseBC existing. We found QseF participated in swarming regulation and deletion of *qseE* or *qseG* also had decreased swarming motility. Previous study has revealed the phosphorylation state of QseE and QseF is governed by interaction with QseG in response to epinephrine for post-transcriptional regulation of virulence genes through GlmY [46]. In addition, QseEGF has been shown to modulate transcription of *phoPQ*, linking to the virulence regulation [47]. The similarity of the transcriptome profiles of *qseE*, *qseF* and *qseG* mutants also indicates that these proteins work together [47]. Therefore, *P. mirabilis* QseEF likely should exert functions other than swarming to affect virulence. Accordingly, our preliminary data showing *qseF* mutant had a significantly impaired ability to colonize mouse bladders and kidneys compared to the wild-type. The investigation of virulence traits such as cytotoxicity, urothelial cell invasion and survival in macrophages is underway. We found that urea could serve as a negative signal of QseEF for swarming (Figure We found that urea could serve as a negative signal of QseEF for swarming (Figure 5). In this way, the expression of QseF will be inhibited in the urine (rich in urea) and the role of QseF in facilitating swarming will be neglected. In view of our preliminary data showing a significant difference between the wild-type and *qseF* mutant in colonization of mouse bladders and kidneys, it is tempting to surmise there is other cues in the urine to increase expression of *qseF*. In this aspect, our unpublished transcriptome data reveal mRNA levels of fatty acid synthetic genes increase but those of degradative genes decrease in the *qseF* mutant relative to the wild-type. This indicates QseF may sense fatty acid and be involved in fatty acid metabolism. Interestingly, we found oleic acid (0.01%) is a positive signal for *qseF* expression (data not shown) and swarming motility of *P. mirabilis* [48]. The oleic acid concentration used is an attainable concentration in the urine according to the Human Metabolome Database (https://hmdb.ca/metabolites?utf8=%E2%9C%93&quantified=1& urine=1&filter=true, accessed on 29 November 2021). Given QseE is also phosphatase [23], it is reasonable to infer that the phosphatase function of QseE prevents QseF activation under the presence of urea, whereas its kinase function triggers QseF regulon expression under appropriate concentrations of oleic acid. A finely tuned balance in these opposing activities should determine the regulon response of QseF as the case of the bifunctional DevS kinase responsive to environmental oxygen in *Mycobacterium tuberculosis* [49].

5). In this way, the expression of QseF will be inhibited in the urine (rich in urea) and the role of QseF in facilitating swarming will be neglected. In view of our preliminary data showing a significant difference between the wild-type and *qseF* mutant in colonization of mouse bladders and kidneys, it is tempting to surmise there is other cues in the urine to increase expression of *qseF*. In this aspect, our unpublished transcriptome data reveal Since catheter-associated UTI (CAUTI) is a major health concern, research directed at understanding the pathogenesis is warranted and should lead to improved diagnosis, prevention and treatment. *P. mirabilis* is notorious for causing CAUTIs. In this work, we demonstrated a new regulatory pathway involving an sRNA, GlmY, participating in swarming motility of *P. mirabilis* during which expression of several virulence genes is

mRNA levels of fatty acid synthetic genes increase but those of degradative genes decrease in the *qseF* mutant relative to the wild-type. This indicates QseF may sense fatty

*rabilis* [48]. The oleic acid concentration used is an attainable concentration in the urine according to the Human Metabolome Database (https://hmdb.ca/metabolites?utf8=%E2%9C%93&quantified=1&urine=1&filter=true, accessed on 30 November 2021). Given QseE is also phosphatase [23], it is reasonable to infer that the phosphatase

acid and be involved in fatty acid metabolism. Interestingly, we found oleic acid (0.01%) is a positive signal for *qseF* expression (data not shown) and swarming motility of *P. mi-*272

increased. It is believed that *P. mirabilis* swarming up catheters is primed to infect the urinary tract [50], so elucidating the swarming mechanisms could provide new approaches in the development of intervention strategies and facilitate the discovery of novel therapeutics.

#### **4. Materials and Methods**

#### *4.1. Bacterial Strains, Plasmids, Reagents and Growth Conditions*

The bacterial strains and plasmids used in this study are listed in Table S1 in the Supplementary Materials. The bacterial strains used are a clinical isolate from a patient of UTI (the wild-type N2), its derived mutants and respective complemented strains. All chemicals were obtained from the Sigma-Aldrich unless otherwise indicated, and primer sequences are given in Table S2. Bacteria were stored at −80 ◦C and routinely cultured in Luria-Bertani (LB) broth at 37 ◦C. The LSW- agar plate [15] was used to prevent the phenotypic expression of swarming motility for selecting mutant clones and colony counting.

#### *4.2. Construction of P. mirabilis Mutants and Complemented Strains*

Sequences flanking the *qseF* gene was amplified by PCR using the primer pairs qseFupF/XbaI-qseF-upR and XbaI-qseF-dnF/qseF-dnR for *qseF* mutant, and cloned into pGEM-T Easy (Promega) to generate pGqseF-up and pGqseF-dn. pGqseF-up was digested with SalI/XbaI, and the *qseF* upstream sequence-containing fragment was ligated to SalI/XbaIdigested pGqseF-dn to produce the pGqseF-updn plasmid, which contains both upstream and downstream sequences of *qseF*. A Kmr cassette was inserted in the XbaI-digested pGqseF-updn plasmid to generate pGqseF-updn-Km. The DNA fragment containing the Kmr cassette-disrupted combined upstream and downstream sequences of *qseF* was cleaved by SalI/SphI from pGqseF-updn-Km, and ligated into SalI/SphI-cleaved pUT-Km1 to generate pUTqseF-Km. For *glmY* mutant, pUTglmY-Km was constructed in a similar way except using primer pairs glmY-upF/XbaI-glmY-upR and glmY-dnF/glmY-dnR. For gene inactivation by homologous recombination, pUTqseF-Km or pUTglmY-Km was transferred to wild-type *P. mirabilis* N2 by conjugation. Transconjugants were spread on LSW- plates containing tetracycline (20 µg/mL) and kanamycin (100 µg/mL), and confirmation of mutants with double-crossover events by colony PCR and Southern blot hybridization were performed. For complementation of mutants, the fragments containing full-length *qseF* gene or *glmY* was amplified by PCR using primer pairs qseF-comF/qseF-comR or glmY-comF/glmY-comR, and cloned into pGEM-T easy to generate the plasmid pGEMqseF or pGEM-glmY. *qseF* is driven by *lac* promoter. *glmY* is driven by its own promoter (382 bp), divergent from *lac* promoter, to ensure expression of the correct transcript. pGEMqseF or pGEM-glmY was then transformed into the respective mutant to generate the complemented strain.

#### *4.3. Swarming and Swimming Assays*

The swarming migration assays were performed as described previously [48]. The overnight bacterial cultures (5 µL) were inoculated onto the center of LB swarming plates containing 1.5% (wt/vol) agar, which were then incubated, and the swarming migration distance was measured by monitoring the swarm fronts of the bacterial cells at 1-h intervals. For swimming assays, the overnight culture was stabbed into the center of the swimming plates containing 0.3% agar and migration distance was recorded after incubation for 16 h at 37 ◦C.

#### *4.4. Measurement of the Haemolysin Activity and Cell Length*

The overnight bacterial cultures (120 µL) were inoculated onto the surface of LB swarming plates, which were then incubated at 37 ◦C for 5 h. The hemolysin activity and cell length were determined as described previously [48]. Cell-associated haemolytic activity was determined by incubation of a 20-µL cell suspension (OD600nm = 0.5) in a 980-µL solution of 0.85% NaCl, 20 mM CaCl2 and 2% washed sheep erythrocytes at 42 ◦C for 15 min. After centrifugation, the amount of haemoglobin released by lysis was measured by the increase of the optical density of the assay supernatant at 543 nm.

#### *4.5. Measurement of the Flagellin Level*

Flagellin levels were determined as described previously by SDS-PAGE and Coomassie brilliant blue staining [15]. Western blotting was performed to confirm the flagellin band. The flagellin samples on the SDS-PAGE gel were transferred to an Hybond-P membrane (GE Healthcare, Chicago, IL, USA). The blot was incubated with mouse polyclonal antiserum against FlaA, followed by sheep anti-mouse IgG conjugated with horseradish peroxidase (GE Healthcare, Chicago, IL, USA), and then developed using enhanced chemiluminescence detection reagents (PerkinElmer, Waltham, MA, USA).

#### *4.6. Transmission Electron Microscopy*

Transmission electron microscopy (TEM) was performed as previously described [3] by using 1% phosphotungstic acid (PTA)-stained bacteria on a carbon-coated grid and TEM pictures were obtained with a Hitachi H-7100 electron microscope (Hitachi High-Tech America, Pleasanton, CA, USA).

#### *4.7. Real Time Reverse Transcription PCR (RT-PCR)*

To study the effect of *qseF* or *glmY* deletion on the mRNA amount of motility-related genes, overnight LB cultures (100 µL) of the wild-type, mutants and the complemented strains were spread on the LB agar and incubated for 5 h at 37 ◦C. Total RNA was extracted, and real time RT-PCR was carried out as described previously [3] to measure the mRNA level using primer pairs listed in Table S2. The levels of RNAs were normalized against housekeeping gene *gyrB* mRNA.

#### *4.8. Transcriptional and Translational Reporter Assays*

For transcriptional reporter assay, the promoter region of the gene was amplified by SphI and PstI-included primers and cloned into pGEM-T Easy. These promoter-containing plasmids were cut by SphI and PstI, and the promoter-containing fragment was ligated, respectively to the *xylE*-containing pACYC184-xylE digested by SphI and PstI to construct the transcriptional reporter plasmid. The XylE activity of the transcriptional reporter plasmid-transformed wild-type and mutants was measured as described previously [3]. For translational reporter assay, the 50UTR (520, 253, 138 or 61 bp) with 27 bp coding region (50UTR-27) of *cheA* gene was amplified by SacI and BglII-included primers and cloned into pGEM-T Easy to generate the plasmid pGcheA-50UTR. The fragment containing fulllength *xylE* gene was amplified by PCR using primer pairs BglII-xylE-F and xylE-R, and cloned into pGEM-T Easy to generate the plasmid pGxylE. pGcheA-50UTR was digested with SacI and BglII, and the *cheA* 5 0UTR-27 sequence-containing fragment was ligated to SacI/BglII-digested pGxylE to produce the translational reporter plasmid pcheA of in-frame *cheA-xylE* fusion driven by *lac* promoter. The wild-type and mutants transformed with the translational reporter plasmid were grown overnight in LB broth containing ampicillin (100 µg/mL). Then, the cultures (100 µL) were spread onto the LB agar plate and incubated for 5 h at 37 ◦C before the XylE activity was measured.

#### *4.9. Site-Directed Mutagenesis*

The translational reporter plasmid pGEM-*cheA*-*xylE* (pcheA) and the *glmY*-harboring pBAD plasmid (pBglmY)) containing mutations in the putative GlmY- and *cheA* 5 0UTRinteracting site, respectively, were generated (pcheAm and pBglmYm) by a KOD-Plusmutagenesis kit (Toyobo, Osaka, Japan) using primers listed in Table S2 according to the manufacturer's protocol. The inverse PCR products of pcheA and pBglmY were digested by DpnI to remove the template plasmid DNA, self-ligation of PCR products was performed using T4 polynucleotide kinase and ligase and then DNA sequencing was performed to confirm the DNA sequence of the mutated sites. We introduced pcheA or pcheAm into

the wild-type and compared the XylE activity to confirm the essential GlmY-interacting site in *cheA* 5 0UTR. Additionally, we transformed the combinations of plasmids pBglmYpcheA, pBglmY-pcheAm, pBglmYm-pcheA or pBglmYm-pcheAm into the *glmY* mutant and determined the XylE activity to evaluate the direct interaction of *cheA* 5 0UTR and GlmY.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijms23010487/s1.

**Author Contributions:** Conceptualization, S.-J.L.; investigation, W.-Y.L., Y.-L.T., P.-Y.S. and T.-S.L.; formal analysis, W.-Y.L., Y.-L.T. and P.-H.Y.; resources, Y.-J.L.; data curation, W.-Y.L., Y.-L.T., Y.-J.L., P.-H.Y. and S.-J.L.; supervision, S.-J.L.; writing—original draft preparation, S.-J.L. and W.-Y.L.; writing review and editing, S.-J.L.; project administration, S.-J.L.; funding acquisition, S.-J.L. and Y.-J.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Science and Technology (MOST), Taiwan (grant number: MOST 105-2320-B-002-050-MY3 to S.-J.L.).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We thank the Graduate Institute of Anatomy and Cell Biology, College of Medicine, National Taiwan University for providing assistance in TEM examination.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analyses or interpretation of data, in the writing of the manuscript or in the decision to publish the results.

#### **References**

## *Article* **Ultrasensitive Detection of** *Bacillus anthracis* **by Real-Time PCR Targeting a Polymorphism in Multi-Copy 16S rRNA Genes and Their Transcripts**

**Peter Braun, Martin Duy-Thanh Nguyen, Mathias C. Walter and Gregor Grass \***

Bundeswehr Institute of Microbiology (IMB), 80937 Munich, Germany<EMAIL_ADDRESS>(P.B<EMAIL_ADDRESS>(M.D.-T.N<EMAIL_ADDRESS>(M.C.W.) **\*** Correspondence<EMAIL_ADDRESS>Tel.: +49-992692-3981

**Abstract:** The anthrax pathogen *Bacillus anthracis* poses a significant threat to human health. Identification of *B. anthracis* is challenging because of the bacterium's close genetic relationship to other *Bacillus cereus* group species. Thus, molecular detection is founded on species-specific PCR targeting single-copy genes. Here, we validated a previously recognized multi-copy target, a species-specific single nucleotide polymorphism (SNP) present in 2–5 copies in every *B. anthracis* genome analyzed. For this, a hydrolysis probe-based real-time PCR assay was developed and rigorously tested. The assay was specific as only *B. anthracis* DNA yielded positive results, was linear over 9 log10 units, and was sensitive with a limit of detection (LoD) of 2.9 copies/reaction. Though not exhibiting a lower LoD than established single-copy PCR targets (*dhp61* or *PL3*), the higher copy number of the *B. anthracis*–specific 16S rRNA gene alleles afforded ≤2 unit lower threshold (Ct) values. To push the detection limit even further, the assay was adapted for reverse transcription PCR on 16S rRNA transcripts. This RT-PCR assay was also linear over 9 log10 units and was sensitive with an LoD of 6.3 copies/reaction. In a dilution series of experiments, the 16S RT-PCR assay achieved a thousand-fold higher sensitivity than the DNA-targeting assays. For molecular diagnostics, we recommend a real-time RT-PCR assay variant in which both DNA and RNA serve as templates (thus, no requirement for DNase treatment). This can at least provide results equaling the DNA-based implementation if no RNA is present but is superior even at the lowest residual rRNA concentrations.

**Keywords:** anthrax; *Bacillus anthracis*; 16S rRNA; detection; identification; real-time PCR; RT-PCR

## **1. Introduction**

Within the genus *Bacillus*, the notorious anthrax pathogen *Bacillus anthracis* poses the greatest risk for humans, mammal livestock, and wildlife [1]. Other *Bacillus* spp. such *as B. cereus* or *B. thuringiensis*, which are typical soil bacteria, may also have pathogenic traits related to food poisoning, infections in immunocompromised persons, or production of insecticides [2]. Yet, only obligatory pathogenic *B. anthracis* (and a few *B. anthracis*-like bacilli) features a unique suite of pathogenicity factors rendering the endospore-forming bacterium a first-rate biothreat agent. These factors are encoded on two plasmids called pXO1 and pXO2. Plasmid pXO1 encodes the anthrax toxin genes producing the lethal toxin (gene products of *pagA* and *lef*) and edema toxin (gene products of *pagA* and *cya*) [1]. These toxins damage host cells on various levels [3]. Plasmid pXO2 harbors the capsule genes endowing the pathogen with a poly-glutamyl capsule which helps evade host immune response [1,4]. Phylogenetically, *B. anthracis* belongs to the very closely related *Bacillus cereus sensu lato* group. Besides the better-known species *B. cereus sensu stricto*, *B. anthracis*, or *B. thuringiensis*, the group also comprises several other familiar species such as *B. weihenstephanensis*, *B. mycoides*, *B. cytotoxicus*, and a variety of lesser-characterized members [5].

**Citation:** Braun, P.; Nguyen, M.D.-T.; Walter, M.C.; Grass, G. Ultrasensitive Detection of *Bacillus anthracis* by Real-Time PCR Targeting a Polymorphism in Multi-Copy 16S rRNA Genes and Their Transcripts. *Int. J. Mol. Sci.* **2021**, *22*, 12224. https://doi.org/10.3390/ ijms222212224

Academic Editors: Amelia Casamassimi, Alfredo Ciccodicola and Monica Rienzo

Received: 14 October 2021 Accepted: 10 November 2021 Published: 12 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In the past, the high degree of genetic relatedness to several *B. cereus s.l.* strains has rendered molecular diagnostics of *B. anthracis* challenging (e.g., by polymerase chain reaction assays, PCR). One would think it should be straightforward to identify *B. anthracis* by detecting genetic marker genes (typically *pagA*, *lef*, *cya, capB*, or *capC*) [6–8] on one or both of its virulence plasmids. Identifying these genes comprising constituents of toxin or capsule biosynthesis (*cap*-genes), however, only verifies the presence of these plasmids. This is relevant because several *B. cereus s.l.* isolates are documented to possess very similar virulence plasmids, but not necessarily all of these belong to the species *B. anthracis*. Further, there are *B. anthracis* strains that lack one or both virulence plasmids. Species-specific molecular identification of *B. anthracis* is achieved by targeting a small number of validated chromosomal targets. These targets comprise sections of genes such as *dhp61* (*BA\_5345*; [9]), *PL3* (*BA\_5358*; [6]), or mutations characterized as single nucleotide polymorphisms (SNPs), e.g., in the *rpoB* [7] or the *plcR* [10] gene. A comprehensive overview of suitable and less ideal specific markers for *B. anthracis* has been provided previously [11]. Notwithstanding, the advantage of assaying for pXO1 or pXO2 markers over chromosomal ones is that the plasmid markers occur as multi-copy genes (since the virulence plasmids are present in more than one copy per cell) [12]. Large-scale genomic sequencing revealed that in *B. anthracis* plasmids, pXO1 and pXO2 (with their respective PCR-marker genes) are present on average in 3.86 and 2.29 copies, respectively [13]. Conversely, no multi-copy chromosomal marker has been employed for *B. anthracis* detection thus far.

Likewise, ribosomal RNA (particularly 16S rRNA) has not yet been routinely used for identification and detection of *B. anthracis* even though rRNA molecules are generally the most abundant ribonucleic acid entities in cells constituting up to approximately 80% of total RNA [14]. In fact, copies of 16S rRNA transcripts per cell as constituents of ribosomes number in many thousands (e.g., in *E. coli*, the number of ribosomes per cell ranges from 8 × 103 at a doubling time of 100 min to 7.3 × 104 at a doubling time of 20 min) [15]. Even in stationary culture, a single *E. coli* bacterium contains about 6.5 × 103 copies of ribosomes [16]. Phylogenetically closer to *B. anthracis* than *E. coli* is *Bacillus licheniformis.* For this bacillus, the average number of ribosomes per cell was calculated at 1.25 × 104 , 3.44 × 104 , or 9.2 × 104 in cultures growing at 37 ◦C with doubling times of 120, 60, and 35 min, respectively [17]. While these numbers are well in agreement, somewhat lower numbers of 9 × 103 ribosomes have been determined for exponentially growing cells of *Bacillus subtilis* [18]. While unexplored for *B. anthracis*, bacterial detection using rRNA genes and transcripts has been successfully harnessed to challenge previous limits of detection (LoD) for other pathogens [19–21].

In this study, we introduce a species-specific multi-copy chromosomal PCR marker of *B. anthracis*. This marker is represented by a unique SNP within a variable number of loci of the multi-copy 16S rRNA gene in this organism. Though the 16S rRNA gene sequences feature a very high degree of identity among the *B. cereus s.l*. group species [22], this SNP has previously been identified as unique and present in all publicly available *B. anthracis* genomic data [23–25]. Since all 16S rRNA gene copies harboring the SNP have 100% sequence identity, this specific sequence variation represents a distinct 16S rRNA gene allele named 16S-BA-allele. For simplification, all other 16S rRNA gene alleles lacking the sequence variation were named 16S-BC-allele. The relative abundance of these 16S-BAand -BC-alleles were recently quantified in 959 *B. anthracis* isolates [25]. Here, we also harnessed this SNP to develop a *B. anthracis* specific reverse transcription (RT) real-time PCR assay. This approach brings the multi-copy marker concept for *B. anthracis* up to a new level owing to the excess numbers of ribosomes (and thus 16S rRNA moieties) in relation to chromosomes within a *B. anthracis* cell.

#### **2. Results**

#### *2.1. Set-Up and Optimization of a New 16S rRNA Gene Allele-Specific PCR Assay*

The "16S SNP BA probe" for hybridization to the *B. anthracis* specific sequence variation in 16S-BA-alleles in the *B. anthracis* genome was designed so that the SNP position was

located centrally. In order to increase the fidelity of this probe, six locked nucleic acid (LNA) bases were introduced (Table 1). Similarly, five LNA positions were added to the alternative "16S SNP BC probe", recognizing the non-*B. anthracis* specific 16S-BC-alleles of *B. anthracis* (Table 1). The 16S SNP BA probe was verified in silico against the NCBI database to be highly specific for *B. anthracis*, i.e., all *B. anthracis* genomes showed a 100% match, and only genomes of a few other bacterial isolates exhibited identical sequences. Among these was, e.g., a small number of *Sphingomonas* spp. Others, such as a few genomes annotated as *Staphylococus aureus*, had the same one-base-pair mismatch at the SNP-position (relative to *B. anthracis*) and were thus identical to other *B. cereus s.l.* genomes, hybridizing perfectly against the alternative "16S SNP BC probe" (Figure S1).

**Table 1.** Primers and probes.

Locked nucleic acids are designated by prepositioned (+); 6FAM—6-Fluorescein phosphoramidite; HEX— Hexachloro-fluorescein; BHQ1—Black Hole Quencher-1. \* The expected amplicon length of the PCR reaction is 57 bp. \*\* blocked with a C3-spacer in 30 -position. Hairpin Tm: Primer 16S SNP R: 37.7, else: none; self-dimer Tm: Primer 16S SNP F: 11.5, else: none.

Initially, the 16S SNP BC probe, which deviates only by the one central SNP base from the 16S SNP BA probe, also carried a fluorescent dye/quencher pair. However, since this probe was found to be not entirely specific for recognizing 16S rRNA fragments of *B. cereus s.l.* members, we decided to additionally design this SNP-competing probe as a fluorescently "dark" probe in order to reduce costs of synthesis (Table 1). Thus, the 16S rRNA SNP-PCR may be considered a pseudo-duplex assay (see below for details). All PCR runs were performed with both probes, typically with the 6FAM-labeled 16S SNP BA probe and the dark 16S SNP BC probe.

In silico analysis against the NCBI nt database confirmed that the PCR amplification primers 16S SNP F and 16S SNP R (Table 1) were not species-specific for *B. anthracis*. Indeed, besides DNA from other members of the *B. cereus s.l.* group, these primers would also amplify genome sequences of various other bacteria, such as *Paenibacillus* spp., or the reverse primer would bind to sequences of *Alkalihalobacillus clausii* or *Bacillus licheniformis*, among others. This ambiguity is not surprising for primers hybridizing against 16S rRNA gene sequences. Conversely, the pivotal factor for the detection assay introduced here is that only the 16S SNP BA probe hybridizes without any mismatch against 16S-BA-allele in *B. anthracis* (Figure S1). Thus, the specificity of the PCR assay is uniquely and entirely governed by the LNA-enhanced 16S SNP BA probe.

The 16S rRNA SNP-PCR was robust for deviations from the optimum annealing temperature (62 ◦C; Table S1). Additionally, primer (Table S2), probe (Table S3), and MgCl2 (Table S4) concentrations and pipetting errors (Table S5) were tolerated quite well. Intra- and inter-assay (Tables S6 and S7) variability was determined with positive, weakly positive, and negative template DNA. The average PCR variations were at 0.0–1.1% (intraassay) and 1.1–1.2% (inter-assay), respectively (Tables S6 and S7), indicating high precision of the PCR. Melt point analysis of the 16S-BA-allele PCR product vs. the 16S-BC-allele PCR product (Figure S2) indicated specific amplification of each allele fragment.

#### *2.2. Competitive Amplification-Inhibition of the 16S-BA-Allele Fragment-PCR by Excess of the Alternative 16S-BC-Allele*

Though the new 16S rRNA SNP-PCR assay was tested very robust and precise, we wondered to which degree the assay would be inhibited by large excesses of the alternative 16S-BC-allele fragment featuring a single mismatch at the SNP located centrally

in the hybridizing 16S-BA-allele-specific PCR probe (Figure S1). For testing this, we first evaluated which probe ratio (16S rRNA SNP BA vs. BC probe) would yield the lowest residual fluorescence values (in the 6FAM-channel of the 16S-BA-allele-specific probe) when providing only 16S-BC-allele containing DNA as PCR template. In these tests, the concentration of the 16S-BA-allele-specific probe was kept constant at 0.25 µM. The resulting 6FAM-fluorescence values were very low compared to regular amplification (Table S8), signals were weakly linearly increasing, and no Ct values were detected. The lowest fluorescence, barely above the negative control level, was recorded at a ratio of 0.25/0.75 µM (16S rRNA SNP-BA probe/-BC probe). Thus, this ratio was used for all following tests.

Next, a constant 100 template copies of the 16S-BA-allele fragment per reaction were titrated against increasing copy numbers of the alternative 16S-BC–allele fragment. Figure S3 and Table S9 show that an excess of 16S-BC-allele to BA-allele fragments of 106 , 105 , 104 , or 103 to 1 (Table S9; assay #1–4) inhibits detection of the 16S-BA-allele fragment. This is because there was neither any *bona fide* sigmoidal PCR amplification nor were there any fluorescence signals with values meaningfully above the 16S-BC-allele-only controls (assays #11 and #12). Starting with 7.5 × 104 copies of competing 16S-BC-alleles (vs. 100 16S-BA-allele copies, i.e., 750 to 1; assay #5), both a regular Ct value was provided, and fluorescence started to markedly increase above the base level. At a ratio of 500 to 1 (16S-BCto BA-alleles), *B. anthracis* detection became possible (assays #6 vs. #12; #7). Latest at a surplus of equal or less than 100 to 1 (assay #8), detection of 16S-BA-allele among BC-alleles was robustly possible. Thus, at the very least, a single copy of 16S-BA-allele can be detected in the presence of 100 BC-alleles.

#### *2.3. Sensitivity and Specificity of the 16S rRNA SNP-PCR Assay*

Similar to earlier work [26], we sought to harness the specificity of SNP-interrogation without assaying the alternative SNP state (i.e., the 16S-BC-allele here). Because detecting the 16S-BC-allele was not of interest for the assay at hand, the respective labeled 16S SNP BC probe was replaced by an unlabeled, fluorescently "dark" probe (i.e., a BA allele SNPcompetitor probe; Table 1). In effect, primers would still amplify both alleles; however, the fluorescent probe for the 16S-BA-allele would be outcompeted by the dark probe on 16S-BC-allele targets, and the fluorescent 16S rRNA BA SNP probe would only generate signals in the presence of cognate 16S-BA-allele sequences. Thus, this approach using a dark competing probe would diminish the inadvertent generation of unspecific fluorescence generated by mishybridization of 16S rRNA BA SNP probes to 16S-BC-allele sequences.

To formally validate the sensitivity of the 16S rRNA SNP-PCR assay, a panel of 14 different *B. anthracis* DNAs was employed. These *B. anthracis* strains represent all major branches A, B, and C [27], including prominent sub-branches [28] of the global *B. anthracis* phylogeny (Table S10). All DNAs produced positive PCR results. Similarly, we tested a "specificity panel" of potentially cross-reacting organisms (Table S11). This panel included 13 DNAs of non-*anthracis B. cereus s.l.* strains. Additionally included were DNAs of common animal host organisms such as cattle, goat, sheep, and human. Neither of these DNAs yielded any positive PCR results. Finally, DNAs of organisms relevant for differential diagnostics and other prominent microbial pathogens were also assayed by the new *B. anthracis* specific 16S rRNA SNP-PCR (Table S12). Again, none of these DNAs resulted in false-positive PCR results. Of note, *Sphingomonas zeae* JM-791 [29] harboring 16S rRNA genes 100% identical in the region of the 16S SNP BA probe but different in the primer binding sites yielded negative PCR results. These results clearly indicated that the new PCR is both sensitive and specific for *B. anthracis*.

### *2.4. Linear Dynamic Range, Efficiency, and Limit of Detection of the B. anthracis Specific 16S rRNA SNP-PCR Assay*

The linear dynamic range of the new PCR was determined based on measurements of serial DNA dilutions using recombinant 16S-BA-allele fragments or genomic DNA of *B. anthracis* Ames, respectively, as templates (Figure 1). Linearity was observed over a

range from 101 to 109 copies per reaction for cloned template DNA (Figure 1A; Table S13). In nine out of nine PCR replicates, positive signals were obtained down to 101 copies per reaction. At 100 , two out of nine reactions were negative, thus defining the lower limit of the linear dynamic range. The coefficient of determination (R2 ) was calculated as >0.999. From the slope of the linear regression, the efficacy of the PCR was derived as 2.0 (which is 100.1% of the theoretical optimum). Thus, the 16S rRNA SNP-PCR assay performed very well over a wide 9 log10 concentration range of template DNA. range from 101 to 109 copies per reaction for cloned template DNA (Figure 1A; Table S13). In nine out of nine PCR replicates, positive signals were obtained down to 101 copies per reaction. At 100, two out of nine reactions were negative, thus defining the lower limit of the linear dynamic range. The coefficient of determination (R2) was calculated as >0.999. From the slope of the linear regression, the efficacy of the PCR was derived as 2.0 (which is 100.1% of the theoretical optimum). Thus, the 16S rRNA SNP-PCR assay performed very well over a wide 9 log10 concentration range of template DNA.

*Int. J. Mol. Sci.* **2021**, *22*, x FOR PEER REVIEW 5 of 17

PCR is both sensitive and specific for *B. anthracis*.

*rRNA SNP-PCR Assay*

binding sites yielded negative PCR results. These results clearly indicated that the new

*2.4. Linear Dynamic Range, Efficiency, and Limit of Detection of the B. anthracis Specific 16S*

The linear dynamic range of the new PCR was determined based on measurements of serial DNA dilutions using recombinant 16S-BA-allele fragments or genomic DNA of *B. anthracis* Ames, respectively, as templates (Figure 1). Linearity was observed over a

**Figure 1.** Linearity of the 16S rRNA SNP-PCR. Serial dilutions of DNA of (**A**) a fragment comprising the 16S-BA-allele or (**B**) *B. anthracis* strain Ames were serially diluted 1:10, PCR-tested, and template copies (**A**) or genome equivalents (**B**) plotted against Ct values. Indicated in the graphs are the slopes of the linear regressions and the coefficients of determination (R2 ). Individual data points represent average values from *n* = 3 × 3 PCR-tests.

The *B. anthracis* Ames genome harbors four copies of the 16S-BA-allele and seven copies of the BC-allele. Linear range parameters were very similar to that of cloned 16S-BAallele DNA fragment (Figure 1B; Table S13). Because of the upper concentration limit of our *B. anthracis* Ames DNA preparations, the highest value in the linear range was 106 genome copies. Thus, here the linear range covered target concentrations from 100 to 106 copies per reaction. R2 was determined as >0.999 and the efficacy of the PCR as 1.99 (which is 98.7% of the theoretical optimum). This indicated that the 16S rRNA SNP PCR assay yielded very similar results in these experiments, whether recombinant target DNA or authentic *B. anthracis* DNA was used as templates. Note, though, a single *B. anthracis* Ames genome carries four copies of the 16S-BA-allele. This explains why all PCRs yielded positive signals with DNA template at 100 copies (genome equivalents), whereas PCRs using single copy recombinant template did not. The *B. anthracis* Ames genome harbors four copies of the 16S-BA-allele and seven copies of the BC-allele. Linear range parameters were very similar to that of cloned 16S-BA-allele DNA fragment (Figure 1B; Table S13). Because of the upper concentration limit of our *B. anthracis* Ames DNA preparations, the highest value in the linear range was 106 genome copies. Thus, here the linear range covered target concentrations from 100 to 106 copies per reaction. R2 was determined as >0.999 and the efficacy of the PCR as 1.99 (which is 98.7% of the theoretical optimum). This indicated that the 16S rRNA SNP PCR assay yielded very similar results in these experiments, whether recombinant target DNA or authentic *B. anthracis* DNA was used as templates. Note, though, a single *B. anthracis* Ames genome carries four copies of the 16S-BA-allele. This explains why all PCRs yielded positive signals with DNA template at 100 copies (genome equivalents), whereas PCRs using single copy recombinant template did not.

*Int. J. Mol. Sci.* **2021**, *22*, x FOR PEER REVIEW 6 of 17

**Figure 1.** Linearity of the 16S rRNA SNP-PCR. Serial dilutions of DNA of (**A**) a fragment comprising the 16S-BA-allele or

nation (R2). Individual data points represent average values from *n* = 3 × 3 PCR-tests.

Next, we determined the LoD for the 16S rRNA SNP-PCR assay by probit analysis (Figure 2; numerical data in Table S14). The assay had a limit of detection of 2.9 copies per reaction. This calculates to about 0.6 copies/µL with a probability of success of 95% with a confidence interval of 2.4–4.5 copies/assay. Next, we determined the LoD for the 16S rRNA SNP-PCR assay by probit analysis (Figure 2; numerical data in Table S14). The assay had a limit of detection of 2.9 copies per reaction. This calculates to about 0.6 copies/µL with a probability of success of 95% with a confidence interval of 2.4–4.5 copies/assay.

**Figure 2.** Limit of detection (LoD) of the 16S rRNA SNP-PCR (analytical sensitivity). DNA fragments comprising the 16S-BA-allele were diluted to the indicated copies per reaction (numerical data in Table S14) and subjected to real-time PCR (12 replicates for each data point). Probit analysis (plot of fitted model) was performed to determine the LoD by fitting template copies against the cumulative fractions of positive PCR observations (blue squares and line) and used for calculating the lower and upper 95% confidence limits (red lines). **Figure 2.** Limit of detection (LoD) of the 16S rRNA SNP-PCR (analytical sensitivity). DNA fragments comprising the 16S-BA-allele were diluted to the indicated copies per reaction (numerical data inTable S14) and subjected to real-time PCR (12 replicates for each data point). Probit analysis (plot of fitted model) was performed to determine the LoD by fitting template copies against the cumulative fractions of positive PCR observations (blue squares and line) and used for calculating the lower and upper 95% confidence limits (red lines).

#### *2.5. Comparison of the New 16S rRNA SNP-PCR Assay with Existing PCR Assays 2.5. Comparison of the New 16S rRNA SNP-PCR Assay with Existing PCR Assays*

In order to further assess the performance of the 16S rRNA SNP-PCR assay, we compared it with other established PCR assays for *B. anthracis* identification currently used in our laboratory. These assays target the single-copy genes *dhp61* [9] or *PL3* [6] that have been individually validated before and compared to other commonly used *B. anthracis* PCRs [11]. Using log10 dilutions of *B. anthracis* Ames DNA, the 16S rRNA SNP-PCR exhibited markedly, at least three units, lower Ct values (27.9 ± 0.4; 31.7 ± 0.1; 35.4 ± 0.7) than *dhp61* (32.1 ± 0.0; 35.4 ± 0.6; 38.9 ± 1.5) or *PL3* (31.8 ± 0.2; 36.1 ± 0.7; >40) at 1000, 100, or 10 genome equivalents, respectively (Figure 3). *B. cereus* DNA did not result in amplification In order to further assess the performance of the 16S rRNA SNP-PCR assay, we compared it with other established PCR assays for *B. anthracis* identification currently used in our laboratory. These assays target the single-copy genes *dhp61* [9] or *PL3* [6]that have been individually validated before and compared to other commonly used *B. anthracis* PCRs [11]. Using log10 dilutions of *B. anthracis* Ames DNA, the 16S rRNA SNP-PCR exhibited markedly, at least three units, lower Ct values (27.9 ± 0.4; 31.7 ± 0.1; 35.4 ± 0.7) than *dhp61* (32.1 ± 0.0; 35.4 ± 0.6; 38.9 ± 1.5) or *PL3* (31.8 ± 0.2; 36.1 ± 0.7; >40) at 1000, 100, or 10 genome equivalents, respectively (Figure 3). *B. cereus* DNA did not result in amplification by any PCR assay. This result strongly suggested that the multi-copy 16S rRNA SNP-PCR assay performs competitively when compared back-to-back with established PCR assays for the detection of *B. anthracis*.

**Figure 3.** Comparison of the new 16S rRNA SNP-PCR assay with existing PCR assays. Different 2 quantities of *B. anthracis* Ames template DNA (1000, 100, or 10 genome equivalents per reaction), nontarget DNA (105 templates of *B. cereus* DNA), or water (negative) control were subjected to real-time PCR using the new 16S rRNA SNP assay (**A**), published *dhp61* gene assay [9] (**B**) or published *PL3* gene assay [6] (**C**). Representative amplification curves (from *n* = 3 with similar results) are shown.

#### *2.6. Challenge of the New 16S rRNA SNP-PCR Assay with Samples from a Ring Trial*

Along this line of reasoning, we next challenged the 16S rRNA SNP-PCR assay with samples from a previous ring trial for *B. anthracis* nucleic acid detection [30]. Again, the test was performed in comparison with the established PCR assays for *B. anthracis* identification, *dhp61* [9], and *PL3* [6]. Each of the assays was able to correctly identify the two positives out of four samples (Figure S4). Similar to evaluating known concentrations (Figure 3), the 16S rRNA SNP-PCR assay performed the best. It yielded the lowest Ct values (Figure S4), about two units lower than that of *dhp61* or *PL3* PCR. The 16S rRNA SNP-PCR assay may thus be ideally suited for this kind of analysis in which low target DNA quantities can be expected.

#### *2.7. Challenge of the New 16S rRNA SNP-PCR Assay with Total DNA from Spiked Soil Samples*

Since the 16S rRNA SNP-PCR assay performed well thus far, even in the presence of *E. coli* and human (Figure S4) or competing *B. cereus* (Figure S3) DNA, we evaluated to what extent the assay would be able to detect target DNA in spiked soil samples. These samples were spiked with cells of *E. coli* and *F. tularensis* and cells or endospores of *B. anthracis* and/or *B. thuringiensis* and were subjected to DNA purification. As above, the 16S rRNA SNP-PCR assay was conducted in comparison with the established PCR assays for *B. anthracis* identification *dhp61* [9] and *PL3* [6]. Figure S5 shows the PCR amplification curves. Samples #1, #2, and #4 were samples spiked with *B. anthracis*; sample #3 only contained *E. coli* and *B. thuringiensis*. Sample #4 had a large excess of *B. thuringiensis* over *B. anthracis* (a factor of 104 ). The 16S rRNA SNP-PCR assay detected *B. anthracis* in samples #1 and #2 but not in #4. Conversely, *dhp61* or *PL3* assays detected all three positive samples. The failure to detect *B. anthracis* by the 16S rRNA SNP-PCR assay in sample #4 is in line with our initial tests using massive excess of *B. cereus* DNA competing with *B. anthracis* detection (Figure S3; Table S9). Notably, the 16S rRNA SNP-PCR exhibited markedly, about three units, lower Ct values (23.6 ± 0.7 or 16.4 ± 0.0) than *dhp61* (26.2 ± 0.1 or 19.9 ± 0.1) or *PL3* (25.7 ± 0.0 or 19.5 ± 0.1) for samples #1 and #2, respectively. This result confirmed our preceding findings that the 16S rRNA SNP-PCR assay can reach a lower detection limit than the established assay as long as there is no large excess of other *B. cereus s.l*. DNA competing for amplification primers.

#### *2.8. The New 16S rRNA SNP-PCR Assay also Functions as an RT-PCR Assay*

We reasoned that the real-time 16S rRNA SNP-PCR assay targeting *B. anthracis* DNA might be converted into an RT-PCR assay targeting RNA in the form of 16S-BA-allele transcripts that harbor the *B. anthracis*-specific SNP. In order to test this, cells of *B. anthracis* Sterne or *B. cereus* 10987 were grown to exponential growth phase, inactivated, and total nucleic acids (including genomic DNA) were isolated alongside parallel preparations of DNA only. The one-step RT-PCR reaction was thus run with a mixture of genomic DNA and RNA, which can both be targeted by the assay. For comparison, the above-validated 16S rRNA real-time SNP-PCR was conducted in parallel with genomic DNA as the only template (no RT-reaction). When using identical samples, RT-PCR reactions (with templates consisting of total RNA and DNA) resulted in intensely lower Ct values than without reverse transcription (since only genomic DNA served as a template; Figure 4). Notably, differences in Ct values (RT-PCR vs. PCR) were in the range between 9 and 10 units. This translates to an about 1000-fold improvement using RT-PCR over DNA-only PCR. This result indicated that the 16S rRNA SNP-PCR assay functions both for DNA- and RNA-based (RT) PCR.

**Figure 4.** Comparison of the 16S rRNA SNP-PCR assay (DNA-only) with the RT-16S rRNA SNP-PCR assay (DNA+RNA). Total DNA or total DNA+RNA isolated from exponentially growing cells of *B. anthracis* or *B. cereus*, respectively, were used for PCR amplification of 16S-BA-allele DNA (**A**) or additionally after reverse transcription of 16S rRNA (ribosomal RNA) (**B**). Representative amplification curves (from *n* = 3 with similar results) are shown.

#### *2.9. Linear Dynamic Range, Efficiency, and Limit of Detection of the B. anthracis 16S rRNA SNP RT-PCR Assay*

To further characterize the RT-PCR, we determined the linear dynamic range and determined the LoD (Probit) of the 16S rRNA SNP RT-PCR using total RNA/DNA of *B. anthracis* Sterne (similar to DNA-only templates, see above). The RT-PCR was linear over a range from 100 to 108 template rRNA+DNA per reaction (Figure 5A; Table S15). The coefficient of determination (R2 ) was 0.9982, and the efficacy of the RT-PCR was 1.92 (which is 92.3% of the theoretical optimum). Thus, the 16S rRNA SNP RT-PCR assay performed well over a wide 9 log10 concentration range of template RNA+DNA (higher template numbers than 1.5 × 108 were not tested).

1 The LoD for the 16S rRNA SNP RT-PCR assay as determined by probit analysis (Figure 5B; numerical data in Table S16) was 6.3 copies per reaction. This calculates to about 1.3 copies/µL with a probability of success of 95% with a confidence interval of 5.0–8.9 copies/assay. Thus, the RT-PCR reaction performed similarly well as the PCR reaction. Mindful of the about 3 log10 units higher number of 16S rRNAs in cells than genomes, detection of *B. anthracis* with the rRNA-directed RT-PCR is superior to the respective real-time PCR assay and all other *B. anthracis* PCR assays tested.

**Figure 5.** Linearity and LoD of the 16S rRNA SNP RT-PCR. Serial dilutions of RNA (with DNA) of *B. anthracis* strain Sterne were serially diluted 1:10, RT-PCR-tested and template copies plotted against Ct values (**A**). Indicated in the graph is the slope of the linear regression and the coefficients of determination (R2). Individual data points represent average values from *n* = 3 × 3 PCR-tests. Analytical sensitivity of the 16S rRNA SNP RT-PCR was determined by diluting samples from (**A**) to the indicated copies per reaction (numerical data in Table S16) and subjected to RT-PCR (12 replicates for each data point). To determine the LoD, probit analysis (plot of fitted model, blue squares, and line) was performed (as in Figure 2), and the lower and upper 95% confidence limits (red lines) were determined (**B**). The LoD for the 16S rRNA SNP RT-PCR assay as determined by probit analysis (Figure 5B; numerical data in Table S16) was 6.3 copies per reaction. This calculates to about **Figure 5.** Linearity and LoD of the 16S rRNA SNP RT-PCR. Serial dilutions of RNA (with DNA) of *B. anthracis* strain Sterne were serially diluted 1:10, RT-PCR-tested and template copies plotted against Ct values (**A**). Indicated in the graph is the slope of the linear regression and the coefficients of determination (R2 ). Individual data points represent average values from *n* = 3 × 3 PCR-tests. Analytical sensitivity of the 16S rRNA SNP RT-PCR was determined by diluting samples from (**A**) to the indicated copies per reaction (numerical data in Table S16) and subjected to RT-PCR (12 replicates for each data point). To determine the LoD, probit analysis (plot of fitted model, blue squares, and line) was performed (as in Figure 2), and the lower and upper 95% confidence limits (red lines) were determined (**B**).

#### 1.3 copies/µL with a probability of success of 95% with a confidence interval of 5.0–8.9 **3. Materials and Methods**

#### copies/assay. Thus, the RT-PCR reaction performed similarly well as the PCR reaction. *3.1. Bacterial Culture, Inactivation, and DNA Samples for Quality Assessment*

Mindful of the about 3 log10 units higher number of 16S rRNAs in cells than genomes, detection of *B. anthracis* with the rRNA-directed RT-PCR is superior to the respective realtime PCR assay and all other *B. anthracis* PCR assays tested. **3. Materials and Methods**  *3.1. Bacterial Culture, Inactivation, and DNA Samples for Quality Assessment B. anthracis* strains and other Bacilli were cultivated at 37 °C on tryptic soy agar plates (TSA, Merck KGaA, Darmstadt, Germany). Bacteria comprising the negative panel (Table *B. anthracis* strains and other Bacilli were cultivated at 37 ◦C on tryptic soy agar plates (TSA, Merck KGaA, Darmstadt, Germany). Bacteria comprising the negative panel (Table S1) were grown on appropriate agar media (with 10% CO2 atmosphere where required) at 37 ◦C until colonies emerged. Risk group 3 (RG-3) *B. anthracis* strains were cultivated in the biosafety level 3 (BSL-3) facilities at the Bundeswehr Institute of Microbiology (IMB) and then chemically inactivated by resuspending a loop of colony material in aqueous peracetic acid solution (4% (*v/v*) Terralin PAA, Schülke & Mayr GmbH, Norderstedt, Germany) before further use [31]. RG-2 strains of endospore formers were inactivated

likewise. All other bacterial cultures were inactivated by 70% (*v*/*v*) ethanol. Ring trial *B. anthracis* DNA samples published in [30] were obtained from Instant (Düsseldorf, Germany).

#### *3.2. Isolation of DNA, RNA, and Nucleic Acid Quantification*

Bacterial DNA and RNA were isolated using MasterPure™ Gram Positive DNA Purification kit (Lucigen, Middleton, WI, USA). For RNA (+DNA) isolation, RNase treatment was omitted. DNA and RNA concentrations were quantified using the Qubit dsDNA HS Assay or RNA HS Assay kits (ThermoFisher Scientific, Darmstadt, Germany) according to the manufacturer's protocols. DNA and RNA (+DNA) preparations were stored at −20 ◦C and −80 ◦C, respectively, until further use.

#### *3.3. Design and in Silico Bioinformatic Analysis of Primer and Probe DNA Sequences*

All relevant DNA sequence data for oligonucleotide design were retrieved from public databases (NCBI). Primer and probe DNA oligonucleotides [25] were designed with Geneious Prime (version 2021.1.1; Biomatters, Auckland, New Zealand). In silico specificity analysis was performed by probing each primer and probing nucleotide sequences against the NCBI nt databases using BLASTN for short input sequences (Primer BLAST) [32]. The two amplification oligonucleotide primers target a consensus region within the 16S rRNA genes on the chromosome of *B. cereus s.l.* species (Table 1), including *B. anthracis*. The two oligonucleotide probes (Table 1) feature the centrally located discriminatory SNP (pos. 1110 in *B. anthracis* strain Ames Ancestor, NC\_007530) [23,24]. These probes thus either match the allele unique for *B. anthracis* (named 16S-BA-allele; with an adenine, A at the SNP position) or the general 16S-BC-allele (guanine, G at the SNP position), respectively (the two alleles are depicted in Figure S1). Due to placement and length restrictions related to another non-discriminatory SNP (pos. 1119), each probe was amended with locked nucleic acids (LNA). LNA are modified nucleic acids in which the sugar is conformationally locked. This rigidity causes exceptional hybridization affinity through stable duplexes with DNA and RNA [33], eventually improving mismatch discrimination in SNP genotyping studies. Similar to unmodified ssDNA probes, the LNA-containing probes (Table 1) are susceptible to 50 -nuclease attack during PCR. LNA probes as well as primers were purchased from TIB MolBiol (Berlin, Germany).

#### *3.4. Real-Time and Reverse Transcription PCR Conditions*

All (pseudo) duplex real-time PCR amplifications were performed in reaction mixtures of a final volume of 20 µL containing 2 µL LightCycler® FastStart DNA Master HybProbe mix (Roche Diagnostics, Mannheim, Germany), 5 mM MgCl2, 0.5 µM of each primer, 0.25 µM of 16S SNP BA probe, 0.75 µM of (dark) 16S SNP BC probe, and various quantities of template DNA template. All reactions were performed on a LightCycler 480 real-time PCR system fitted with color compensation (Roche Diagnostics, Mannheim, Germany). The optimized amplification conditions were 95 ◦C for 10 min, and then 45 consecutive cycles of first 15 s at 95 ◦C and then 20 s at 62 ◦C, followed by 20 s at 72 ◦C.

Reverse transcription PCR reaction mixtures contained 7.4 µL LightCycler® 480 RNA Master Hydrolysis Probes mix, 1.3 µL Activator, 1 µL Enhancer (Roche Diagnostics, Mannheim, Germany), 0.5 µM of each primer, 0.25 µM of 16S SNP BA probe, 0.75 µM of (Dark) 16S SNP BC probe, a variable volume of RNA and/or DNA template. Finally, nuclease-free water (Qiagen, Hilden, Germany) was added to a final volume of 20 µL. Using the LightCycler 480 real-time PCR system (Roche Diagnostics, Mannheim, Germany), reverse transcription was performed at 63 ◦C for 3 min followed by an activation step at 95 ◦C for 30 s and 45 cycles of 95 ◦C for 15 s, 62 ◦C for 20 s and 72 ◦C for 1 s.

A fluorescent signal 10-fold higher than the standard deviation of the mean baseline emission was counted as a positive detection. Samples were tested in triplicate (unless noted otherwise) and data recorded as Cycle thresholds (Ct) with Ct defined as the PCR cycle at which the fluorescent intensity raised above the threshold [34].

#### *3.5. Droplet Digital PCR (ddPCR) and Reverse Transcription (RT) ddPCR*

All DNA and RNA templates used for real-time and reverse transcription PCR were quantified by ddPCR and RT ddPCR, respectively. A 20 µL ddPCR reaction mixture consisted of 10 µL ddPCR Supermix for Probes (Bio-Rad Laboratories, Munich, Germany), 0.9 µM of each primer, 0.15 µM of each probe, and 5 µL of template DNA. RT-ddPCR reaction mixtures comprised of 5 µL One-Step RT-ddPCR Advanced Supermix for Probes (Bio-Rad, Munich, Germany), 2 µL of Reverse Transcriptase (Bio-Rad, Munich, Germany; final concentration 20 U/µL), 0.6 µL of DTT (Bio-Rad, Munich, Germany; final concentration 10 nM), 0.9 µM of each primer, 0.15 µM of each probe, and 5 µL of template RNA. Droplets were generated using a QX200 ddPCR droplet generator (Bio-Rad, Munich, Germany). PCR amplification for both assays was performed on the Mastercycler Gradient (Eppendorf, Hamburg, Germany) with the following conditions.

Initial reverse transcription was carried out at 48 ◦C for 60 min (only for RT-ddPCR). Enzyme activation at 95 ◦C for 10 min was followed by 40 cycles of denaturation at 94 ◦C for 30 s and annealing/extension at 58 ◦C for 1 min. Before the samples were cooled to 4 ◦C, a final enzyme inactivation was carried out at 98 ◦C for 10 min. The cooling and heating ramp rate was set to 2 ◦C/s for all steps. After PCR runs, droplets were analyzed using the QX100 Droplet Reader (Bio-Rad, Munich, Germany), and absolute target concentrations of each sample were calculated using Quantasoft Pro Software (Bio-Rad, Munich, Germany).

### *3.6. Generation of PCR Positive Controls from Reference Plasmids Harboring 16S-BA- or BC-Allele Fragments*

Though we generally used genomic DNA from *B. cereus* or *B. anthracis*, respectively, for PCR testing and validation, generic positive control reference plasmids for either allele, the *B. anthracis*-specific 16S-BA-allele or the *B. cereus*-specific 16S-BC-allele were constructed. For this, a PCR-amplicon was generated from *B. anthracis* Ames DNA with primers 16S SNP F and 16S SNP R using Platinum™ Taq DNA Polymerase High Fidelity (ThermoFisher Scientific, Darmstadt, Germany). This DNA comprises a mixture of both alleles in a ratio of 4 to 7 [25]. The PCR-amplicon was analyzed on agarose gel electrophoresis, a band of the expected size (57 bp) cut from the gel and gel-purified using QIAquick Gel Extraction kit (QIAGEN, Hilden Germany). PCR products were ligated into pCR2.1 TOPO vector (ThermoFisher Scientific, Darmstadt, Germany) using TOPO TA Cloning kit (Thermo Scientific, Darmstadt, Germany) and transformed into One Shot TOP10 chemically competent cells (ThermoFisher Scientific, Darmstadt, Germany) according to the manufacturer's protocol. Several recombinant plasmids isolated from different clones were sequenced (Eurofins Genomics Germany, Ebersberg, Germany) in order to obtain plasmids harboring either the 16S-BA-allele or the 16S-BC-allele. From these plasmids, PCR products were generated using primers M13 F and M13 R, which contained the target region for the 16S rRNA SNP-PCR with either the 16S-BA- or BC-allele. After purification with QIAquick PCR purification kit (QIAGEN, Hilden Germany), PCR products were quantified using digital PCR and diluted as required.

### *3.7. Determination of the Specificity (Inclusivity/Exclusivity) of the B. anthracis 16S rRNA Allele Assay*

PCR specificity for the 16S rRNA SNP assay was assessed by verifying the amplification of DNA containing or lacking respective markers. "Inclusivity" was evaluated by (exponential) amplification above threshold levels obtained with template DNA comprising the markers' sequences. Vice versa, "exclusivity" was confirmed by lack of amplification of genomic DNA from *B. cereus s.l.* strains reported to lack the particular 16S-BA-allele but also may harbor the alternative 16S-BC-allele or include no-template negative controls (NTC). Positive PCR results were further analyzed via agarose gel electrophoresis, demonstrating a single band with a molecular weight corresponding to the predicted size of the 16S rRNA SNP-PCR amplicon (note: this cannot differentiate between the two alternative SNP states in the 16S rRNA gene alleles).

#### *3.8. Dynamic Linear Range, PCR Efficiency, and Limit of Detection*

The dynamic linearity of the PCR assay was determined over a 9 log10 concentration range for DNA (real-time PCR) and RNA (RT-PCR) templates. Each dilution was assayed 6-fold, and analysis for linearity and PCR-efficiency (E) was performed from the plot of the Ct's versus the logarithm of the target concentrations [35]. The sensitivity of the PCR assay was expressed as the limit of detection (LoD) of 16S rRNA SNP genome or transcript copies. LoD was formally defined as the concentration permitting detection of the analyte at least 95% of the time. For this, DNA fragments comprising the 16S rRNA SNP were diluted to between 10 and 0 copies per reaction, subjected to real-time PCR with 12 replicates for each dilution step. Probit analysis (plot of the fitted model) was performed [36] using StatGraphics Centurion XVI.I (16.1.11; Statgraphics Technologies, The Plains, VA, USA) to determine the LoD by fitting template copies against the cumulative fractions of positive PCR observations and used for calculating the lower and upper 95% confidence limits. The LoD of the 16S rRNA SNP RT-PCR was determined likewise using samples with 0–15 rRNA copies per reaction (12 replicates for each dilution step).

#### **4. Discussion**

The use of SNPs as reliable markers for the identification of *B. anthracis* among its closest relatives of the *B. cereus* group is not a novel approach. This has previously been achieved with high specificity and sensitivity for nucleotide position 640 in the *plcR* gene [10] or at position 1050 in the *purA* gene [26], and diverse assays were thoroughly evaluated in [11]. Likewise, ribosomal gene sequences and intergenic transcribed spacers (ITS) between 16S and 23S rRNA genes have also been employed for *B. anthracis* identification in the past [37–40]. However, while these authors focused on the specific identification of *B. anthracis*, they neglected the potential of developing a sensitive assay making use of the multi-copy nature of their targets. An interesting exception is a study on fluorescent DNA-heteroduplex detection of *B. anthracis* [41]. Herein detection was preceded by general PCR-amplification of a fragment of the 16S rRNA gene region of *B. cereus s.l.* group strains containing a presumably specific SNP (pos. 980). This SNP, however, is neither specific for *B. anthracis* nor for the *B. cereus s.l.* group [24]. Anyway, Merrill et al. succeeded in establishing a LoD for their PCR of approximately 0.05 pg of purified *B. anthracis* genomic DNA (which can be calculated to represent 10–20 cell equivalents per reaction) [41]. This is higher than the LoD of about 1–2 cell equivalents per reaction found in our study. More importantly, Merrill et al. also took the effort to determine the detection limit of their presumably specific SNP in mixtures of 16S rRNA gene amplicons from *B. anthracis* and *B. cereus* [41]. The authors observed a detection limit of 1 out of 50 for *B. anthracis* DNA mixed with *B. cereus* DNA. They explained this limit as narrowed by methodological constraints and from competitive hybridization dynamics during probe annealing [41]. This finding can be compared with our results. The PCR assay developed here was able to detect at least one *B. anthracis* 16S-BA-allele target among 100 BC-allele targets (Figure S3 and Table S9). At higher alternative (16S-BC-allele) concentrations, these templates will outcompete the 16S-BA-allele for primer binding. Thus, the higher the fraction of 16S-BC-allele, the lower the relative amplification of 16S-BA-allele resulting in increasingly non-exponential amplification of the latter. In contrast, for a SNP in the DNA target *plcR* used for the differentiation of *B. anthracis* from *B. cereus*, a 20,000-fold excess of the alternate *B. cereus* allele did not preclude the detection of the *B. anthracis* allele [42]. With *B. cereus* spore counts in soils spanning a wide range of from 1 × 101 to 2.5 × 104 CFU per g soil [43], the *plcR* SNP-PCR should be able to detect *B. anthracis* in practically any sample. Here, the new 16S rRNA SNP-PCR on DNA as target molecule would fall short with only covering up to medium *B. cereus*-loaded soils. However, when targeting ribosomal RNA, the sensitivity (LoD) of the 16S rRNA SNP RT-PCR would be at least three orders of magnitude increased. Then, it should be possible to challenge the LoD values achieved by the *plcR* SNP-PCR (25 fg DNA or about 5 genome equivalents) [42].

A potential limitation of the multi-copy nature of the 16S-BA-allele may be the variable abundance of this allele in different *B. anthracis* strains. Previously, we could show that most *B. anthracis* strains harbor 3 (58.39%) 16S-BA-alleles. There are, however, also a number of isolates only possessing 2 (23.04%), 4 (17.10%), 5 (1.15%), and a single one with only 1 (0.31%) 16S-BA-alleles [25]. Thus, in most cases, this multi-copy gene allele can be harnessed nevertheless. A more typical multi-copy marker for detection of bacterial biothreat agents (and of other pathogens) constitute insertion sequence (IS) elements, which are widespread mobile genetic entities. For instance, in *Brucella* spp. *IS711* occurs in multiple genomic copies, and thus, the detection of this *IS711* is very sensitive. *B. melitensis* and *B. suis* contain seven complete copies, *B. abortus* carries six complete and one truncated *IS711* copies, *B. ovis*, *B. ceti*, and *B. pinnipedialis* even more than 20 copies [44]. Consequently, the lowest concentration of *Brucella* sp. DNA that could be detected was about ten times lower for *IS711* than, e.g., for single-copy genes *bcsp31* (*Brucella* cell surface 31 kDa protein) or *per* (perosamine synthetase), respectively [45]. Similarly, in *Coxiella burnetii*, the detection sensitivity of specific *IS1111* was compared to that of the single-copy *icd* gene (isocitrate dehydrogenase) [46]. While both PCRs for *icd* and *IS1111* had similar LoDs of 10.38 and 6.51, respectively, the sensitivity of *IS1111* was still superior because of its multiple-copy nature. Between 7 and 110 copies of this mobile element were found in various *C. burnetii* isolates [46].

The differences in threshold values (∆Ct = 9.96 ± 0.65) of identical samples obtained from (RT)-PCR using 16S-BA-allele DNA-only vs. 16S-BA-alleleDNA+RNA is enormous. There is an approximate factor of about 1000 (29.96) times more templates in the DNA+RNA sample than in the DNA-only sample. This factor favorably agrees with the numbers of genome copies and 16S rRNA transcripts in cells [17,18]. Similar to the work at hand, earlier work employed a combination of a DNA multi-copy marker and sensitive detection of rRNA transcript targets in *Mycobacterium ulcerans* [20]. The authors determined an LoD of six copies of the 16S rRNA transcript target sequence. For comparison, an LoD of two target copies of the high-copy insertion sequence element *IS2404*, which is present in 50–100 copies in different *M. ulcerans* strains, was calculated from parallel experiments [20]. Ribosomal RNA detection was also utilized for *Mycobacterium leprae* diagnosis by the same research team. Here, an LoD of three *M. leprae* target copies was achieved for a novel 16S rRNA RT-PCR assay; the same value as determined for the *M. leprae* specific multi-copy repetitive DNA target assayed in parallel [21]. At first glance, these values do not especially speak in favor of querying for 16S rRNA transcripts; however, one has to consider the high numbers of these molecules per cell in comparison to DNA markers (including the high-copy ones). Thus, the chance of capturing one of the more abundant rRNA molecules should be higher than that of the more limited DNA molecules. Indeed, this idea was explored, e.g., for *Escherichia coli*, *Enterococcus faecalis*, *Staphylococcus aureus*, *Clostridium perfringens*, and *Pseudomonas aeruginosa* by [19]. Comparative quantitative detection of these bacteria by RT-PCR (16S rRNA) vs. PCR (16S rRNA genes) revealed that the rRNAdetecting assay was from 64- to 1024-fold more sensitive than the one detecting DNA. Similarly, work on pathogenic spirochete *Leptospira* spp. found that 16S rRNA-based assays were at least 100-fold more sensitive than a DNA-based approach [47]. These authors also found that Leptospiral 16S rRNA molecules remain appreciably stable in blood. From this insight, the authors then highlighted the potential use of 16S RNA targets for the diagnosis of early infection. Nevertheless, potential limitations of this approach were also noted. Efficacy of the required reverse transcription reaction has to be considered, RNA molecules are notoriously less stable than other biomarkers, and their cellular abundance (and as a consequence, their detection) can be expected to be variable [47]. Finally, though not required for qualitative detection, absolute quantification of microbial cells based solely on enumeration of RNA molecules is complicated because of these variations in transcript numbers depending, e.g., on the growth phase [47]. However, the cell numbers determined by RT-PCR were similar when compared alongside standard methods such as cell counts, PCR, or fluorescence in situ hybridization (FISH) [48]. Yet, in certain instances, there might

be an additional advantage of performing PCR on rRNA directly (via RT-PCR) instead of targeting DNA (including DNA of rRNA genes). Because DNA is more stable than RNA, DNA may originate from both live and dead bacterial cells. In contrast, rRNA molecules may be considered to be more closely associated with viable bacteria [49]. Though this might also be possible with the new PCR assays introduced in the work at hand, we chose to combine DNA and rRNA detection in a single test tube for the sake of simplicity (no troublesome DNase treatment of purified RNA required) and depth of detection.

#### **5. Conclusions**

In this work, we designed and validated a new PCR-based detection assay for the biothreat agent *B. anthracis*. This assay can be run as a real-time PCR with solely DNA as a template or as an RT-real-time version using both cellular nucleic acid pools (DNA and RNA) as a template. This assay was found to be highly species specific, yielding no false positives, and was sensitive with a LoD of about 0.6 copies/µL (DNA-only) and about 1.3 copies/µL (DNA+RNA). With the high abundance of 16S rRNA moieties in cells, this assay can be expected to facilitate the detection of *B. anthracis* by PCR. While standard PCR assays are well established for the identification of *B. anthracis* from pure culture, the exceptional sensitivity of the new 16S rRNA-based assay might excel in clinical and public health laboratories when detection of minute residues of the pathogen is required.

**Supplementary Materials:** The following are available online: https://www.mdpi.com/article/10.3 390/ijms222212224/s1.

**Author Contributions:** Conceptualization, G.G. and P.B.; investigation, P.B., M.D.-T.N. and M.C.W.; methodology, M.D.-T.N. and P.B.; formal analysis and validation, P.B., M.D.-T.N. and G.G.; resources, G.G. and M.C.W.; data curation, P.B., M.D.-T.N. and M.C.W.; writing—original draft preparation, G.G. and P.B.; writing—review and editing, P.B., M.D.-T.N., M.C.W. and G.G.; visualization, M.D.-T.N., P.B. and G.G.; supervision and project administration, G.G. and P.B.; funding acquisition, G.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by funds from the Medical Biological Defense Research Program of the Bundeswehr Joint Medical Service.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors thank Mandy Knüper (Bundeswehr Institute of Microbiology) for support with probit-analysis and the gift of DNA from spiked soil samples, Stefanie Gläser (Justus-Liebig-University Gießen, Germany), Erwin Märtlbauer (Ludwig-Maximilians-University, Munich, Germany) and Monika-Ehling Schultz (University of Veterinary Medicine, Vienna, Austria) for several *B. cereus* strains, as well as Paul Keim (Northern Arizona University, Flagstaff, AZ, USA) and Wolfgang Beyer (Hohenheim University, Stuttgart, Germany) for the gift of *B. anthracis* strains/DNA. Thanks are due to Rahime Terzioglu for technical assistance and Olfert Landt (TIB Molbiol, Berlin) for support in LNA-probe design.

**Conflicts of Interest:** The authors declare no conflict of interest. Opinions, interpretations, conclusions, and recommendations are those of the authors and are not necessarily endorsed by any governmental agency, department, or other institutions. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*International Journal of Molecular Sciences* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/ijms

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-7736-4