Housekeeping gene gyrA, a potential molecular marker for Bacillus ecology study

doi:10.21203/rs.3.rs-1557860/v1

PDF

Research Article

Housekeeping gene gyrA, a potential molecular marker for Bacillus ecology study

https://doi.org/10.21203/rs.3.rs-1557860/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 20 Apr, 2022

You are reading this latest preprint version

Background: Bacillus is a ubiquitous microorganism and is one of the most commercially important species widely used in industry, agriculture and healthcare. Bacillus is relatively well understood at the single-cell level; however, molecular tools that determine diversity and ecology of Bacillus community are limited, which limits our understanding of how the Bacillus community works.

Results: In this study, we investigated the potential of the housekeeping gene gyrA as a molecular marker for determining the diversity of Bacillus species, then designed a new primer pair that specifically target gyrA (gyrA3). The rate of variation of the gyrA gene was higher than 16S rRNA gene in inter- and intraspecific Bacillus. And the amplification efficiency for Bacillus species diversity could be greatly improved by primer design, the novel primer pair gyrA3 that can detect at least 92 Bacillus species and related species and showed its advantage in mock Bacillus community compared to 16S rRNA gene. For Bacillus amyloliquefaciens, Bacillus pumilus, and Bacillus megaterium, the high variability of the gyrA gene allows for more detailed clustering at the subspecies level that cannot be achieved by the 16S rRNA gene.

Conclusion: Collectively, the gyrA gene provides better phylogenetic resolution than 16S rRNA gene and informs on the diversity of the Bacillus community, therefore gyrA gene may have broad application prospects in the study of Bacillus ecology.

Molecular marker

16S rRNA

Housekeeping gyrA gene

SNP

Bacillus diversity

Microbial communities in soil are known to be one of the largest reservoirs of biological diversity and have been extensively studied[1]. Current advances in high-throughput DNA sequencing of a portion of the small subunit of ribosomal RNA (16S and 18S rRNA) form the backbone of most studies of soil microbial ecology[2]. For bacteria, the 16S rRNA gene is usually preferred, because it contains both highly conserved and hypervariable regions[3], and especially because comprehensive reference databases have been compiled for comparison[4–7]. However, 16S rRNA amplicon sequencing also has many shortcomings: first, 16S rRNA evolves slowly and is highly conserved, making it a poor marker for distinguishing between closely related strains. Second, chimera formation during PCR is high because 16S rRNA variability is very low[8, 9]. Third, the number of 16S rRNA copies in different species is highly variable, and single nucleotide polymorphisms (SNPs) at the single-cell level may result in an overestimation of diversity[10]. Fourth, the similarity between species can be very high, making it difficult to delineate species in cluster analysis, and different clustering levels lead to different results[11]. Therefore, new complementary taxonomic markers for genetic and bioinformatic analysis need to be developed to study microbial diversity in more detail, especially at the subspecies level.

Bacillus is one of the most intensely studied bacterial genus comprising at least 200 species[12]. It is a heterogeneous bacterial taxon that is ubiquitous in various ecological niches and widely used in medicine, industry and agriculture. Although there have been many in-depth studies on Bacillus model species, community-level studies on Bacillus in soil and other habitats lag behind. Because sequence of 16S rRNA within Bacillus species is often similar, the definition and delineation of bacterial species based on the 16S rRNA gene comparison among related species in the genus Bacillus is unclear. Therefore, identification and typing of Bacillus isolates based on the 16S rRNA gene alone cannot provide accurate results and it is important to explore and use other genes as molecular markers to assess the diversity of the Bacillus community[12].

Housekeeping genes are potential candidates for assessing microbial diversity because they have been shown to elicit higher phylogenetic resolution than the 16S rRNA gene, such as the rpoB gene, gyrA gene, and gyrB gene, etc.[13–21]. Typically, there are only one or two copies of housekeeping genes per genome, and the use of low-copy number genes compared to high-copy number genes could lead to more accurate diversity analysis by avoiding overestimation of diversity due to SNPs in different gene copies. Indeed, several studies have tested the rpoB gene and the gyrB gene as molecular markers to analyze the diversity of bacterial communities by amplicon sequencing. The results showed that housekeeping gene sequencing provided a more accurate description of bacterial community composition than 16S rRNA sequencing under certain conditions[22–24].

The housekeeping gene gyrA, encoding DNA gyrase subunit A, is essential for DNA replication and present in all bacteria[25]. Analyses of gyrA identity provided higher phylogenetic resolution than the 16S rRNA gene for tested Bacillus isolates[13, 21]. Specifically, partial gyrA gene sequences was used for phylogenetic analysis and species identification of seven Bacillus strains, including B. amyloliquefaciens, B. atrophaeus, B. licheniformis, B. mojavensis, B. subtilis, B. subtilis subsp. Spizizenii, and B. vallismortis[13]. Moreover, the gyrA gene sequences provided a good marker for B. subtilis and B. amyloliquefaciens and showed better discriminatory potential between these two closely related species than the rpoB gene[13]. The gyrA gene has been also successfully used to detect intraspecific diversity of B. subtilis isolates from soil microscale[17] and tomato rhizosphere isolates[26]. Overall, these works suggest that gyrA has a good potential to be used as a molecular marker for microbial ecology studies of Bacillus genus and related species, however, to the best of our knowledge, the available primer pairs used in studies indicated above, have not been used for amplicon-based community analyses of Bacillus species.

In this study, we first compared the rate of variation of the 16S rRNA and gyrA genes between 20 Bacillus species and then designed a new primer pair that specifically target gyrA (gyrA3). We then evaluate their performance by using in silico PCR, testing their efficiency on Bacillus isolates and performed SNP analysis of 16S rRNA and gyrA genes for selected species that were available in the NCBI database. Finally, verified gyrA3 primers to differentiate species and strains of Bacillus mock community, and compared the obtained results with those targeting 16S rRNA. Our results suggest that the gyrA gene is a useful molecular marker for identification of Bacillus isolates and describing the diversity of the Bacillus community.

