Introduction {#s1}
============

Bacterial pathogens continue to cause problems for humans by the continuous evolution of known pathogens and emergence of new pathogens (Jackson et al., [@B9]). Approximately 160 new bacterial infectious diseases have been discovered between 1940 and 2010, and emerging infectious diseases constitute a significant burden on global economies and public health (Jones et al., [@B10]). While the mechanisms by which many pathogenic bacteria cause disease are relatively well-described, the evolutionary origin of virulence is not fully understood (Wilson et al., [@B34]; Persson et al., [@B20]). Bacterial virulence determinants causing disease in humans may evolve through complex host-pathogen interactions (Wilson and Salyers, [@B33]; Mikonranta et al., [@B17]) but their origin may also lie in non-human reservoirs amongst the environmental microbiota (Martínez, [@B15]).

The large majority of studies on virulence in the environment have focused on the detection and survival of clinically relevant bacteria in non-host environments (Garcia-Aljaro et al., [@B5]; Vezzulli et al., [@B30]; Davidson et al., [@B4]). Only few studies (e.g., Persson et al., [@B20]; Søborg et al., [@B23]) have targeted the environmental gene pool of pathogenic traits, hence evidence of the presence of virulence genes in environmental microbiomes is sparse.

By PCR and cultivation-based approaches we have previously investigated the presence of virulence factors in soil and freshwater environments (Søborg et al., [@B23], [@B24]). However, to better understand (i) the significance of virulence gene homologs for bacterial survival in non-human environments and (ii) the role these environments play for the evolution of bacterial human pathogens, a much broader investigation across multiple environmental microbiomes is required. The increasing availability of environmental metagenomes provides an ideal opportunity for interrogating the occurrence of bacterial human virulence genes in natural environments. Hence, we performed an *in silico* analysis targeting 24 bacterial human virulence genes in 46 environmental bacterial metagenomes. The analysis showed that not only were the bacterial virulence gene homologs widespread in natural (non-human) environments, evidence of environmental transcription was also observed as was a high degree of sequence conservation between some of the environmental and clinical genes. Thus, the results provide documentation of an important role of virulence gene homologs in environmental microbiomes and suggest they may be precursors of virulence genes seen in clinically relevant pathogenic bacteria.

Materials and methods {#s2}
=====================

Blast searches for virulence determinants in metagenomic datasets
-----------------------------------------------------------------

Complete translated sequences of 24 bacterial virulence genes were used as queries in TBLASTN searches in 46 environmental metagenomes of the Community Cyberinfrastructure for Advanced Microbial Ecology Research and Analysis (CAMERA) database (Sun et al., [@B25]; note that the database has been transferred to iMicrobe (<http://data.imicrobe.us>). The genes covered several classes of virulence determinants including toxin genes (*hlyA, hylA, ply, stx1, stx2, vacA, hlgB*), adhesin genes (*papA, papH, papC, fimH, focG*), secretion genes (*invA, spiA, sipB, sipC, ssaD*), inflammatory genes (*wbdI, rfbE*), regulatory genes (*phoP, phoQ*) and resistance genes (*pqaB, mecA, yfbI*) originally described from *Salmonella* spp., *Escherichia coli/Shigella, Helicobacter* spp., *Streptococcus pneumoniae*, and *Staphylococcus aureus* (Table [1](#T1){ref-type="table"}).

###### 

**Bacterial human virulence genes used for TBLASTN analysis of environmental metagenomic datasets**.

  **Virulence factor**                                                            **Gene**                 **Protein size (aa)**   **Reference**
  ------------------------------------------------------------------------------- ------------------------ ----------------------- ----------------
  Shiga-like toxin 1                                                              *stx1*                   315                     ACU32680.1
  Shiga-like toxin 2                                                              *stx2*                   319                     YP_003234845.1
  Vacuolating cytotoxin                                                           *vacA*                   1296                    YP_002301516.1
  Haemolysin A                                                                    *hlyA*                   1024                    NP_755445.1
  Hyaluronidase                                                                   *hylA*                   1066                    NP_344851.1
  Pneumolysin                                                                     *ply*                    471                     NP_346351.1
  Gamma-haemolysin component C precursor                                          *hlgB*                   325                     YP_041861.1
  P fimbriae                                                                      *papC*                   839                     NP_755465.1
  P fimbriae                                                                      *papA*                   199                     NP_755467.1
  P fimbriae                                                                      *papH*                   195                     NP_755466.1
  Type 1 fimbriae                                                                 *fimH*                   303                     NP_757248.1
  F1C fimbriae                                                                    *focG*                   167                     NP_753158.1
  Type III, secretion apparatus                                                   *invA*                   685                     YP_151923.1
  Type III, secretion apparatus (SPI[^a^](#TN1){ref-type="table-fn"}-2)           *spiA*                   497                     YP_150710.1
  Type III, secretion apparatus (SPI[^a^](#TN1){ref-type="table-fn"}-2)           *ssaD*                   403                     YP_150709.1
  Type III, translocators and effectors (SPI[^a^](#TN1){ref-type="table-fn"}-1)   *sipB*                   593                     YP_151912.1
  Type III, translocators and effectors (SPI[^a^](#TN1){ref-type="table-fn"}-1)   *sipC*                   409                     YP_151911.1
  Regulator of virulence determinants                                             *phoP*                   224                     YP_150858.1
  Regulator of virulence determinants                                             *phoQ*                   487                     YP_150859.1
  Part of O-antigen 111                                                           *wbdI*                   149                     AAC44881.1
  Part of O-antigen 157                                                           *rbfE*                   364                     NP_310868.1
  L-Ara4N transferase                                                             *pqaB* (*arnT* region)   548                     YP_149876.1
  L-Ara4N transferase                                                             *yfbI* (*arnT* region)   550                     NP_754685.1
  Penicillin-binding protein                                                      *mecA*                   668                     YP_039515.1

*Salmonella Pathogenicity Island (SPI)*.

The metagenomes represented terrestrial, marine, limnic as well as several extreme environments (Supplementary Table [1](#SM1){ref-type="supplementary-material"}). Environments potentially impacted by human or livestock discharges were excluded from the analysis as were metagenomes restricted to specific bacterial groups.

Default NCBI blastall parameters (expected threshold: 10; word size: 3; matrix: BLOSUM62; gap open cost: 11; gap extend cost: 1) were applied. Hits with an expect (*E*) value \< 10^−6^ and a sequence similarity \>30% were considered to be significant (Persson et al., [@B20]). We further applied the criterium of a minimal length of 35 amino acids (aa) to weed-out the shortest 454 pyrosequencing sequences. Longest TBLASTN hits were aligned against the query sequences, and query coverage and percent aa identity determined. As a final criterium to confirm that the environmental hits indeed matched known classes of virulence genes, NCBI BLASTP searches (default parameters) were done using environmental hits as queries.

Blast searches for clinically relevant bacteria in metagenomic datasets
-----------------------------------------------------------------------

Bacterial groups from which the virulence genes were originally described (*Salmonella* spp., *E. coli/Shigella, Helicobacter* spp., *Streptococcus pneumoniae* and *Staphylococcus aureus*) were targeted by BLASTN (default parameters, *E*-value \< 10^−6^) in the metagenomes. Partial 16S/23S rRNA gene sequences usually amplified with published PCR primers were used as queries (Supplementary Table [2](#SM2){ref-type="supplementary-material"}). BLASTN hits were checked for presence of the primer sequences as well as for similarity to the 16S or 23S gene sequences of the clinical bacteria.

Phylogenetic analysis of virulence gene homologs
------------------------------------------------

Phylogenetic trees showing the relationships between the translated environmental sequences and the corresponding sequences of the clinical genes, were prepared with MEGA Version 5.01 (Tamura et al., [@B26]) using the Minimum Evolution algorithm (pairwise deletion). In addition, the Maximum Parsimony method of MEGA Version 5.01 was employed to create phylogenetic trees (complete deletion; 500 replicates) to infer the evolutionary history of the homologs in comparison with the clinically relevant genes.

Statistical analysis
--------------------

To look for trends that potentially could explain differences between metagenomes with respect to presence of virulence gene homologs, a Principal Component Analysis (PCA) was carried out. The PCA was performed with the Matlab based chemometric software package Latentix 2.0 (Latent5 Aps, Copenhagen N, Denmark, [www.latentix.com](http://www.latentix.com)) using the data set shown in Table [2](#T2){ref-type="table"}. A PCA model, *X*, was built as the product of an I × K column-wise orthogonal score matrix T and the transpose of the column-wise orthogonal J × K loading matrix P, where I denotes the number of samples (metagenomes), J the number of variables (virulence gene homologs), and K the number of principal components. E is the matrix of model residuals of size I × J (Equation 1). Zeros were inserted in *X* if genes were not detected in a specific environment.

X

=

T

P

T

\+

E

###### 

**Frequency of putative bacterial human virulence genes in environmental metagenomic datasets**.

  **Metagenome**                           **Toxin**                               **Adhesin**   **Secretion**                        **Regulatory**   **Inflammatory**                      **Resistance**                          **Genes detected**                                                                                                                                                                                                                                                                                                                                                 
  ---------------------------------------- --------------------------------------- ------------- ------------------------------------ ---------------- ------------------------------------- --------------------------------------- --------------------------------------- -------------------------------------- ------------------------------------- --------------------------------------- -------------------------------------- -------------------------------------- ------------------------------------- ------------------------------------- ------------------------------------- ----- ----
  Microbialites (Marine microbial mat)                                                                                                                                                                                                                                       20                                                                                                                                                                                                                                                                                                                         1
  FLAS (Saltern)                                                                                                                                                                                                                                                                                                                                                                                  420                                    840                                                                                                                                                            2
  LineIsland (Seawater)                                                                                                                                                                                                                                                                                                                                   4[^\*^](#TN2){ref-type="table-fn"}                                             4[^\*^](#TN2){ref-type="table-fn"}                                                                                                                             2
  PeruMarginSediment (Seawater)                                                                                                                                                                                                                                              16                                                                                                                                                          16                                                                                                                                                             2
  Mountain lake (Freshwater)                                                                                                                                                                                                                                                 37                                                                           74                                                                             258                                                                                                                                                            3
  GeneExpression (Seawater)                                                                                                                                                                                                                                                  81[^\*^](#TN2){ref-type="table-fn"}                                          162[^\*^](#TN2){ref-type="table-fn"}                                           162[^\*^](#TN2){ref-type="table-fn"}                                                                                                                           3
  Contaminated soil                        23                                                                                                                                                                                                                                                                                                             69                                                                             46                                                                                                                                                             3
  AmazonRiverPlume (Estuarine)                                                                                                                                                                                                                                               58                                                                           43                                      1                                      41[^\*^](#TN2){ref-type="table-fn"}                                                                                                                            4
  Sapelo2008 (Seawater)                                                                                                                                                                                                                                                      186                                                                          127                                     17                                     214                                                                                                                                                            4
  SargassoSea (Seawater)                                                                                                                                                                                                                                                     498                                                                          134                                                                            197[^\*^](#TN2){ref-type="table-fn"}   8[^\*^](#TN2){ref-type="table-fn"}                                                                                      4
  DeepMed (Seawater)                       555                                                                                                                                                                                                                                                                                                            416                                                                            1527                                                                                                                                                     139   4
  Yellowstone (Hot spring)                                                                                                                                                                                                                                                                                                                                3554                                    135                                    779                                                                                                                                                      34    4
  AcidMine (Biofilm)                                                                                                                                                                                                                                                         212                                                                          734                                     28                                     295                                                                                                                                                      34    5
  PacificOcean (Seawater)                  4[^\*^](#TN2){ref-type="table-fn"}                                                                                                                                                                                                106                                                                          85                                      1                                      51[^\*^](#TN2){ref-type="table-fn"}                                                                                                                            5
  PBSM (Marine beach sand)                 167                                                                                                                                                                                                                               167                                                                          834                                                                            500                                                                                                                167                                         5
  MontereyBay transect (Seawater)                                                                                                                                                                                                                                            33                                                                           101                                     2                                      512                                    2                                                                                                                 6     6
  HOT:454 (Seawater)                                                                                                                                                                                                                                                         37                                                                           177                                     4                                      956                                    5                                                                                                                 10    6
  N. Pac. Line67 (Seawater)                1                                                                                                                                                                                                                                 15                                                                           40                                                                             29                                     2                                                                           1                                           6
  BotanyBay:Sanger (Seawater)                                                                                                                                                                                                                                                13                                     3                                     764                                     81                                     527[^\*^](#TN2){ref-type="table-fn"}                                                                                                                     61    6
  AlvinellaPompejana (Hydrothermal vent)                                                                                                                                                                                                                                     10                                     4                                     286                                                                            1078                                                                                                               28                                    65    6
  HydrothermalVent                         101                                                                                                                                                                                                                               182                                                                          706                                     40                                     948                                                                                                                                                      20    6
  TermiteGut                                                                                                                                           16                                                                                                                    246                                                                          410                                     66                                     689                                                                                                                                                      98    6
  SalternMetagenome                                                                                                                                    5[^\*^](#TN2){ref-type="table-fn"}                                                                                    7[^\*^](#TN2){ref-type="table-fn"}                                           35[^\*^](#TN2){ref-type="table-fn"}     14[^\*^](#TN2){ref-type="table-fn"}                                           2[^\*^](#TN2){ref-type="table-fn"}    5[^\*^](#TN2){ref-type="table-fn"}    14[^\*^](#TN2){ref-type="table-fn"}         7
  MontereyBay (Seawater)                                                                                                                                                                                                                                                     123                                                                          337                                     7                                      88                                     4                                                                           4                                     21    7
  GuaymasBasin (Hydrothermal Vent)         1                                                                                                                                                                                                                                 49                                                                           122                                     27                                     729                                    8                                                                                                                 34    7
  Bacterioplankton (Seawater)              3                                                                                                                                                                                                                                 110                                    3                                     250                                     3                                      503                                    3                                                                                                                 7     8
  BATS (Seawater)                          1                                                                                                           1                                                                                                                     51                                                                           151                                     3                                      445                                    2                                                                                                                 6     8
  BermudaOceanic (Seawater)                4                                                                                                                                                                                                                                 92                                                                           454                                     33                                     1141                                                                         2                                     4                                     122   8
  IceMetagenome (Freshwater ice)                                                                                                                       4                                                                             4                                       138                                                                          960                                     13                                     605                                                                          13                                    21                                          8
  HypersalineMat (Seawater)                83[^\*^](#TN2){ref-type="table-fn"}                                                                                                                                                                                               142[^\*^](#TN2){ref-type="table-fn"}                                         1460[^\*^](#TN2){ref-type="table-fn"}   107[^\*^](#TN2){ref-type="table-fn"}   973[^\*^](#TN2){ref-type="table-fn"}   12[^\*^](#TN2){ref-type="table-fn"}                                         12[^\*^](#TN2){ref-type="table-fn"}   71    8
  CellCapture (Marine sediment)                                                                                                                        63                                                                                                                    333                                    10                                    1553                                    135                                    1053                                                                         156                                   177                                   42    9
  Drifting-ESP (Seawater)                  0.45                                                                                                        0.45                                                                          0.23                                    21                                                                           58                                      1                                      152                                    2                                                                                                                 2     9
  GutlessWorm                              16                                                    3                                                                                                                                                                           92                                                                           658                                     48                                     705                                                                          10                                    32                                    86    9
  HOT:Sanger (Seawater)                    8                                                                                                                                                                                                                                 61                                     1                                     653                                     26                                     1000                                   11                                    1                                     11                                          9
  WesternChannelOMM (Seawater)                                                                                                                         5[^\*^](#TN2){ref-type="table-fn"}                                                                                    71                                     1                                     332                                     27                                     787[^\*^](#TN2){ref-type="table-fn"}   14[^\*^](#TN2){ref-type="table-fn"}   4                                     4[^\*^](#TN2){ref-type="table-fn"}    56    10
  MILOCO:454Shotgun (Seawater)             0.25[^\*^](#TN2){ref-type="table-fn"}   0.06                                                                                                                                                                                      85                                     1                                     275                                     16                                     842[^\*^](#TN2){ref-type="table-fn"}   10[^\*^](#TN2){ref-type="table-fn"}   0.19                                  1[^\*^](#TN2){ref-type="table-fn"}    71    11
  Washington Lake (Freshwater)             88                                      5                                                                                                                                                                                         126                                    14                                    1643                                    144                                    824                                    19                                    23                                    42                                    84    11
  BotanyBay:454 (Seawater)                 2[^\*^](#TN2){ref-type="table-fn"}                                                                          1[^\*^](#TN2){ref-type="table-fn"}                                                                                    33                                     1                                     414                                     19                                     352[^\*^](#TN2){ref-type="table-fn"}   1[^\*^](#TN2){ref-type="table-fn"}    1                                     1[^\*^](#TN2){ref-type="table-fn"}    54    11
  AntarcticaAquatic:Shotgun (Seawater)     1                                                                                          1                                                                                                                                      41                                     2                                     36                                      7                                      661                                    3                                     1                                     3                                     24    11
  YLAKE (Hot spring)                       0.25                                    6                                                                   1                                                                                                                     69                                     0.12                                  427                                     7                                      543                                    1                                     4                                     8                                     47    13
  AntarcticaAquatic:454 (Seawater)         0.44                                    0.04                                                                0.13                                                                                                                  60                                     1                                     250                                     15                                     783                                    4                                     4                                     4                                     52    13
  BisonMetagenome (Hot spring)             11[^\*^](#TN2){ref-type="table-fn"}                   2[^\*^](#TN2){ref-type="table-fn"}                    95[^\*^](#TN2){ref-type="table-fn"}   8[^\*^](#TN2){ref-type="table-fn"}                                              144[^\*^](#TN2){ref-type="table-fn"}   11[^\*^](#TN2){ref-type="table-fn"}   1242[^\*^](#TN2){ref-type="table-fn"}   47[^\*^](#TN2){ref-type="table-fn"}    679[^\*^](#TN2){ref-type="table-fn"}   6[^\*^](#TN2){ref-type="table-fn"}    28[^\*^](#TN2){ref-type="table-fn"}   83[^\*^](#TN2){ref-type="table-fn"}   44    13
  MILOCO:454 (Seawater)                    0.08[^\*^](#TN2){ref-type="table-fn"}                 1[^\*^](#TN2){ref-type="table-fn"}                    18[^\*^](#TN2){ref-type="table-fn"}   2[^\*^](#TN2){ref-type="table-fn"}      4[^\*^](#TN2){ref-type="table-fn"}      19                                     2                                     80                                      7                                      136[^\*^](#TN2){ref-type="table-fn"}   1[^\*^](#TN2){ref-type="table-fn"}    8                                     10[^\*^](#TN2){ref-type="table-fn"}   5     14
  GOS (Seawater)                           4[^\*^](#TN2){ref-type="table-fn"}      0.08          4[^\*^](#TN2){ref-type="table-fn"}                    19[^\*^](#TN2){ref-type="table-fn"}   0.16[^\*^](#TN2){ref-type="table-fn"}   0.16[^\*^](#TN2){ref-type="table-fn"}   92                                     5                                     685                                     53                                     902[^\*^](#TN2){ref-type="table-fn"}   4[^\*^](#TN2){ref-type="table-fn"}    3                                     6[^\*^](#TN2){ref-type="table-fn"}    5     16
  DayNight (Seawater)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   0
  SapeloIsland (Seawater)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               0
  No. of metagenomes                       24                                      5             4                                    1                13                                    3                                       4                                       40                                     16                                    42                                      34                                     43                                     22                                    17                                    22                                    29    

*The prevalence of the individual genes in the different metagenomes were normalized to the metagenomic library sizes, i.e., reported values are frequencies per billion bases sequenced. Virulence genes not detected in any of the metagenomic datasets were ply, vacA, stx1, stx2, hlgB, sipB, sipC, and ssaD. Please refer to Supplementary Table [1](#SM1){ref-type="supplementary-material"} for a description of the metagenomes*.

*Detection (BLASTN; E-value \< 10^−6^) of partial 16S/23S rRNA genes of clinical strains known to carry the virulence genes*.

The data set was autoscaled across metagenomes (i.e., the mean value of a variable across all metagenomes was subtracted from the data set (mean-centered) and the mean-centered data divided by the standard deviation) prior to PCA. Autoscaling was needed to allow each virulence gene to contribute equally to the PCA model, independent of the overall abundance of each gene. The model was validated by cross-validation (Martens and Martens, [@B14]).

Only metagenomes consisting of subsamples from similar environments were included in the analysis, hence, 10 metagenomes were excluded (MontereyBay transect, BotanyBay:Sanger, GuaymasBasin, MILOCO:454Shotgun, HOT:Sanger, AntarcticaAquatic:454, BotanyBay:454, AntarcticaAquatic:Sanger, MILOCO:454, GOS). Furthermore, one of the observed 16 virulence genes (*papH*) was not included the PCA as it was not found in any of the 34 metagenomes included in the analysis.

Results and discussion {#s3}
======================

Homologs of 16 of the 24 investigated bacterial human virulence genes were identified in the metagenomic datasets (Table [2](#T2){ref-type="table"}). The genes represented toxin, adhesin, secretion system, regulatory, inflammatory, and resistance virulence groups involved in urinary tract infections, gastrointestinal diseases, skin and foreign body infections, and systemic infections. To confirm that environmental homologs matched known virulence factors in the databases, we took the longest environmental homologs with the highest identity to the virulence genes and did a NCBI BLASTP search using the environmental sequences as queries. The environmental sequences matched relevant virulence groups in 95% of the cases (Supplementary Table [3](#SM3){ref-type="supplementary-material"}).

A positive relationship between numbers of detected genes and metagenome size was observed (Figure [1](#F1){ref-type="fig"}). Highest numbers of observed genes were seen in some of the largest metagenomes, but sequencing of more than 3--4 Gb generally did not result in the detection of further genes. Several of the investigated metagenomes were obtained by 454 pyrosequencing with average sequence reads of 100--400 bp. Some metagenomes were based on Sanger sequencing or on assembled sequences with average reads of 600--1000 bp (22). Since short reads (≤ 400 bp) tend to miss a significant amount of the BLAST homologs detected by long reads (Wommack et al., [@B35]), this may explain the high number of homologs detected in some of the relatively small Sanger-based metagenomes such as the "HydrothermalVent" (50 Mb) (Table [2](#T2){ref-type="table"} and Supplementary Table [1](#SM1){ref-type="supplementary-material"}).

![**Rarefaction plot showing numbers of detected virulence genus as function of size (numbers of bases sequenced) of metagenomic datasets**.](fmicb-07-01712-g0001){#F1}

The most frequently observed virulence gene homologs were the inflammatory O-antigen gene, *rbfE*, the regulator gene, *phoP*, and the Type III secretion apparatus gene, *invA*, which were detected in ≥85% of the metagenomes (Table [2](#T2){ref-type="table"}). These genes also showed the highest relative abundance within the individual metagenomes (Table [2](#T2){ref-type="table"}). Homologs to the regulator gene *phoQ*, the toxin gene *hlyA*, the fimbrial gene *papC*, the Type III secretion apparatus gene *spiA*, the inflammatory O-antigen gene *wbdI*, and the resistance genes *pqaB, yfbI* and *mecA* were also found widely distributed (28--72% of the metagenomes). On the other hand, homologs to the fimbrial genes *papA, papH, fimH*, and *focG* and the toxin gene *hylA* were only found in 2--10% of the metagenomes (Table [2](#T2){ref-type="table"}).

To assess if the origin of the virulence gene homologs was truly environmental or due to the presence of the clinically relevant bacteria from which the virulence factors were originally described, partial 16S or 23S rRNA gene sequences of *Salmonella* spp., *E. coli, Helicobacter* spp., *Streptococcus pneumoniae* and *Staphylococcus aureus* (Supplementary Table [2](#SM2){ref-type="supplementary-material"}) were targeted by BLASTN in the metagenomes. Although we specifically avoided metagenomes influenced by livestock and human discharges, relevant bacterial groups and their respective virulence genes co-occurred 66 out of the 313 times a virulence gene homolog was observed (Table [2](#T2){ref-type="table"}). Only *Salmonella* spp., *E. coli*, and *Helicobacter* spp. were observed. The presence of the bacterial groups may be explained by contamination during sample handling, survival and spreading from farms and wastewater treatment plants, or by the natural occurrence of environmental strains of the bacteria. We argue that the presence of the bacterial groups most likely is due to their natural occurrence. For instance, environmental survival of *E. coli* and *Salmonella* is well documented (Maugeri et al., [@B16]; Miller et al., [@B19]) and may be occurring in e.g., coastal environments. However, *E. coli* and *Salmonella* were also identified in oceanic, hot spring, and hypersaline environments (Table [2](#T2){ref-type="table"}), and albeit contamination by these bacteria cannot be ruled out, it does suggest that they naturally inhabit "unexpected" environments. If, on the other hand, it is assumed that the observation of clinically relevant bacteria is caused either by contamination or by livestock and human discharges, and all cases of dual presence of virulence genes and bacteria are excluded from the analysis, 80% of the detections of virulence gene homologs would still be unbiased by the presence of clinically relevant bacteria. Since 16S/23S rRNA genes contrary to virulence genes often are found in multiple copies per cell (Vos et al., [@B31]), the risk that 16S/23S rRNA genes of clinically relevant bacteria were overlooked in the metagenomic datasets is miniscule. Thus, the large majority of our data supports the existence of a natural environmental reservoir of genes normally associated with human pathogenic bacteria.

Distribution of virulence gene homologs across environmental metagenomes
------------------------------------------------------------------------

The bacterial virulence gene homologs were present in virtually all the metagenomes representing soil, seawater, freshwater, marine sediment, hot spring, the deep-sea, hypersaline mats, microbialites, gutless worms, and glacial ice environments (Table [2](#T2){ref-type="table"}). Hence, environments were characterized by highly contrasting environmental conditions in terms of pH, pressure, temperature, salinity, and nutrient availability. Moreover, the virulence gene homologs were globally distributed on six continents and in four oceans (Figure [2](#F2){ref-type="fig"}), i.e., all continents and oceans for which relevant metagenomes were available. The maximal number of different genes detected within one metagenome was 16 (Table [2](#T2){ref-type="table"}). Only in two metagenomes were no genes detected (Table [2](#T2){ref-type="table"}). These metagenomes originated from surface water off Hawaii, and coastal sea water off the Sapelo Island (USA) (Table [2](#T2){ref-type="table"} and Supplementary Table [1](#SM1){ref-type="supplementary-material"}). The fact that more than half of the sequences of the Hawaiian metagenome were cyanobacterial may explain the lack of virulence gene homologs in this metagenome. With respect to the Sapelo Island metagenome, homologs to virulence genes were observed in another and larger metagenome from Sapelo Island (the Sapelo2008 metagenome; Table [2](#T2){ref-type="table"}). Thus, the lack of virulence gene homologs in the SapeloIsland metagenome was most likely due to its small size.

![**Global richness map of putative bacterial human virulence genes in environmental metagenomic datasets**. Circles represent findings of genes, and areas of circles the number of virulence genes detected.](fmicb-07-01712-g0002){#F2}

Metagenomes representing comparable environments (e.g., seawater) did not share a set of typical or common virulence gene homologs. To discern relative differences between the metagenomes with respect to occurrence and frequency of the gene homologs, a principal component analysis (PCA) was performed. The variation explained by the first two principal components was relatively low (43%) yet useful for describing correlations between metagenomes and virulence gene homologs (Figure [3](#F3){ref-type="fig"}). No general trends were observed as the bulk of the metagenomes clustered relatively closely together with low, and in some cases, negative scores for both PC1 and PC2. However, the CellCapture and Bison metagenomes had relatively high positive PC1 and PC2 scores, while a high positive PC1 and a large negative PC2 score were observed for Washington Lake. The DeepMed metagenome, on the other hand, had a low PC1 score and a very large negative PC2 score. The high positive PC1 scores of the adhesion genes *papC* and focG, the resistance genes *pqaB* and *yfbI*, and the secretion gene *spiA*, indicated that these genes were relatively over-represented in BisonMetagenome and CellCapture as compared to the other metagenomes (Figure [3](#F3){ref-type="fig"} and Table [2](#T2){ref-type="table"}). The Washington Lake metagenome had above average gene frequencies of the regulatory gene *phoP*, the resistance gene *mecA* and the inflammatory gene *rbfE*, whereas the DeepMed metagenome had a relative over-representation of the toxin gene *hlyA* in comparison to the other metagenomes (Figure [3](#F3){ref-type="fig"} and Table [2](#T2){ref-type="table"}).

![**Bi-plot of principal component 1 (PC1; 26.4% explained variance) and 2 (PC2; 14.2% explained variance) of the PCA model**. Metagenomes are marked in blue, while virulence gene homologs are marked in red.](fmicb-07-01712-g0003){#F3}

Homologs to eight of the investigated virulence genes (*vacA, stx1, stx2, sipB, sipC, hlgB, ply, ssaD)* were not identified in any of the metagenomes. According to the databases, these genes have narrow phylogenetic host ranges, limited to only one or two bacterial classes/families (Søborg et al., [@B24]; unpublished data). This narrow phylogenetic host range, however, is unlikely to be the reason why homologs were not identified in the metagenomes since e.g., *phoQ* and *yfbI* also have rather limited host ranges (*Gammaproteobacteria*), yet were frequently observed (Table [2](#T2){ref-type="table"}). By PCR-based techniques we have previously observed *vacA, stx1*, and *stx2* in various soils and a freshwater biofilm (Søborg et al., [@B23]) demonstrating that these genes are present amongst non-human environmental microbiomes.

Potential ecological functions of virulence gene homologs
---------------------------------------------------------

The ecological function of the virulence gene homologs in environmental bacteria, and why certain homologs seem to be overrepresented in some microbiomes, is uncertain. The ubiquitous presence of the genes suggests they carry important functions, and it has been proposed that they are house-keeping genes essential for survival and proliferation in outer environments (Søborg et al., [@B23]) for instance by enabling the bacteria to cope with their physico-chemical environment, avoid predation and/or interact with host organisms. Bacterial DNA isolation procedures usually involve a size fractionation step of the environmental sample prior to DNA extraction. Thus, with the exception of the termite gut, the gutless worm and the Alvinella pompejana metagenomes it is not likely that DNA of bacteria associated with large multicellular organisms is represented in the metagenomes. This could suggest that the primary function of the virulence gene homologs is unrelated to interactions with multicellular host organisms. However, since free-living and host-associated subpopulations may coexist, one cannot rule out the possibility that the identified virulence gene homologs also play a role in host-pathogen interactions in the outer environment.

To be important for the survival of the bacteria in the environment the virulence gene homologs must be both functional and expressed in the non-host environments. Full-length or almost full-length (\>95%) sequences of seven genes (*papA, papC, focG, invA, phoP, rfbE*, and *wbdI*) were observed by single TBLASTN hits. In the case of the *wbdI* gene in the Washington Lake metagenome, 98% of the sequence, including the GDP-mannose binding site, was covered by a single TBLASTN hit (Supplementary Figure [1A](#SM4){ref-type="supplementary-material"}). In some cases (e.g., *phoP* and *wbdI*) complete or almost complete sequences were found in several metagenomes. For the genes *papH, spiA, phoQ, pqaB, yfbI* and *mecA*, 65--80% of the complete gene sequence were covered by single TBLASTN hits in at least one metagenome (exemplified by *pqaB*; Supplementary Figure [1B](#SM4){ref-type="supplementary-material"}). The fact that full-length genes were frequently observed suggests they are potentially functional and do not just occur as fragmented virulence genes within the microbiomes. As further support we have previously demonstrated *in vitro* expression of *fimH* in natural soil bacterial isolates and the expression of a virulent phenotype in soil bacterial isolates carrying homologs to the hemolysin gene, *hblA*, (Søborg et al., [@B24]). Some of the metagenomes (DayNight, Contaminated soil, MILOCO454, PacificOcean, and WesternChannel) were metatranscriptomes. Hence, the 15 different virulence gene homologs detected in these metatranscriptomes provide evidence of environmental transcription.

The BisonMetagenome represented samples collected along a thermal gradient (93--56⋅C) of a hot spring (Bison Pool) in Yellowstone National Park (Supplementary Table [1](#SM1){ref-type="supplementary-material"}). Theoretically, the bacterial community in Bison Pool may have been contaminated by bacteria from livestock as bison frequenting the pool during winter months may potentially have been in contact with livestock outside the borders of the park. However, since samples were taken at water temperature between 93 and 56⋅C, the risk of contamination is miniscule. Hot spring environments are often characterized by dense biofilms, hence, a relative overrepresentation of adhesin genes (*papC* and *focG*) is intuitively understandable, as fimbriae on the outer surface of the bacteria may play an important ecological role in biofilm formation (Lasaro et al., [@B11]). The resistance genes *pqaB* and *yfbI* also appeared to be relatively overrepresented in BisonMetagenome. They encode L-Ara4N transferases that add L-Ara4N moieties to lipid A of the cell membrane of Gram-negative bacteria. The modification of lipid A interferes with the recognition by the Toll-like receptor TLR-4 on macrophages/neutrophils and promote resistance toward cationic antimicrobial peptides in phagosomes (Baker et al., [@B1]; Trent et al., [@B27]; Miller et al., [@B18]). Thus, an over-representation of these genes may be explained as an adaptation to resist phagocytosis by thermophilic amebae which in some cases have been shown to be abundant in hot water environments (Tyndall et al., [@B29]; Reeder et al., [@B21]). In a similar manner, *spiA* may play a role for the bacteria in their interaction with amebae. In *Salmonella, spiA* encodes a Type III secretion system outer membrane protein, SpiA, involved in the transport of virulence proteins into host cells (Cirillo et al., [@B3]). Therefore, SpiA, in combination with the L-Ara4N transferases encoded by *pqaB* and *yfbI*, may be involved in the infection (internalization) of amebae.

Homologs to *papC, focG, pqaB, yfbI*, and s*piA* were also overrepresented in the CellCapture metagenome based on samples of Pacific deep-sea sediment (Supplementary Table [1](#SM1){ref-type="supplementary-material"}). A plausible explanation for an over-representation of adhesion genes in deep-sea sediment is not straightforward but the resistance genes and the Type III secretion gene may help the bacteria avoid predation by protozoa the same way as in hot springs.

The Washington Lake metagenome had above average frequencies of homologs to the regulatory gene *phoP*, the inflammatory gene *rbfE* and the resistance gene *mecA*. The two-component regulatory system phoP/Q, encoded by *phoP* and *phoQ*, governs adaptation to low Mg^2+^ and Ca^2+^ environments and to intramacrophage survival (Groisman, [@B6]; Tu et al., [@B28]), while *rbfE* encodes a perosamine synthetase involved in the synthesis of the O antigen of lipopolysaccharide in the outer membrane of Gram-negative bacteria (Iguchi et al., [@B8]). The relative over-representation of *phoP* and *rfbE* in the WashingtonLake metagenome suggests they provide specific selective advantages for the bacteria in these two environments. On the other hand, since they encode general house-keeping functions, their presence may not be exceptional. *mecA* encodes a peptidoglycan transpeptidase (with low affinity to ß-lactam antibiotics) involved in cell wall synthesis (Wielders et al., [@B32]). Although *mecA* may be considered a basic house-keeping gene, the finding of environmental homologs is noteworthy. *mecA* is highly conserved among clinical methicillin-resistant strains of *Staphylococcus aureus* (MRSA) (\>90% sequence identity between strains) and encodes a peptidoglycan transpeptidase unlike any of the peptidoglycan transpeptidases normally produced by *S. aureus* (\< 21% sequence identity) (Lim and Strynadka, [@B12]). It consists of two regions with homology to a penicillinase gene and a penicillin-binding protein gene, and has been hypothesized to have evolved by homologous recombination and transferred to MRSA by horizontal gene transfer from an unidentified species (Hiramatsu, [@B7]; Lim and Strynadka, [@B12]). Our data suggest that this unidentified bacterial species may have been of environmental origin.

The ecological role of the above average frequency of homologs to the toxin gene *hlyA* in the DeepMed seawater metagenome is somewhat perplexing. *hlyA* encodes a hemolysin (HlyA) involved in the lysis of (red blood) cells. However, marine opportunistic fish and shellfish pathogens such as Vibrios may produce hemolysins, which share common features with the HlyA hemolysin (Zhang and Austin, [@B36]; Ruwandeepika et al., [@B22]). Hence, the bacteria carrying homologs to *hlyA* may possibly be involved in the pathogenesis of marine animals. Another potential role of hemolysin could be as an antipredator defense strategy to subvert phagocytosis similar to the group B streptococcal β-hemolysin/cytolysin encoded by the *cylE* gene (Liu et al., [@B13]).

Phylogenetic relationships between virulence genes and environmental homologs
-----------------------------------------------------------------------------

Similarities \>50% between clinical and environmental sequences were observed for eleven genes (*hlyA, papC, fimH, invA, phoP, phoQ, pqaB, rfbE, wbdI, yfbI* and *mecA*) and for 6% of all hits. In four cases (*fimH, phoP, phoQ* and *yfbI*) were an identity of more than 90% observed between the environmental and clinical sequences (Supplementary Table [3](#SM3){ref-type="supplementary-material"}). The results concur with a previous study in which sequence similarities of 90--100% between environmental and clinical virulence genes were found (Søborg et al., [@B23]), and taken together emphasize that sequence conservation between environmental and clinical genes may be high.

Phylogenetic analyses underlined the close relationship between clinical genes and environmental homologs, and demonstrated relatively large sequence variations between environmental homologs with no clear relationship between environment and phylogeny (Figures [4A](#F4){ref-type="fig"}, [5A](#F5){ref-type="fig"} and Supplementary Figure [2](#SM5){ref-type="supplementary-material"}). For instance, the phylogenetic distance between the clinical sequence of WbdI and its most closely related environmental homolog (HypersalineMat metagenome), and the distance between YfbI and its most closely related environmental neigbor (GOS Sargasso Station 13 metagenome), were both 0.7 (Figures [4A](#F4){ref-type="fig"}, [5A](#F5){ref-type="fig"}). In contrast, phylogenetic distances up to 1.8 and 2.2 were observed between the environmental homologs to WbdI and YfbI, respectively. Furthermore, distances between the environmental sequences were larger than the largest distances between the clinical sequences and their environmental homologs (Figures [4A](#F4){ref-type="fig"}, [5A](#F5){ref-type="fig"} and Supplementary Figure [2](#SM5){ref-type="supplementary-material"}).

![**Minimum evolution distance tree (A) and maximum parsimony tree (B) of environmental sequences homologous to the protein sequence of *wbdI.*** Only environmental sequences with a coverage \>95% (\>141 aa) are included. The clinical sequence of wbdI is shown in red boldface type. Bootstrap values above 50% are displayed.](fmicb-07-01712-g0004){#F4}

![**Minimum evolution distance tree (A) and maximum parsimony tree (B) of environmental sequences homologous to the protein sequence of *yfbI***. Only environmental sequences with a coverage \>40% (\>220 aa) are included. The clinical sequence of yfbI is shown in red boldface type. Bootstrap values above 50% are displayed.](fmicb-07-01712-g0005){#F5}

Examination of the evolutionary history of the virulence genes indicated that the environmental homologs in some cases represented ancient lineages in the phylogeny of the clinical genes (Figures [4B](#F4){ref-type="fig"}, [5B](#F5){ref-type="fig"} and Supplementary Figure [2](#SM5){ref-type="supplementary-material"}). In the case of YfbI, all environmental sequences were evolutionary basal to the clinical sequence (Figure [5B](#F5){ref-type="fig"}), while several of the environmental homologs to WbdI, PapA and HlyA showed deep branching positions relative to the clinical sequences (Figure [4B](#F4){ref-type="fig"} and Supplementary Figure [2](#SM5){ref-type="supplementary-material"}). This suggests that the environmental genes in some cases evolved relatively early in the phylogeny of the clinical genes and that the origin of virulence genes in clinically relevant bacteria may be environmental.

Conclusion {#s4}
==========

Our findings highlight the widespread occurrence and transcription of genes often associated with bacterial human pathogens across diverse natural (non-human) environments. Environmental bacteria carrying homologs to virulence genes are probably not opportunistic human pathogens. Most likely, the homologs to the virulence factors are maintained because of the advantages they provide for the survival of the bacteria outside the host environment (Brown et al., [@B2]). Pathogenesis may, therefore, evolve by "coincidental" selection in non-host environments because traits that confer virulence in a host environment may also be beneficial for bacteria in non-host environments (Mikonranta et al., [@B17]). The high degree of sequence conservation between some of the environmental and clinical genes, and the phylogenetically deep branching positions of the environmental sequences, suggest that the environmental genes are precursors of virulence genes found in clinically relevant bacteria. Hence, our results suggests that coincidal selection in outer (non-host) environments may be a source of virulence traits for opportunistic pathogenic bacteria.

Author contributions {#s5}
====================

DS, NH, MK, and NK conceived and designed the research project and contributed with the interpretation of results. JC performed the statistical analysis. DS and NK wrote the manuscript which was approved by all authors.

Conflict of interest statement
------------------------------

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The work was supported by a grant \[645-08-0110\] from The Danish Agency for Science, Technology and Innovation to the Postgraduate Research School of Environmental Chemistry and Toxicology (RECETO). We wish to thank the staff of the CAMERA Portal for their assistance during the data analysis. In particular, we would like to express our gratitude to M. Chiu for his speedy response to all our many queries. G. Levin is thanked for his help with the preparation of the global richness map, and L. H. Hansen for his valuable comments on the manuscript. Finally, we wish to thank the reviewers for their constructive criticism that substantially improved the manuscript.

Supplementary material {#s6}
======================

The Supplementary Material for this article can be found online at: <http://journal.frontiersin.org/article/10.3389/fmicb.2016.01712/full#supplementary-material>

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

[^1]: Edited by: Frank T. Robb, University of Maryland, Baltimore, USA

[^2]: Reviewed by: Carl James Yeoman, Montana State University, USA; Marla Trindade, University of the Western Cape, South Africa; Aaron Conrad Ericsson, University of Missouri, USA

[^3]: This article was submitted to Evolutionary and Genomic Microbiology, a section of the journal Frontiers in Microbiology
