Indication of spatially random occurrence of Chlamydia-like organisms in Bufo bufo tadpoles from ponds located in the Geneva metropolitan area.

Occurrence of bacteria belonging to the order Chlamydiales was investigated for the first time in common toad (Bufo bufo) tadpole populations collected from 41 ponds in the Geneva metropolitan area, Switzerland. A Chlamydiales-specific real-time PCR was used to detect and amplify the Chlamydiales 16S ribosomal RNA-encoding gene from the tails of 375 tadpoles. We found the studied amphibian populations to host Chlamydia-like organisms (CLOs) attributable to the genera Similichlamydia, Neochlamydia, Protochlamydia and Parachlamydia (all belonging to the family Parachlamydiaceae), Simkania (family Simkaniaceae) and Estrella (family Criblamydiaceae); additionally, DNA from the genus Thermoanaerobacter (family Thermoanaerobacteriaceae) was detected. Global autocorrelation analysis did not reveal a spatial structure in the observed CLOs occurrence rates, and association tests involving land cover characteristics did not evidence any clear effect on CLOs occurrence rates in B. bufo. Although preliminary, these results suggest a random and ubiquitous distribution of CLOs in the environment, which would support the biogeographical expectation 'everything is everywhere' for the concerned microorganisms.

Introduction

The order Chlamydiales consists of strict intracellular bacteria that replicate within animal eukaryotic cells [1]. Molecular evidence suggests the existence of two lineages within Chlamydiales, which diverged between 700 and 1400 million years ago [2]: the family Chlamydiaceae and the Chlamydia-like organisms (CLOs) belonging to the families Piscichlamydiaceae, Clavichlamydiaceae, Simkaniaceae, Rhabdochlamydiaceae, Waddliaceae, Parachlamydiaceae, Criblamydiaceae and Parilichlamydiaceae [1], [3].

Chlamydiaceae were discovered in the 1960s and are associated with important diseases [1], [4]. Among the representatives of the family are Chlamydia trachomatis, aetiologic agent of the human visually impairing trachoma [1]; C. psittaci, responsible for lung infections in birds and humans [5] and equine infections [6]; C. abortus, causing abortion in domestic animals [4]; and C. pneumoniae, triggering several pathologies including respiratory diseases in humans, rhinitis in koalas and horses, as well as conjunctivitis in reptiles [1], [3], [4]; notably, C. pneumoniae was also reported in the African clawed frog, the great barred frog, the blue mountains tree frog and the common frog, which was also positive for C. abortus and C. suis [7], [8], [9]. Furthermore, Amphibiichlamydia ranarum was found at high prevalence in invasive bullfrog populations, possibly contributing to the amphibian biodiversity crisis [10].

CLOs were discovered in the 1980s, when Waddlia chondrophila was isolated from a bovine foetus [1], [3]. Most research has focused on identifying emerging pathogens among CLOs [3], [4], [11], with (a) W. chondrophila being suspected to trigger human miscarriage [12] and ruminant abortion [13], [14], (b) Parachlamydia acanthamoebae human miscarriage [15], (c) Simkaniaceae, Parachlamydiaceae and Rhabdochlamydiaceae respiratory diseases in humans and cattle [16], [17], [18], [19] and ocular infections in cats [20], (d) Parachlamydia species to concur to the massive mortality affecting an endangered midwife toad population [21] and (e) bacteria belonging to Clavichlamydiaceae, Parachlamydiaceae, Parilichlamydiaceae, Piscichlamydiaceae, Rhabdochlamydiaceae and Simkaniaceae having a recognized role in epitheliocystis [3], [4].

The majority of CLOs were first isolated from heterogeneous environmental sources [3]. Notably, Parachlamydiaceae, Simkaniaceae, Criblamydiaceae and Waddliaceae are observed as obligate endosymbionts of free-living amoebae of the genera Acanthamoeba and Hartmanella [7], [22], [23], [24], [25], which are therefore expected to convey CLOs dispersal [26]. Possibly reflecting the ecological tolerance and dispersal capabilities of such vectors, CLOs are commonly considered ubiquitous [1], [3], [4], [27], with a diversified set of hosts including mammals like humans and the fruit bat, marsupials, reptiles like chelonians, lizards and snakes, and fish like the leafy sea dragon and the Atlantic salmon [1], [3], [4], [28].