The gyrA gene of the Bacillus genus shows higher variation rates than 16S rRNA

The housekeeping gene gyrA is considered to be more variable than 16S rRNA and has been used as a molecular tool for classification and identification of Bacillus subtilis species[13, 27]. To determine the specific variation divergence of 16S rRNA and gyrA gene sequences in the genus Bacillus, two alignments with sequences of the entire 16S rRNA and gyrA genes from 20 Bacillus species (361 genomes) (Table S1) were used to analyze nucleotide diversity. The nucleotide diversity (Pi) of 16S rRNA and gyrA gene calculated by the DnaSP software was 0,039 and 0.491, respectively (Fig. 1A,B blue line). The results showed that the degree of interspecies variation was significantly higher for the gyrA gene than for the 16S rRNA gene.

Intraspecific variation analysis of the two genes was performed using the same approach but focusing on three Bacillus species: Bacillus amyloliquefaciens (88 genomes), Bacillus licheniformis (71 genomes) and Bacillus pumilus (97 genomes) (Table S1). The intraspecific nucleotide diversity (Pi) of 16S rRNA was again lower than the intraspecific nucleotide diversity of gyrA gene sequences in all three species: B. amyloliquefaciens (Pi_16S = 0.0014; Pi_gyrA=0.0244), B. licheniformis (Pi_16S = 0.00024; Pi_gyrA=0.0021) and B. pumilus (Pi_16S = 0.00136; Pi_gyrA=0.0344) (Fig. 1A,B nonblue lines). In B. amyloliquefaciens, B. licheniformis and B. pumilus, the degree of variation of the gyrA gene was significantly higher than that of 16S rRNA. In conclusion, the Bacillus gyrA gene shows higher variation rates than 16S rRNA, hence we propose that gyrA represents a promising molecular marker for analyses of Bacillus community diversity analyses and the diversity of Bacillus isolates.

First comparative tests of three primer pairs for the detection of Bacillus species

As indicated above Bacillus isolates have been already analyzed by primers targeting gyrA, however the specificity of these primers has not been investigated broadly[28, 29]. In order to satisfy the amplicon sequencing requirements, we designed a new primer pair (gyrA3) (Fig. 2A), and compared its amplification potential in colony PCR and virtual PCR with the previously designed primers, referred to here as gyrA1[17] and gyrA2[30].

First, we selected seven strains of different Bacillus species: Lysinibacillus fusiformis, Paenibacillus polymyxa, Bacillus pumilus, Bacillus velezensis, Bacillus megaterium, Bacillus cereus and Bacillus subtilis (Fig. 2B) to perform PCR amplification with primer pairs gyrA1, gyrA2 and gyrA3. The PCR amplification results showed that gyrA1 detected only B. subtilis; gyrA2 detected B. subtilis, B. velezensis and L. fusiformis; whereas gyrA3 performed much better and detected all Bacillus species included in the analysis (Fig. 2B and Fig. S1).

For in-silico PCR analysis, we constructed a gyrA gene library containing 5062 full-length gyrA gene sequences from 226 Bacillus species obtained from the National Center for Biotechnology Information (NCBI) (Table S2). Our criteria for positive in-silico gyrA PCR amplification were annealing of the primer pair with at least 18 bases match (forward primer and reverse primer). The results showed that only 8 Bacillus species were amplified in-silico by gyrA1, 9 Bacillus species were amplified by gyrA2 (Fig. S2), while 55 Bacillus species were amplified by gyrA3 (Fig. 2C). The majority of sequences amplified by gyrA1 and gyrA2 belonged to B. subtilis, whereas gyrA3 demonstrated broader diversity as evidenced by the amplification of seven species and in-silico PCR (Fig. 2B,C and Table S3).

Specificity range of the gyrA3 primer pair by using PCR and in silico PCR

Because the gyrA3 primer pair performed better than the previously reported primers, we next combined analysis of the in-silico amplified gyrA genes with PCR analysis of Bacillus isolates from our laboratory culture collection. Virtual gyrA3 PCR amplicons from 55 different Bacillus species from the NCBI database were evenly distributed among the branches of the phylogenetic tree (Fig. 3 orange and green).

Next, we used 127 Bacillus strains and related genera (Paenibacillus, Lysinibacillus, Aneurinibacillus, Virgibacillus, Brevibacillus, Halobacillus, Fictibacillus) from our laboratory culture collection to amplify their gyrA genes with a gyrA3 primer pair (Table S4). The results showed that 28 Bacillus species and 16 Bacillus-related species could be amplified by the gyrA3 primers (Fig. 3 blue and green). These Bacillus species, marked in green in the phylogenetic tree, 7 species were detectable by both methods (Fig. 3 green). In summary, the primer pair gyrA3 can potentially detect 76 Bacillus species and as many as 16 species from related genera (Fig. 3).

The gyrA gene provides better intraspecific phylogenetic resolution than the 16S rRNA gene among certain species

Compared to 16S rRNA, the molecular evolution rate of gyrA gene sequences is faster[1], so we hypothesized that gyrA might provide better phylogenetic resolution at the subspecies level. To analyze the intraspecific variability of the 16S rRNA and the gyrA gene, we downloaded the genomes of four Bacillus species, each of which had more than 100 genomes available in the NCBI database, including B. amyloliquefaciens (116 genomes), B. pumilus (140 genomes), B. megaterium (117 genomes) and B. anthracis (226 genomes) (Table S5). We did not include analysis the of B. subtilis genomes because templates for gyrA primers have already been developed and applied for analyses of this species[17, 28–30].

In B. amyloliquefaciens, the 480 bp long 16S rRNA region (V3-V4 region) contained 12 variable base sites (Fig. 4A and Fig. S3A), which were detected in only 4 of 116 genomes of this species (Fig. 4A), and the SNP frequency at variable sites within the 4 genomes ranged from 0.21%-1.46% (Fig. S3C red column). In contrast, the 500-nucleotide gyrA region (positions 350–850) contained 59 variable sites (Fig. 4B and Fig. S3B). The variable sites were detected in 109 of 116 genomes (Fig. 4B), and the frequency of SNPs at the variable sites ranged from 0.4–6.4% (Fig. S3C blue column).

In B. pumilus, alignment of the 16S rRNA V3-V4 region revealed 21 variable base sites (Fig. 4C and Fig. S4A), but again only in 11 out of 140 genomes (Fig. 4C). The frequency of SNPs at variable sites ranged from 0.21–3.54% (Fig. S4C red column). In contrast, in the gyrA gene 130 of the base sites were variable (Fig. 4D and Fig. S4B) and these were found in 87 of 140 genomes, with SNP frequencies at variable sites ranging from 0.2–15.8% (Fig. S4C blue column). However, in B. pumilus nearly 50% of the genomes examined had 100% identity in the gyrA gene, indicating high degree of relatedness between genomes that may require sequencing of additional marker genes for clonality verification and strain typing.

We also observed that the variation of the gyrA gene is quite different in different Bacillus species, which could be a bias of the NCBI database or a property of certain species. For example, B. megaterium, showed lower diversity with 3 bases of variation in the 16S rRNA alignment region (V3-V4 region) (Fig. S5A and Fig. S6A). Although 48 of 117 genomes showed polymorphism, the maximum SNP frequency at variable sites of B. megaterium genomes was only 0.42% (Fig. S6C red column). The gyrA gene was again more polymorphic with 58 bases of variation (Fig. S5B and Fig. S6B) occurring in 99 of 117 genomes, with SNP’s frequencies at variable sites ranging from 0.2–2.8% (Fig. S6C blue column).

The variation divergence between the two genes of Bacillus anthracis was much lower than in the three Bacillus species described above. Although we identified 13 variable base sites in the 16S rRNA V3-V4 region and 65 variable base sites in the gyrA gene region (Fig. S5C,D and Fig. S7A,B), these single nucleotide polymorphisms (SNPs) occurred in only 7 and 18 of 226 genomes, respectively. Moreover, the SNP frequencies at variable sites in 7 and 18 of B. anthracis genomes were also low: 0.21%-1.04% and 0.2%-7.4%, respectively (Fig. S7C).

Overall, our results showed that within Bacillus species the frequency of SNPs in the gyrA gene was consistently much higher than in the 16S rRNA (Fig. 4 and Fig. S5). We therefore suggest that the gyrA gene provides better resolution than 16S rRNA for identification and typing of Bacillus isolates at the subspecies level. This is particularly true for B. amyloliquefaciens, B. pumilus and B. megaterium but less so for B. anthracis (Table 1).

Table 1

Polymorphisms of the 16S rRNA and *gyrA* gene
Species	Gene	Gene length (bp)	Number of variable sites	Total number of genomes	Genomes number with variable bases		Range of genomes variable sites (%)
Bacillus amyloliquefaciens	16S	480	12	116	4	1–7	0.21–1.46
Bacillus amyloliquefaciens	gyrA	500	59	116	109	2–32	0.4–6.4
Bacillus pumilus	16S	480	21	140	11	1–17	0.21–3.54
Bacillus pumilus	gyrA	500	130	140	87	1–79	0.2–15.8
Bacillus megaterium	16S	480	3	117	48	1–2	0.21–0.42
Bacillus megaterium	gyrA	500	58	117	99	1–14	0.2–2.8
Bacillus anthracis	16S	480	13	226	7	1–5	0.21–1.04
Bacillus anthracis	gyrA	500	65	226	18	1–37	0.2–7.4

The resolution power of Bacillus mock community gyrA amplicon sequencing

Our results above suggest that the primer pair gyrA3 could be used as a molecular marker for diversity analysis of Bacillus. Next, we aimed to design a mock community to test the efficacy of the gyrA3 primer and used the general 16S rRNA primers (V3-V4) as a positive control. For better comparison, we designed two additional primer pairs, gyrA4 and gyrA5, which are very close to the position of the gyrA3 in the gyrA gene (Fig. S8A). We selected eight strains belonging to four species (B. altitudinis, B. licheniformis, B. velezensis and Lysinbacillus pakistanensis) and successfully amplified their gyrA gene by gyrA3, gyrA4 and gyrA5 primer pairs in a routine PCR for selected strains (Fig. 5A). Next, we artificially mixed eight indicated strains in equal proportions to obtain a mock community and then retested the resolution power of the three gyrA primer pairs and the 16S rRNA-specific primers. Our goal was to test whether the primer pair is suitable for determining the diversity of mock community (Fig. 5). Sequencing of the 16S rRNA amplicons showed that 16S rRNA primers could distinguish only five units, as strains LY1 and LY18, LY37 and LY43, and LY39 and LY48 had identical V3-V4 nucleotide sequence (Fig. 5B). The gyrA primer pairs were capable of resolving six units, but gyrA4 and gyrA5 produced amplicons of variable abundance and preferentially amplified LY35 and LY2, respectively(Fig. 5C-E). In comparison, primers gyrA3 also amplified six fragments but the relative abundance of these amplicons was more uniform (Table S6) with the exception of LY2 strain (Fig. 5C). Our data suggest that gyrA3 has potential for Illumina amplicon sequencing of more complex Bacillus communities.