Ecologic studies have been conducted to elucidate possible links between environmental conditions and composition of protozoa communities and distributions, leading to the rejection of the ‘everything is everywhere, but the environment selects’ paradigm associated with free-living protozoa in some species [29]. Particularly, precipitation [30] and soil moisture, temperature, pH, dissolved oxygen and land cover showed an association with testate [31], [32], [33] and protosteloid amoebae [30] occurrence. Such evidences would suggest a possible environmental influence on CLO vectors and their endosymbionts at a local geographical scale.

Here, we investigated the presence of CLOs in common toad (Bufo bufo) tadpole populations from the Geneva metropolitan area (Switzerland) for the first time and tested the ‘everything is everywhere, but the environment selects’ principle on the basis of the occurrence patterns observed [34]. Furthermore, in a public health perspective, we derived human population density around sampling sites and studied a possible relationship with bacteria occurrence, given the CLOs' ability to infect humans from the environment [3], [18], [35].

Materials and methods

The protocol for this study was designed according to the regulations implemented in the state of Geneva wildlife law (RFaune M 5 05.01), and the corresponding fieldwork was approved by the Direction of Biodiversity for the state of Geneva under authorization 20140429/03 AS.

Sampling

In the context of the URBANGENE project, 145 ponds were identified and then inspected in the state of Geneva on the basis of the Marville (http://campus.hesge.ch/mareurbaine/) and the Karch (http://www.karch.ch/) databases, and also by means of a crowdsourcing campaign (http://urbangene.heig-vd.ch). Tadpoles were sampled from 9 to 22 April 2015 in a subset of 41 ponds (Fig. 1) that differed in size and degree of naturalness in the surrounding environments. Overall, 375 tadpoles were sampled (9.2 ± 2.7 tadpoles per pond, on average). In order to characterize the whole tadpole population in a pond, sampling privileged tadpoles coming from different frog spawns, whenever present; in such a case, tadpoles were collected shortly after they hatched from their frog spawn to reduce the chance of sampling siblings.

Fig. 1

Spatial representation of sampled ponds. Ponds are represented in red, with numbers corresponding to pond IDs in Table 1. Contextual information in background represents Geneva metropolitan area and highlights forests, agricultural, water and urbanized areas.

Table 1

Information reported for each pond

Pond ID	Pond name	Municipality	Longitude	Latitude	No. of inhabitants	N	Pa	OR	Genusb
1	Étang de l'ancien château de St-Victor	Avully	5.98	46.17	394c	10	6 (1)	0.600	?
2	Étang de la pépinière Jacquet	Satigny	6.05	46.22	107	13	4	0.308
3	Étang du signal de Bernex	Bernex	6.07	46.17	2038	12	3 (1)	0.250	Neochlamydia
4	Grand Etang des Mouilles	Bernex	6.08	46.19	149	10	4	0.400
5	Étang des Evaux	Confignon	6.09	46.19	1441	10	9	0.900
6	Étang de la ferme Lignon (PAF)	Vernier	6.09	46.21	2725	12	2	0.167
7	Étang Autrichien	Onex	6.11	46.18	4750	10	2 (1)	0.200	Simkania
8	Bassin du Parc Louis Bertrand	Lancy	6.12	46.19	6340	4	0	0.000
9	Étang Zimmermann	Lancy	6.12	46.19	5740	5	0	0.000
10	Étang Perfetta	Lancy	6.12	46.19	5829	4	1	0.250
11	Étang Lescaze	Genève	6.12	46.2	6393	10	3	0.300
12	Bassin du cimetière de St Georges	Genève	6.12	46.2	5963	9	2	0.222
13	Étang Spring	Lancy	6.12	46.19	5048	10	0	0.000
14	Bassin du Parc Chuit	Lancy	6.12	46.19	4609	10	4	0.400
15	Étang du Pré-d’œufeuf	Plan-les-Ouates	6.14	46.16	1574	10	6 (1)	0.600	Simkania
16	Étang Pinchat	Carouge	6.15	46.18	7000	15	2	0.133
17	Le Marais	Grand-Saconnex	6.11	46.23	7116	15	10 (1)	0.667	Protochlamydia
18	Étang du chemin des Prjins	Grand-Saconnex	6.12	46.23	3635	7	4	0.571
19	Bassin du château de Penthes	Pregny-Chambésy	6.14	46.23	467	5	0	0.000
20	Étang Berthier	Pregny-Chambésy	6.15	46.23	380	7	0	0.000
21	Étang du Jardin Botanique 2 (Serres)	Genève	6.15	46.23	2137	10	6	0.600
22	Étang du jardin botanique	Genève	6.15	46.23	2336	9	2	0.222
23	Étang Vieux-Clos	Chêne-Bougeries	6.18	46.18	911	12	10 (2)	0.833	Estrella; Thermoanaerobacter
24	Étang Flory	Chêne-Bougeries	6.18	46.19	946	7	1	0.143
25	Bassin Waldvogel	Chêne-Bougeries	6.18	46.19	2560	6	5 (2)	0.833	Estrella; Thermoanaerobacter
26	Bassin de la Station de Zoologie	Chêne-Bougeries	6.18	46.19	2641	11	2	0.182
27	Étang Richard	Chêne-Bougeries	6.18	46.19	1345	5	5	1.000
28	Étang Paradis Haake	Chêne-Bougeries	6.18	46.2	4661	8	2	0.250
29	Étang route de Chêne	Chêne-Bougeries	6.19	46.2	4342	8	1	0.125
30	Étang Broud	Chêne-Bougeries	6.19	46.2	3161	10	7 (3)	0.700	Parachlamydia
31	Étang Hutins	Chêne-Bourg	6.2	46.2	6043	8	7	0.875
32	Étang Loutan	Thônex	6.2	46.2	6349c	6	1	0.167
33	Étang de la clinique Bel-Air	Thônex	6.21	46.21	1543	9	1	0.111
34	Étang de Miolan	Choulex	6.21	46.23	318	8	2 (1)	0.250	?
35	Bassin du Centre horticole Lullier 2	Jussy	6.25	46.23	216	8	4 (1)	0.500	Similichlamydia
36	Étang des Dolliets	Jussy	6.28	46.23	23	12	6 (1)	0.500	?
37	Bois du Faisan amont	Versoix	6.15	46.28	436	10	1	0.100
38	Étang Bon-séjour	Versoix	6.16	46.28	3396	10	5 (1)	0.500	Simkania
39	Étang des Douves	Versoix	6.14	46.29	51	12	7	0.583
40	Étang Est Pré-Béroux	Versoix	6.14	46.3	3c	6	4	0.667
41	Étang de Combes-Chapuis	Versoix	6.12	46.3	13c	12	4	0.333

Information includes: experimental and actual name of ponds (pond ID and pond name, respectively), corresponding municipality, geographical coordinates, number of inhabitants within 1 km radius, number of B. bufo tadpoles sampled (N), number of Chlamydia-like organism–positive samples (P), observed occurrence rate (OR) and Chlamydiales taxonomic assignments at genus level.

In parentheses are number of samples (among positive ones) for which a sequence attributable to Chlamydiales was obtained.

Question marks indicate unsolved taxonomic assignations in samples for which a sequence attributable to Chlamydiales was obtained.

Reported number of inhabitants can be partial when pond is located close to French or to Canton of Vaud border (Fig. 1), for which demographic data were not considered.

DNA extraction

Sample preparation

After sampling, tadpoles were individually put in a water dish with Tricaine methane sulphonate (MS-222), causing rapid anaesthesia and death. The apical part of the tail was carefully clipped to investigate CLO occurrence, as this is one of the suggested options for retrieving tissue from amphibian juveniles [36]. After freeze-drying the tails, DNA was extracted at the LGC laboratories (Berlin, Germany), using the Sbeadex tissue kit (LGC, Teddington, UK) and following the manufacturer's instructions.

Pan-Chlamydiales real-time PCR assay

A Chlamydiales-specific real-time PCR [37] was used to detect and amplify the DNA fragment from 207 to 215 bp belonging to the Chlamydiales 16S ribosomal RNA–encoding gene. Quantification was performed using a plasmidic tenfold-diluted positive control tested in duplicate. Amplification reactions were performed in a final volume of 20 μL, containing: (a) iTaq Universal Probes Supermix with ROX (Bio-Rad, Reinach, Switzerland); (b) 0.1 μM concentration of primers panCh16F2 (5′-CCGCCAACACTGGGACT-3′) (the underlined bases representing locked nucleic acids) and panCh16R2 (5′-GGAGTTAGCCGGTGCTTCTTTAC-3′) (Eurogentec, Seraing, Belgium); (c) 0.1 μM concentration of probe panCh16S (5′-FAM [6-carboxyfluorescein]-CTACGGGAGGCTGCAGTCGAGAATC-BHQ1 [black hole quencher 1]-3′) (Eurogentec); (iv) molecular-biology-grade water (Five Prime, Hilden, Germany); and (v) 5 μL of sample DNA. Amplification started with an initial step of activation and denaturation at 95°C for 3 minutes, followed by 40 cycles at 95°C, 67°C and 72°C, each lasting 15 seconds, and was performed in a StepOne Plus real-time PCR system (Applied Biosystems, Zug, Switzerland). Samples with a threshold cycle value (C) of <35 were finally sequenced, as this is the observed limit for amplicon sequencing (S. Aeby and G. Greub, unpublished data).

DNA sequencing of the PCR-positive samples

According to the manufacturer's instructions, amplicons from positive samples were purified using the MSB Spin PCRapace (Stratec Molecular, Berlin, Germany). The sequencing PCR assay was performed using a BigDye Terminator v1.1 cycle sequencing kit (Applied Biosystems, Zug, Switzerland) and with specific inner primers panFseq (5′-CCAACACTGGGACTGAGA-3′) and panRseq (5′-GCCGGTGCTTCTTTAC-3′). Amplification was performed after an initial denaturation step at 96°C for 1 minute, followed by 25 cycles at 96°C for 10 seconds and 60°C for 4 minutes. Purification of the sequencing PCR products was done using the SigmaSpin Sequencing Reaction Clean-up (Sigma-Aldrich, Buchs, Switzerland), and the sequencing was performed in a 3130xL genetic analyser (Applied Biosystems). Sequences were blasted using the Geneious software [38].

Global spatial autocorrelation

After taxonomic assignments, the rates of CLO occurrence (ORs) were computed in each pond. A global spatial autocorrelation analysis was then conducted to investigate the presence of clusters or dissimilarities (i.e. the existence of spatial groups with similar ORs values or the tendency of similar values to stay apart in space) under the null expectation of a random distribution of occurrence. Moran scatterplots [39] were constructed testing different weighting criteria (i.e. the mean OR from the first two, four, six, eight and ten nearest ponds), and Moran's I [40] was estimated in each scenario as the slope of the linear regression between weighted and observed ORs. Observed and weighted ORs were centred before the analysis. Under the null hypothesis (h) of Moran's I = 0 and a significance threshold (α) of 0.05, statistical significance was derived by permuting ORs 9999 times, reestimating Moran's I in each permutation to obtain a reference distribution and computing the pseudo p-value associated with the observed Moran's I in each scenario. Analysis was performed using a script written in R [41].

Human population density around ponds

The number of inhabitants residing within 1 km radius from each pond (surface area: ∼3.14 km2) was derived from the Federal Statistical Office database for the Republic and Canton of Geneva (https://www.bfs.admin.ch/) for the year 2013. To test for the existence of a spatial relationship between human and CLO occurrence, a correlation test relying on Pearson's r was performed between the number of inhabitants and the observed ORs through the R function cor.test [41] (h: r = 0; α = 0.05). Ponds 1, 32, 40 and 42 were discarded, given the partial information about the number of inhabitants in the surrounding area.

Group comparison analysis

Thirty-two categories describing land cover were derived from the territorial information system in the Geneva database (https://ge.ch/sitg/sitg_catalog/geodataid/1133) at a 10 m resolution. Proportions were computed for each land cover category around the ponds as a function of a selected buffer. Particularly, buffers were defined by radii from 20 m to 3 km (total, 299), and the R function extract [42] was used to extrapolate the categories present in the buffers.

To group ponds sharing similar characteristics, hierarchical clustering was performed on the resulting proportions with the R function hclust [41] using the ‘average,’ ‘complete,’ ‘single,’ ‘Ward 1’ and ‘Ward 2’ clustering methods [43]. The Silhouette method [44] was used to identify the optimal number of clusters.

When two groups of ponds were identified, independent t tests or Wilcoxon rank sum tests were run with the R functions t.test or wilcox.test; one-way ANOVAs or Kruskal-Wallis rank sum tests were performed with aov or kruskal.test [41] in the presence of more than two clusters. Test choice was driven by first checking ORs for normality and homoscedasticity. P-values from the analyses performed with the same cluster method were corrected for multiple testing with the Benjamini-Hochberg method through the R function p.adjust [41]. Buffer scenarios determining at least one group with a single pond were discarded, given the impossibility of assessing normality.

Beta regression models

Association between land cover and ORs was also investigated with univariate beta regression, as implemented by the R function betareg [45]. To ease interpretation of regression coefficients (β), previously obtained land cover proportions were aggregated into five classes (Supplementary Table 1): ‘managed vegetation,’ referring to human-managed green areas; ‘vegetation near water,’ i.e. species communities well-adapted to live in close proximity with water; ‘forests,’ grouping forest species; ‘no vegetation’ encompassing urban environments as well as bare riverbeds; and ‘open fields,’ enclosing natural vegetation different from forests. Association tests (h0: β = 0; α = 0.05) were performed separately for each class and each buffer. Occurrence rates were transformed before the analysis to account for the presence of extreme values (i.e. OR = 0 or 1) [45]. The p-values from the tests involving the same aggregated class were corrected for multiple testing with the Benjamini-Hochberg method [41].

Results

We found 145 tadpoles (i.e. 38.7% of the samples, with the 95% credible region comprised between 33.7% and 43.9%; Supplementary Text) positive for Chlamydiales occurrence. Positive samples occurred across 36 ponds, with five ponds displaying no evidence of Chlamydiales occurrence. Occurrence rates spanned from 0 to 100% (Table 1; Fig. 2). Moran's I indicated the absence of significant spatial autocorrelation among ORs, regardless of the weighting scenario tested (Fig. 3). No significant relationship was observed between human presence and ORs (r = −0.145; p = 0.391) (Supplementary Fig. 1).

Fig. 2

Observed CLO occurrence rates over Geneva metropolitan area. Circles represent sampling sites (see Fig. 1), with size proportional to the number of tadpoles sampled and hue intensity following gradient in observed occurrence rates. CLO, Chlamydia-like organism.

Fig. 3

Global autocorrelation analysis results. Left column reports Moran scatterplots obtained using the first two, four, six, eight and ten nearest neighbours (i.e. ponds), respectively. Right column reports Moran's I reference distributions as obtained for each weighting scenario by permutation tests. Red vertical tick highlights position of observed Moran's I in reference distribution; grey vertical line is drawn to show I = 0 (i.e. null hypothesis). Red horizontal lines pinpoint percentiles 2.5 and 97.5 of reference distributions, underlining range of significant I values.

Only Wilcoxon and Kruskal-Wallis rank sum tests were performed due to non-normality and/or heteroscedasticity. After multiple testing correction, none of the land cover–based group showed any significant difference in ORs (Fig. 5). Likewise, none of the aggregated land cover variable showed a significant association with ORs after multiple testing correction (Fig. 6).

Fig. 5

Results of group comparison tests. P-values are reported on logarithmic scale after multiple testing correction, and as a function of both the radius used to characterize land cover around sampling sites and the clustering method used to classify ponds into environmental groups. In uppermost part of plot, dotted line indicates significance threshold (α = 0.05; thus, −log10α = 1.30). Line discontinuities depict tests not run (i.e. where at least one group was constituted by a single pond).

Fig. 6

Results of beta regression analysis. P-values associated with estimated regression coefficients are reported on logarithmic scale after multiple testing correction, and as a function of both the radius used to characterize land cover around sampling sites and the aggregated land cover category. In uppermost part of plot, dotted line indicates significance threshold (α = 0.05; thus, −log10α = 1.30).

Out of the 145 positive samples, 16 presented C < 35 and were subsequently sequenced at the 16S ribosomal RNA gene. Taxonomic attribution was possible at a family-level lineage for 13 samples. Particularly, six were found to be positive for Parachlamydiaceae, three for Simkaniaceae and two for Criblamydiaceae. The remaining two samples were positive for the genus Thermoanaerobacter (family Thermoanaerobacteriaceae). Among the six samples hosting Parachlamydiaceae, one was positive for the genus Similichlamydia (from pond 35), one for the genus Neochlamydia (pond 3), one for the genus Protochlamydia (pond 17) and three for the genus Parachlamydia (from pond 30, with two sequences displaying less than 1% divergence). All the samples infected by Simkaniaceae were assigned to the genus Simkania (from ponds 7, 15 and 38). Finally, samples positive for Criblamydiaceae and Thermoanaerobacteriaceae were assigned to the genera Estrella and Thermoanaerobacter, respectively, and were retrieved from pond 23 and 25 (Table 1; Fig. 4). As a result of low sequencing quality, no lineage could be identified for the samples coming from ponds 1, 34 and 36.

Fig. 4

Spatial occurrence of observed CLOs and Thermoanaerobacteriaceae. Red circles highlight ponds where each genus was found. Number of infected samples (i.e. tadpoles) is reported for each bacterial genus in brackets (Table 1) and refers to highlighted ponds (e.g. three tadpoles are positive for genus Parachlamydia from same highlighted pond). Grey area in background represents Geneva metropolitan area; size of circles is proportional to number of tadpoles sampled in each sampling site (see Fig. 2). CLO, Chlamydia-like organism.

Discussion

The order Chlamydiales comprises bacterial agents of important human and animal diseases, as well as emerging pathogens which affect a broad spectrum of hosts [1], [3], [4]. To our knowledge, the present study reports the first observation of CLOs occurrence in the species B. bufo.

Among the operational taxonomic units found are the genus Parachlamydia, which characterized the microbiome of a Pyrenean midwife toad population unable to recover from aggression by the fungus Batrachochytrium dendrobatidis [21]. Considering that co-occurrence of Chlamydiae and B. dendrobatidis was also observed in an Xenopus tropicalis population undergoing epizootic disease dynamics [9], an association was proposed between the skin microbiome of amphibians and B. dendrobatidis infection outcome [21]. Therefore, Parachlamydia occurrence might pinpoint a potential vulnerability for the host populations under study which should deserve attention for conservation.

The genera Simkania and Neochlamydia encompass recognized emerging pathogens. Indeed, Simkania species were associated with respiratory deficit in humans and epitheliocystis in fish and Neochlamydia species with ocular diseases in domestic cat and epitheliocystis [3]. There is weak evidence about Estrella involvement as a human pathogen, even if this bacterial genus remains understudied [46].

Such findings candidate B. bufo as a new host reservoir for CLOs, with implications for biodiversity conservation and possible consequences for public health. Nevertheless, no evident relationship seems to exist between the observed CLO occurrence rates and population density in the study area, although further studies are advisable to obtain more robust evidence in this regard; particular consideration should be accorded to Le Marais (pond 17) and Étang Hutins (pond 31), given their combination of high ORs and number of inhabitants (Supplementary Fig. 1).

Several studies support the paradigm of microbial biogeography ‘everything is everywhere, but the environment selects’ for Chlamydiales [1], [3], [34]. Here the absence of a clear spatial pattern and the ubiquitous occurrence of CLOs in the study area support the ‘everything is everywhere’ expectation, even at a local scale (Fig. 2, Fig. 3); nevertheless, no association was found with land cover typologies able to explain how ‘the environment selects’ CLOs distribution (Fig. 5, Fig. 6). In this regard, we believe further studies should focus on environmental variables with proved effects on amoebae distributions like precipitation, temperature, moisture, pH and dissolved oxygen [30], [31], [32], [33], [47]. Ideally, the characterization of CLOs niches might contribute in protecting biodiversity in host populations—and potentially the health of citizens.

Both the Chlamydia-like organisms' sequences and the R scripts developed for data analysis are deposited in the Zenodo public repository (https://zenodo.org) and are freely available online [48], [49].

Conflict of Interest

None declared.