It is believed that the diversity of microorganisms in nature is immense, and therefore its detection still remains a challenge[31]. Molecular tools (e.g., for specific amplification of marker genes) combined with high-throughput sequencing are expected to open the door to the vast diversity of microorganisms[2]. Here, we systematically investigated the potential of the gyrA gene as a marker gene for taxonomic typing of Bacillus isolates and assessment of Bacillus community composition by amplicon sequencing. The novel gyrA3 primer pair is capable of detecting 76 Bacillus species by virtual and colony PCR; hence, our results suggest that the gyrA gene is a good phylogenetic marker for detecting intra- and interspecific diversity of the genus Bacillus.

Sequencing of 16S rRNA amplicons on the Illumina platform has been the mainstream method of studying microbial communities, but this gene also shows shortcomings, especially when we aim to determine diversity within a genus or species[32–37]. This is particularly true for Bacillus species, which exhibit very low interspecific variability[24]. Although 16S rRNA is widely used as a molecular marker for bacterial community analyses[38], its amplicon sequencing can only describe community diversity at the genus level[39]. In contrast, faster evolution and consequently higher diversity of the gyrA gene[1] suggests that this gene might provide higher phylogenetic resolution than the 16S rRNA gene within the genus Bacillus.

In our study, the differential effect of the gyrA gene on Bacillus interspecies and several Bacillus species was better than that 16S rRNA (Fig. 1). Moreover, the results of comparative sequence analysis (including the 16S rRNA and the gyrA gene regions) of the three species: B. amyloliquefaciens, B. pumilus, and B. megaterium, support this prediction, as we find wider range of SNPs in the gyrA genes than in the 16S rRNA (Fig. 4, Fig. S5 and Table 1). Our results are consistent with findings that housekeeping genes, including gyrA, evolve much faster than 16S rRNA genes and are suitable for the identification and typing of closely related species[23] and for the intraspecific resolution of isolates, as previously shown for B. subtilis[28, 29].

To date, there have been only a few reports in which conserved genes (e.g. rpoB and gyrB) have been used as templates for amplicon sequencing of microbial communities[22–24]. These reports show that sequencing of conserved protein coding genes provides a more accurate description of bacterial community composition than 16S rRNA sequencing. Specifically, the rpoB gene has been used in addition to the 16S rRNA molecular marker for high-throughput sequencing studies of species diversity in the Proteobacteria phylum[24].

Although it remains to be investigated whether gyrA is suitable as a molecular marker for the majority of Bacillus species, we identified here its resolution power for 76 Bacillus species (Fig. 3) and propose that the diversity of the Bacillus community can be more accurately assessed by combining 16S rRNA and gyrA amplicon sequencing.

As a molecular marker, the housekeeping gene gyrA also has some limitations. Due to the high variability of gyrA gene sequences, it was difficult to design universal primers that can amplify the entire Bacillus genus, as has been previously described for the Escherichia genus[40]. Although gyrA3 significantly improved the amplification efficiency of Bacillus species compared to previous primers, not all Bacillus species could be detected (Fig. S2 and Fig. 3). On the other hand, the gyrA amplicon region we selected was well suited for resolving subspecies of B. amyloliquefaciens, B. pumilus, and B. megaterium. For B. anthracis, which is known for its low diversity[41] intraspecific resolution was limited (Fig. 4 and Fig. S5). Using amplicon sequencing, gyrA3 also had its own blind area in the detection, for example, it was not good at distinguishing between LY1 and LY18, LY39 and LY48 (Fig. 5). This could be due to the extremely high similarity of the gyrA genes in selected genomes, some of which have even identical gyrA sequences. However, the gyrA3 distinguished very well two strains with highly similar gyrA genes such as LY37 and LY43, which 16S rRNA did not (Fig. 5). Although primers for amplicon sequencing of the gyrA gene had certain limitations, such as not amplifying all selected targets or not reflecting the abundance of added DNA, this study has put forward the advantages that gyrA3 covers the broadest diversity of Bacillus species reported to date.

This study investigated the application of the gyrA gene as a molecular marker in Bacillus subspecies typing and high-throughput sequencing of the Bacillus mock community. The greater ability of gyrA-based analyses to distinguish Bacillus strains at the subspecies level should increase resolution and provide more reliable results for the ecological studies of the genus Bacillus. We believe that the primer pair will have broad applications in the Bacillus research.

Strains and culture condition

The 127 strains used in the study were strains of Bacillus and related genera of Bacillus, most of which were isolated from soil. Some strains labelled with ACCC are deposited in the Agricultural Culture Collection of China (http://www.accc.org.cn) (Table S4).

All strains were grown at 30 ℃ in low-salt Luria-Bertani medium (LB), containing 10 g tryptone, 5 g yeast extract, and 3 g NaCl per litre.

Nucleotide diversity (Pi) analysis

For interspecies analysis, the whole sequences of the 16S rRNA and the gyrA genes of 20 Bacillus species (361 genomes) were downloaded from the NCBI genome database (Table S1). For intraspecies analysis, three species were selected and whole sequences of 16S rRNA and gyrA gene were downloaded from the database: Bacillus amyloliquefaciens (88 genomes), Bacillus licheniformis (71 genomes) and Bacillus pumilus (97 genomes) (Table S1). Alignment of these sequences was conducted using online alignment tool Kalign on EMBL (https://www.ebi.ac.uk/Tools/msa/). Subsequently, nucleotide diversity (Pi) was estimated with DnaSP 6 (v.6) using a window size of 100 bp and a step size of 10 bp[42].

DNA extraction of strains, PCR and gel electrophoresis

Genomic DNA was extracted using the Omega Bacterial DNA Kit D3350 (Omega, Bio-tek, Norcross, GA, USA), and the concentration and quality of DNA were determined using a NanoDrop 2000 spectrophotometer (Wilmington, DE, USA).

The reaction mixture for PCR amplification was prepared in 25 µL containing 1 µL of DNA, 2 µL of each primer (the primer pairs gyrA1 were gyrA1-F: 5′-CAGTCAGGAAATGCGTACGTCCTT-3’[29] and gyrA1-R: 5′- GTATCCGTTGTGCGTCAGAGTAAC-3’[28], the primer pairs gyrA2 were gyrA2-F: 5′-CAGTCAGGAAATGCGTACGTCCTT-3’[29] and gyrA2-R: 5′-CAAGGTAATGCTCCAGGCATTGCT-3’[29], the primer pairs gyrA3 were gyrA3-F: 5′-GCDGCHGCNATGCGTTAYAC-3’ and gyrA3-R: 5′-ACAAGMTCWGCKATTTTTTC-3’, the primers for the 16S rRNA gene were 27F: 5′-AGAGTTTGATCCTGGCTCAG-3’ and 1492R: 5′-GGTTACCTTGTTACGACTT-3’), 12.5 µL Green Taq Mix (http://www.vazyme.com), and 7.5 µL deionized water.

PCR was performed under the following conditions: Predenaturation at 94 ℃ for 5 min; denaturation at 94 ℃ for 30 s; annealing at 50 ℃ for 30 s; elongation at 72 ℃ for 40 s (35 cycles); and elongation at 72 ℃ for 7 min.

In Silico PCR

The gyrA sequence database was constructed as a virtual community, containing 5062 full-length sequences from the National Center for Biotechnology Information (NCBI) belonging to 226 Bacillus species (Table S2) (the method for obtaining the database sequences is described in detail in the SNP analysis section below). First, the degenerate primer pair gyrA3 was converted to primers that do not contain degenerate bases, and the converted gyrA3 is listed in Table S7. Subsequently, the primer pairs gyrA1, gyrA2 and the converted gyrA3 were aligned with the gyrA database using NCBI-blast + software (v.2.9.0). The match at 18 bases of both, the forward and reverse primer, was considered amplifiable by the primer pair.

Phylogenetic analysis

In this study, phylogenetic analysis of genes was performed using MEGA (v.5.05)[43] for Neighbor-Joining and the reliability of clades was tested using 1000 bootstrap replicates. Furthermore, annotation and beautification of trees was performed using programs available at the iTol online site (https://itol.embl.de)[43].

SNP analysis

A total of 734 available genomes of 4 different Bacillus species were downloaded from the NCBI database using the NCBI-genome-download script (https://github.com/kblin/ncbigenome-download/) (Table S5) including 124 genomes of B. amyloliquefaciens, 167 genomes of B. pumilus, 150 genomes of B. megaterium and 293 genomes of B. anthracis. Then, the location information of 16S rRNA and gyrA gene sequences on the genomes were determined by aligning 16S rRNA and gyrA sequences of the same species using NCBI-blast+ (v.2.9.0) and the 16S rRNA / gyrA gene sequences were extracted using the Fasta Extract tool in TBtools (v1.0986853)[44]. After filtration, genomes that did not carry complete 16S rRNA and gyrA genes, were removed and 600 genomes were retained, including 116 genomes of B. amyloliquefaciens, 140 genomes of B. pumilus, 117 genomes of B. megaterium and 226 genomes of B. anthracis (Table S5).

The alignment of the 16S rRNA (base sites:330–810) and gyrA genes (base sites:350–850) within each species was performed using the L-INS-I method of MAFFT (v7.487) (https://mafft.cbrc.jp/alignment/software/). The two gene sequences (16S rRNA, gyrA) on the same genome were selected as representative sequences (the top sequence of each figure was the representative sequence), and the base mismatch sites on other sequences were marked with color after comparison with the representative sequence. The genomes of the representative sequences in different species were GCF_017815555.1 (B. amyloliquefaciens), GCF_014269965.1 (B. pumilus), GCF_002577645.1 (B. megaterium) and GCF_000007845.1 (B. anthracis) (Fig. 4 and Fig. S5).

Construction of the DNA-based Bacillus mock community, amplicon sequencing and data analysis

For the Bacillus mock community, we selected 8 strains with known genome sequences. Genomic DNA from 8 strains was extracted and its quality and quantity determined. Eight genomic DNAs were pooled in equal amounts after being diluted to approximately the same concentrations.

The hypervariable region V3-V4 of the 16S rRNA gene was amplified with the universal primers 338F: 5′-CCTACGGRRBGCASCAGKVRVGAAT-3’ and 806R: 5′-GGACTACNVGGGTWTCTAATCC-3′. The gyrA gene from 8 strains of the mock community was amplified with primer pairs gyrA3 (see above), gyrA4 (F:5′-TAYGCRATGAGYRTHATYGT-3’ and R: 5′- TTBGTNGCCATHCCDACMGC-3’), and gyrA5 (F:5′-GCDGCNGCVATGCGTTAYAC-3’ and R: 5′- CGNAGRTYBGTAATDCCDTC-3’). Sequencing was performed on an Illumina Miseq PE300 instrument.

Raw data were processed using the UPARSE pipeline (http://drive5.com/usearch/manual/uparse_pipeline.html)[37] to obtain the ZOTUs table. This included truncating and assembling, removing double end primers, filtering low-quality sequences, finding unique sequences, denoising, and creating the ZOTUs with 100% similarity, which were then used to create the ZOTUs table. Denoising (error-correction) of the amplicon reads was performed by using the Unoise3 algorithm[36] to identify and correct reads with sequencing errors and remove chimeras.

The eight strains of the mock community with known complete 16S rRNA and gyrA gene sequences were used to annotate the ZOTUs (Dataset S2 and S3). The annotation of ZOTUs sequences was performed with NCBI-blast+ (v.2.9.0). The ZOTUs of the same species were pooled into a single unit after annotation. Finally, we used box plots to show the community structure and the characteristics of the different primer pairs during sequencing. The box plots were drawn using R (v.4.0.3).

SNP: single nucleotide polymorphisms; NCBI: National Center for Biotechnology Information.

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and materials

Data and materials used in the analysis are available upon request from the corresponding authors for the purposes of reproducing or extending the analysis. Amplicon sequencing reads from the 16S rRNA gene and gyrA gene are available at NCBI Sequence Read Archive under accession number PRJNA823863.

Competing Interests

The authors declare that they have no conflicts of interest.

Funding

This work was financially supported by the National Nature Science Foundation of China (31972512 and 42090064), the Fundamental Research Funds for the Central Universities (KYXK202009). The Innovative Research Team Development Plan of the Ministry of Education of China (Grant No. IRT_17R56). PS and IMM were supported by the Program Grant P4-0116 funded by the Slovenian national research agency (ARRS).

Authors’ contributions

YL and ZX designed the study; YL and YX performed the experiments. YL, YM, WX and PS analysed the data and created the figures. YL, PS and ZX wrote the first draft of the paper; PS, ZX, RZ, QS and IM revised the paper.

Acknowledgements

Not applicable

Timmis K, Ramos JL. The soil crisis: the need to treat as a global health problem and the pivotal role of microbes in prophylaxis and therapy. Microb Biotechnol. 2021;14:769–97. https://doi.org/10.1111/1751-7915.13771
Klindworth A, Pruesse E, Schweer T, Peplies J, Quast C, Horn M, et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 2013;41:e1. https://doi.org/10.1093/nar/gks808
Peer Y Van De, Chapelle S, Wachter R De. A quantitative map of nucleotide substitution rates in bacterial rRNA. Nucleic Acids Res. 1996;24:3381–91. https://doi.org/10.1093/nar/24.17.3381
McDonald D, Price MN, Goodrich J, Nawrocki EP, Desantis TZ, Probst A, et al. An improved greengenes taxonomy with explicit ranks for ecological and evolutionary analyses of bacteria and archaea. ISME J. 2012;6:610–8. https://doi.org/10.1038/ismej.2011.139
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 2013;41:590–6. https://doi.org/10.1093/nar/gks1219
Yoon SH, Ha SM, Kwon S, Lim J, Kim Y, Seo H, et al. Introducing EzBioCloud: a taxonomically united database of 16S rRNA gene sequences and whole-genome assemblies. Int J Syst Evol Microbiol. 2017;67:1613–7. https://doi.org/10.1099/ijsem.0.001755
Cole JR, Wang Q, Fish JA, Chai B, McGarrell DM, Sun Y, et al. Ribosomal Database Project: Data and tools for high throughput rRNA analysis. Nucleic Acids Res. 2014;42:633–42. https://doi.org/10.1093/nar/gkt1244
Pinto AJ, Raskin L. PCR biases distort bacterial and Archaeal community structure in pyrosequencing datasets. PLoS One. 2012;7:e43093. https://doi.org/10.1371/journal.pone.0043093
Sun D-L, Jiang X, Wu QL, Zhou N-Y. Intragenomic heterogeneity of 16S rRNA genes causes overestimation of prokaryotic diversity. Appl Environ Microbiol. 2013;79:5962–9. https://doi.org/10.1128/AEM.01282-13
Johnson JS, Spakowicz DJ, Hong B, Petersen LM, Demkowicz P, Chen L, et al. Evaluation of 16S rRNA gene sequencing for species and strain-level microbiome analysis. Nat Commun. 2019;10:5029. http://dx.doi.org/10.1038/s41467-019-13036-1
Eren AM, Maignien L, Sul WJ, Murphy LG, Grim SL, Morrison HG, et al. Oligotyping: differentiating between closely related microbial taxa using 16S rRNA gene data. Methods Ecol Evol. 2013;4:1111–9. https://doi.org/10.1111/2041-210X.12114
Mandic-Mulec I, Stefanic P, van Elsas JD. Ecology of bacillaceae. Microbiol Spectr. 2015;59–85. https://doi.org/10.1128/microbiolspec.TBS-0017-2013
Chun J, Bae KS. Phylogenetic analysis of Bacillus subtilis and related taxa based on partial gyrA gene sequences. Antonie Van Leeuwenhoek. 2000;78:123–7. https://doi.org/10.1023/A:1026555830014
Kasai H, Tamura T, Harayama S. Intrageneric relationships among Micromonospora species deduced from gyrB-based phylogeny and DNA relatedness. Int J Syst Evol Microbiol. 2000;50:127–34. https://doi.org/10.1099/00207713-50-1-127
Yamamoto S, Kasai H, Arnold DL, Jackson RW, Vivian A, Harayama S. Phylogeny of the genus Pseudomonas: intrageneric structure reconstructed from the nucleotide sequences of gyrB and rpoD genes. Microbiology. 2000;146:2385–94. https://doi.org/10.1099/00221287-146-10-2385
Hurtle W, Bode E, Kulesh DA, Kaplan RS, Garrison J, Bridge D, et al. Detection of the Bacillus anthracis gyrA gene by using a minor groove binder probe. J Clin Microbiol. 2004;42:179–85. https://doi.org/10.1128/JCM.42.1.179-185.2004
Stefanic P, Mandic-Mulec I. Social interactions and distribution of bacillus subtilis pherotypes at microscale. J Bacteriol. 2009;191:1756–64. https://doi.org/10.1128/JB.01290-08
Stefanic P, Decorosi F, Viti C, Petito J, Cohan FM, Mandic-Mulec I. The quorum sensing diversity within and between ecotypes of Bacillus subtilis. Environ Microbiol. 2012;14:1378–89. https://doi.org/10.1111/j.1462-2920.2012.02717.x
Levican A, Collado L, Figueras MJ. Arcobacter cloacae sp. nov. and Arcobacter suis sp. nov., two new species isolated from food and sewage. Syst Appl Microbiol. 2013;36:22–7. http://dx.doi.org/10.1016/j.syapm.2012.11.003
Stefanic P, Kraigher B, Lyons NA, Kolter R, Mandic-Mulec I. Kin discrimination between sympatric Bacillus subtilis isolates. Proc Natl Acad Sci U S A. 2015;112:14042–7. https://doi.org/10.1073/pnas.1512671112
Ménard A, Buissonnière A, Prouzet-Mauléon V, Sifré E, Mégraud F. The GyrA encoded gene: a pertinent marker for the phylogenetic revision of Helicobacter genus. Syst Appl Microbiol. 2016;39:77–87. http://dx.doi.org/10.1016/j.syapm.2015.09.008
Ogier JC, Pagès S, Galan M, Barret M, Gaudriault S. rpoB, a promising marker for analyzing the diversity of bacterial communities by amplicon sequencing. BMC Microbiol. 2019;19:171. https://doi.org/10.1186/s12866-019-1546-z
Poirier S, Rué O, Peguilhan R, Coeuret G, Zagorec M, Champomier-Vergès MC, et al. Deciphering intra-species bacterial diversity of meat and seafood spoilage microbiota using gyrB amplicon sequencing: a comparative analysis with 16S rDNA V3-V4 amplicon sequencing. PLoS One. 2018;13:e0204629. https://doi.org/10.1371/journal.pone.0204629
Vos M, Quince C, Pijl AS, de Hollander M, Kowalchuk GA. A comparison of rpoB and 16S rRNA as markers in pyrosequencing studies of bacterial diversity. PLoS One. 2012;7:e30600. https://doi.org/10.1371/journal.pone.0030600
Cozzarelli NR. DNA gyrase and the supercoiling of DNA. Science. 1980;207: 953–60. https://doi.org/10.1126/science.6243420
Oslizlo A, Stefanic P, Vatovec S, Beigot Glaser S, Rupnik M, Mandic-Mulec I. Exploring ComQXPA quorum-sensing diversity and biocontrol potential of Bacillus spp. isolates from tomato rhizoplane. Microb Biotechnol. 2015;8:527–40. https://doi.org/10.1111/1751-7915.12258
Borshchevskaya LN, Kalinina AN, Sineokii SP. Design of a PCR test based on the gyrA gene sequence for the identification of closely related species of the Bacillus subtilis group. Appl Biochem Microbiol. 2013;49:646–55. https://doi.org/10.1134/S0003683813070028
Ansaldi M, Marolt D, Stebe T, Mandic-Mulec I, Dubnau D. Specific activation of the Bacillus quorum-sensing systems by isoprenylated pheromone variants. Mol Microbiol. 2002;44:1561–73. https://doi.org/10.1046/j.1365-2958.2002.02977.x
Roberts MS, Nakamura LK, Cohan FM. Bacillus mojavensis sp. nov., distinguishable from Bacillus subtilis by sexual isolation, divergence in DNA sequence, and differences in fatty acid composition. Int J Syst Bacteriol. 1994;44:256–64. https://doi.org/10.1099/00207713-44-2-256
De Clerck E, Vanhoutte T, Hebb T, Geerinck J, Devos J, De Vos P. Isolation, characterization, and identification of bacterial contaminants in semifinal gelatin extracts. Appl Environ Microbiol. 2004;70:3664–72. https://doi.org/10.1128/AEM.70.6.3664-3672.2004
Widder S, Allen RJ, Pfeiffer T, Curtis TP, Wiuf C, Sloan WT, et al. Challenges in microbial ecology: Building predictive understanding of community function and dynamics. ISME J. 2016;10:2557–68. https://doi.org/10.1038/ismej.2016.45
Callahan BJ, McMurdie PJ, Holmes SP. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 2017;11:2639–43. http://dx.doi.org/10.1038/ismej.2017.119
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3. https://doi.org/10.1038/nmeth.3869
Amir A, Daniel M, Navas-Molina J, Kopylova E, Morton J, Xu ZZ, et al. Deblur rapidly resolves single-nucleotide community sequence patterns. mSystems. 2017;2:e00191-16. https://doi.org/10.1128/mSystems.00191-16
Edgar R. UCHIME2: improved chimera prediction for amplicon sequencing. bioRxiv. 2016;074252. https://doi.org/10.1101/074252
Edgar R. UNOISE2: improved error-correction for Illumina 16S and ITS amplicon sequencing. bioRxiv. 2016;081257. https://doi.org/10.1101/081257
Edgar RC. UPARSE: Highly accurate OTU sequences from microbial amplicon reads. Nat Methods. 2013;10:996–8. https://doi.org/10.1038/nmeth.2604
Clarridge JE. Impact of 16S rRNA gene sequence analysis for identification of bacteria on clinical microbiology and infectious diseases. Clin Microbiol Rev. 2004;17:840–62. https://doi.org/10.1128/CMR.17.4.840-862.2004
Gupta S, Mortensen MS, Schjørring S, Trivedi U, Vestergaard G, Stokholm J, et al. Amplicon sequencing provides more accurate microbiome information in healthy children compared to culturing. Commun Biol. 2019;2:291. http://dx.doi.org/10.1038/s42003-019-0540-1
Johnning A, Kristiansson E, Fick J, Weijdegård B, Larsson DGJ. Resistance mutations in gyrA and parC are common in Escherichia communities of both fluoroquinolone-polluted and uncontaminated aquatic environments. Front Microbiol. 2015;6:1355. https://doi.org/10.3389/fmicb.2015.01355
Lista F, Faggioni G, Valjevac S, Ciammaruconi A, Vaissaire J, Le Doujet C, et al. Genotyping of Bacillus anthracis strains based on automated capillary 25-loci Multiple Locus Variable-Number Tandem Repeats Analysis. BMC Microbiol. 2006;6:1–14. https://doi.org/10.1186/1471-2180-6-33
Rozas J, Ferrer-Mata A, Sanchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, et al. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol. 2017;34:3299–302. https://doi.org/10.1093/molbev/msx248
Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011;28:2731–9. https://doi.org/10.1093/molbev/msr121
Chen C, Chen H, Zhang Y, Thomas HR, Frank MH, He Y, et al. TBtools: An Integrative Toolkit Developed for Interactive Analyses of Big Biological Data. Mol Plant. 2020;13:1194–202. https://doi.org/10.1016/j.molp.2020.06.009

No competing interests reported.

Additionalfile1.zip
Additional file1:FigureS1.The agarose gel electrophoresis of DNA amplification products of three gyrA gene primer pairs: gyrA1 (A), gyrA2 (B) and gyrA3 (C). 1-7 indicates strains LY18, LY25, LY37, FZB42, SXL408, SXL277 and ACCC01043. FigureS2.The amplification ability of the primer pair gyrA1 and gyrA2 in Bacillus species by computer simulation. The eight Bacillus species were amplified by gyrA1 and nine Bacillus species were amplified by gyrA2. The phylogenetic trees were constructed based on gyrA gene (2403bp). FigureS3.Polymorphisms in the 16S rRNA and gyrA gene’s regions of B. amyloliquefaciens (116 genomes). (A) The proportion of variation at different base sites along 16S rRNA gene V3-V4 region. (B) The proportion of variation at the different base positions of the gyrA gene region. (C) The proportion of variants along 16S rRNA (red column) and gyrA gene (blue column) region in different genomes. FigureS4.Polymorphisms in the 16S rRNA and gyrA gene’s region of B. pumilus (140 genomes). (A) The proportion of variation at different base sites along the 16S rRNA gene (B) and the gyrA gene. (C) The proportion of variable base sites in 16S rRNA (red column) and gyrA (blue column) sequences in different genomes. FigureS5.Alignment of the 16S rRNA and gyrA sequences of B. megaterium (117 genomes) and B. anthracis (226 genomes). Sequences were aligned using L-INS-I method of MAFFT (v7.487). Analysis involved positions along 16S rRNA from 330-810 and along gyrA from 350-850. On the left side of the graph GCF reference numbers of genomes were displayed, with reference sequence GCF_002577645.1 displayed at the top for B. megaterium and GCF_000007845.1 for B. anthracis. The display of base site variation was drawn using MEGA (v.5.05). The colored line (red or blue) indicates a variation of specific base sites as compared to the reference sequence. FigureS6.Polymorphisms in the 16S rRNA and gyrA gene’s regions of B. megaterium (117 genomes). (A) The proportion of variation at different base sites along 16S rRNA gene V3-V4 region. (B) The proportion of variation at the different base positions of the gyrA gene region. (C) The proportion of variants along 16S rRNA (red column) and gyrA gene (blue column) region in different genomes. FigureS7.Polymorphisms in the 16S rRNA and gyrA gene’s region of B. anthracis (226 genomes). (A) The proportion of variation at different base sites along the 16S rRNA gene (B) and the gyrA gene. (C) The proportion of variable base sites in 16S rRNA (red column) and gyrA (blue column) sequences in different genomes. FigureS8.Description of gyrA gene primer pairs for amplicon sequencing. (A) The position of gyrA amplicons obtained by primer pairs gyrA3, gyrA4 and gyrA5. (B) PCR amplification of eight Bacillus strains by gyrA gene primer pairs. The orange square indicates positive and white square unsuccessful amplification. The phylogenetic tree was contracted based on the complete gyrA genes (2469bp) (Dataset S3) by using Neighbor-Joining method in MEGA 5.05 software. The reliability of clades was tested by the 1000 bootstrap replications.
Additionalfile2.zip
Additional file2:DatasetS1.The 16S rRNA gene sequences of the seven Bacillus strains. DatasetS2.The 16S rRNA gene sequences of the eight Bacillus strains for mock Bacillus community. DatasetS3.The gyrA gene sequences of the eight Bacillus strains for mock Bacillus community.
Additionalfile3.zip
Additional file3:TableS1.The genomes information for nucleotide variability analysis of gyrA gene and 16S rRNA gene in inter- and intraspecies Bacillus. TableS2.The gyrA gene database information of Bacillus for computer simulation. TableS3.The information of sequences amplified by three gyrA gene primer pairs in computer simulation. TableS4.The information of 127 strains for colony PCR. TableS5.The genomes information for SNP analysis of gyrA gene and 16S rRNA gene in four Bacillus species. TableS6.The analysis results of mock community amplicon sequencing. TableS7.The general primers converted from degenerate primer pair gyrA3.

PDF

Version 1

posted 20 Apr, 2022

You are reading this latest preprint version

Housekeeping gene gyrA, a potential molecular marker for Bacillus ecology study

Status:

Version 1

Abstract

Figures

Background

Results

The gyrA gene of the Bacillus genus shows higher variation rates than 16S rRNA

First comparative tests of three primer pairs for the detection of Bacillus species

Specificity range of the gyrA3 primer pair by using PCR and in silico PCR

The gyrA gene provides better intraspecific phylogenetic resolution than the 16S rRNA gene among certain species

The resolution power of Bacillus mock community gyrA amplicon sequencing

Discussion

Conclusions

Materials And Methods

Strains and culture condition

Nucleotide diversity (Pi) analysis

DNA extraction of strains, PCR and gel electrophoresis

In Silico PCR

Phylogenetic analysis

SNP analysis

Construction of the DNA-based Bacillus mock community, amplicon sequencing and data analysis

Abbreviations

Declarations

References

Competing Interests

Supplementary Files

Status:

Version 1