Microbial community and antimicrobial resistance in fecal samples from wild and domestic ruminants in Maiella National Park, Italy.

This study aimed to provide new insights about antimicrobial resistance genes abundance and microbial communities of wild and domestic ruminants in wildlife-livestock interface. In total, 88 fecal samples were recovered from Apennine chamois, red deer, goat, cattle and sheep, and were collected in pools. The populations under study were selected based on ecological data useful to define sympatric and non-sympatric populations. Samples were screened for commonly used in farms under study or critically important antimicrobial resistance genes (aadA2, TetA, TetB, TetK, TetM, mcr-1). The microbial community composition was found to be different based on the species and land use of animals under study. Indeed, it was mostly characterized by phyla Firmicutes in bovine, Bacteroidota in chamois and Proteobacteria in red deer. Additionally, positive correlations between antibiotic resistance genes and microbial taxa (e.g., Tet genes correlated with Firmicutes and Patescibacteria) were described. Of the antimicrobials investigated, the abundance of mcr-1 gene suggests the importance of monitoring the wildlife in order to detect the emerging resistance genes contamination in environment. This study provides new data that highlight the importance of multidisciplinary and uncultured study in order to describe the spreading of antimicrobial resistance and related contamination in the environment.

Introduction

Antimicrobial resistance (AMR) is a multifaced issue considered one of the the 21st century most serious global threats [1]. The spread of AMR involves many habitats which include environments characterized by low or no human impact [1]. The mechanisms of resistance spreading are mostly related to mobile resistance genes and to selection and maintenance of multi-resistant bacteria in the environment based on the horizontal transfer [2]. Resistant microorganisms and their antimicrobial resistance genes (ARGs) may remain and eventually persist into environmental microbial community [3]. In this view the ARGs have been considered as emerging environmental contaminants [3]. Indeed, the 60–90% of antimicrobials used at farm level may be released in the environment through animal urine and feces [2], leading to the selection of resistance bacteria by means of their long-time persistence in environment.

Previous studies have mainly focused on the detection of ARGs starting from cultured selected bacteria, but it is noteworthy that most microbial species are non-culturable [4,5]. In this regard, alternative PCR- based approaches without preliminary culture have been suggested as a powerful method to detect potential antibiotic resistance in total microbial DNA [6]. This approach resulted mainly applied for soil, residues and water-related samples analysis [7,8] and during the last years has been also realized in wildlife [4,[9], [10], [11]]. Despite the wildlife has been suggested as an indicator for the spread of ARGs [6], there is non-univocal published data about its role in antibiotic resistance dynamics. Species living in different ecological niches are considered to have different roles in AMR dispersion. Indeed, the presence of AMR in wildlife results to be associated to habitat features, trophic characteristics, and anthropic pressure. Sympatric wild animals (that use anthropized habitats as feeding or refuge sites) have higher possibilities of harboring and spreading resistant bacteria [12]. Therefore, the wildlife has been used as sentinels of AMR and related ARGs in the environment and the species mostly investigated resulted wild birds and small mammals [13].

In this view, the aim of this study is to determine the fecal microbial community and occurrence of ARGs of wild and domestic ruminants with different levels of wildlife-livestock interface in the protected area of Maiella National Park (MNP), Italy.

Material and methods

Study area and sampling design

During October and November 2019, were collected a total of 22 fecal pools from red deer (Cervus elaphus), Apennine chamois (Rupicapra pyrenaica ornata) and extensive livestock, living in MNP (Italy). The Apennine chamois is a rare subspecies of chamois living in defined areas of Central Italy. The population size of the MNP is approximately 1300 individuals, and these data are the result of the reintroduction programs carried out in the past years. This species coexists in the territories of the Park with other widespread wild ungulates as the red deer population, of approximately 1500 individuals, and domestic cattle, sheep and goats traditionally raised in small farms by extensive grazing systems. The distribution of wild and domestic animals was previously determined by monitoring programs and georeferencing data useful to define the level of grazing lands sharing [14]. In this way, a group of 100 Apennine chamois sympatric with a farm of 120 goats, and 50 red deer coexisting with a farm of 300 sheep has been selected. Additionally, the non-sympatric populations composed by 70 cattle, 210 goats, 100 Apennine chamois and 20 red deer living in different areas of the Park were included (Fig. 1).

Fig. 1

Distribution of sympatric and non-sympatric animals investigated in grazing lands of the MNP [14].

As previously reported, the samples were collected gathering the freshly deposited feces and following the groups of grazing animals without disturbance of the animal activities [14]. The samples were stored at refrigeration temperature and analyzed by the laboratory of the Faculty of Veterinary Medicine of the University of Teramo within 3 h from the collection.

DNA extraction

Twenty-five grams of sample were homogenized in 250 ml of Phosphate Buffered Saline (PBS) and filtered using a sterile Filtra-bag (280 μm pore size) by Stomacher® for 2 min. After centrifugation at 2903 G for 10 min at 4 °C; total DNA was extracted from 200 mg of the resulting pellet using GeneAll Exgene™ Stool DNA mini kit. DNA quality and quantity were assessed using the spectrophotometer and fluorometer Denovix DS-11 FX (Wilmington, USA).

16S rRNA gene sequencing and data analysis

The V3 and V4 regions of 16S rRNA gene were amplified using the primer V3(5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG-3′) and V4 (5’-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC-3′) according to Klindworth et al. [15] with Illumina adapter overhang nucleotide sequences. The PCR products were pooled and indexed using the Nextera XT Index kit (Illumina, USA) following the 16S Metagenomic Sequencing Library Preparation Guide protocol (Illumina, USA). Libraries were sequenced using the Illumina MiSeq sequencing platform (San Diego, California, USA) with a 2 × 300 bp paired-end approach. The 16S rRNA data analysis was carried out using DADA2 package within the Quantitative Insights into Microbial Ecology 2 (QIIME2 version 2019.4) software [16,17]. The taxonomy categories were assigned by a Naive Bayes classifier and the q2-feature-classifier plugin, and identified using SILVA-Naive Bayes sklearn trained database [18]. The raw sequence reads have been deposited in the NCBI Short Read Archive under the accession number PRJNA832550.

The microbial community characterization and α- and β-diversity statistics were performed by the software Calypso (http://cgenome.net/wiki/index.php/Calypso) [19]. The microbial community composition and quantification within each sample group and among different animals were analyzed as previously described [1].

Determination of antimicrobial resistance genes (ARGs) by quantitative PCR (qPCR) analysis

The panel of 6 ARGs (aadA2, TetA, TetB, TetK, TetM, mcr-1) was selected based on the main antimicrobial classes used in livestock of study areas (tetracycline and streptomycin) along with the colistin, which is considered a critically important antibiotic for human medicine, and it has been particularly monitored in Europe by official surveillance program during the last years [20]. These genes along with the 16S rRNA gene were analyzed by SYBR Green® qPCR as previously described [7,9] in order to normalize the abundance of ARGs in fecal samples [21]. Each run was realized in triplicate. The copy number of target genes was calculated including an 8-point calibration curve obtained using ten-fold diluted positive controls and a good correlation coefficient (0.993 > R2 > 0.999). The normalized copy number of ARGs was calculated as the ratio of ARG copy number to 16S rRNA gene copy number [21]. The sensitivity of the technique was determined using a baseline threshold of −7. Any value ≤ − 8 was considered negative as described by [21].

Statistical analysis

The student's t-test of ARGs copies and Spearman correlation analysis to assess the association between the relative abundance of microbial taxa at phylum levels with ARGs, species or animal's group, were calculated using standard statistical software packages, Stata [22].

Results

General description of DNA sequences

Quality analysis and trimming returned a total of 8432 sequences with 8432 different features, with an average of 746,182 sequences per individual sample based on the analysis at operational taxonomy unit (OTU) level. The sequencing depth was good as indicated by the rarefaction curves for samples (Supplementary data 1).

Composition of bacterial communities

The microbial community composition of each sample was defined by means of the 16S rRNA sequencing analysis. At phylum level, member of the Firmicutes and Bacteroidota were dominant. In detail, Firmicutes phylum was highly abundant in bovine with the average of relative frequency of 56.3%, while a higher percentage of Bacteroidota (48%) respect to Firmicutes (38.7%) was identified in chamois. Proteobacteria phylum (26.3%) was more abundant in red deer samples, while Verrucomicrobiota (bovine: 5.2%; goat: 8.6%; sheep: 4.3%) was mostly detected in domestic animals (Fig. 2). The distribution at class level is consistent with the phylum frequency (Fig. 3). Indeed, the Bacteroidia class (Bacteroidota phylum) and Clostridia class (Firmicutes phylum) resulted prominent. Clostridia class (54.3%) was highly abundant in bovine, while a higher percentage of Bacteroidia (47.9%) respect to Clostridia (30.7%) was described in chamois. Gammaproteobacteria (Proteobacteria phylum) (30.7%) was more predominant in red deer, while Verrucomicrobiae (Verrucomicrobiota phylum) was mostly reported in domestic animals (bovine:3.4%; goat: 5.1%; sheep: 1.7%).

Fig. 2

Distribution of relative frequency of bacterial phyla in wild and domestic animals.

Fig. 3

Heatmap of microbial community distribution in wild and domestic animals at class level.

A strong positive correlation was observed between domestic species and phyla Cyanobacteria (r = 0.81p = 0.001), Fibrobacteria (r = 0.83 p = 0.001), Spirochaetota (r = 0.85 p = 0.001), Verrucomicrobiota (r = 0.77 p = 0.001), while a moderate correlation was described between domestic animals and Desulfomicrobiota (r = 0.66 p = 0.0005) or Elusimicrobiota (r = 0.66 p = 0.0005). Additionally, a moderate positive correlation resulted between wild animals and phyla Patescibacteria (r = 0.41 p = 0.04)) and Proteobacteria (r = 0.51 p = 0.001).

Concerning the group of animals, it was described a moderate positive correlation between sympatric animals and Proteobacteria (r = 0.48 p = 0.01) or Patescibacteria (r = 0.68 p = 0.0003), while it was reported a moderate positive correlation between non-sympatric animals and Firmicutes (r = 0.59 p = 0.003).

B: bovine; C: chamois; G: goat; D: red deer; S: sheep.

B: bovine; C: chamois; G: goat; D: red deer; S: sheep.

Bacterial community α-diversity

The α-diversity was evaluated at OTU level using Shannon index (Supplementary data 2) highlighting a significant difference (p < 0.01) between the species. The lowest alpha diversity was observed. in Apennine chamois and red deer while the highest diversity was noticed in bovine, goats and sheep. Finally, α-diversity appeared comparable between sympatric and non-sympatric populations.

Comparison among bacterial communities (β-diversity)

The analysis of β-diversity at OTU level was realized with principal coordinate analysis.

(PCoA) along with permutational multi- variable analysis of variance (PERMANOVA) showing a significantly difference (p < 0.05) of microbial community in sympatric animals compared to non-sympatric animals. This finding was supported by PCoA graphs, in which the Bray-Curtis distance was used as a distance metrics of β -diversity, and a difference in bacterial communities' distribution between sympatric and non-sympatric animal samples (Fig. 4) was reported. Additionally, the differences in microbial community composition between species were significant (p < 0.05) except for goat and sheep or chamois and red deer, and PCoA graph confirmed these observations (Fig. 4).

Fig. 4

β-diversity between sympatric or non-sympatric animals (A) and between species (B) by means of PCoA analysis according to Bray Curtis distances.

All investigated samples were found positive to at least one ARG.

The aadA2, TetA and TetK were detected in all the samples. The gene TetB was not amplified in samples from sympatric and non-sympatric goats and in non-sympatric red deer, while only one samples from sympatric sheep was resulted negative for TetM gene. Finally, the gene mcr-1 was not reported in one sample from sympatric chamois.

The total abundance of ARGs normalized to 16S rRNA ranged from 4.59 × 10−7 to 4.30 × 101 as reported in Table 1.

Table 1

Normalized copies of ARGs detected in each sample of wild and domestic animals.

Group	ID sample	Animal	aadA2	TetA	TetB	TetK	TetM	mcr-1
Sympatric	S1	Sheep	1.56 × 10−4	2.47 × 10−3	1.96 × 10−3	4.79 × 10−5	4.28 × 10−7	6.36 × 10−2
	S2	Sheep	1.36 × 10−4	7.97 × 10−4	2.24 × 10−3	1.52 × 10−5	Neg	9.38 × 10−3
	S3	Sheep	1.02 × 10−3	2.07 × 10−3	3.42 × 10−3	3.63 × 10−5	4.59 × 10−7	1.50 × 10−1
	D1	Red deer	3.65 × 10−4	8.35 × 10−3	3.62 × 10−3	8.24 × 10−6	3.54 × 10−7	3.05 × 10−3
	D5	Red deer	3.15 × 10−1	3.66 × 10−4	1.92 × 10−3	1.40 × 10−5	1.47 × 10−6	1.97 × 10−2
	D6	Red deer	4.72 × 10−4	2.24 × 10−3	4.47 × 10−4	2.25 × 10−5	1.79 × 10−7	6.42 × 10−4
	G1	Goats	2.31 × 10−4	3.16 × 10−2	Neg	2.26 × 10−5	8.14 × 10−6	5.52 × 10−1
	G2	Goats	1.75 × 10−4	4.14 × 10−2	Neg	1.99 × 10−5	2.74 × 10−6	5.11 × 10−2
	G3	Goats	9.93 × 10−5	2.44 × 10−2	Neg	1.16 × 10−4	6.41 × 10−6	2.85 × 10−1
	C1	Chamois	8.81 × 10−2	1.41 × 10−1	1.20 × 10−2	7.41 × 10−4	2.94 × 10−6	7.51 × 10−1
	C2	Chamois	1.67 × 10−3	1.87 × 10−2	1.01 × 10−3	1.24 × 10−4	1.13 × 10−5	Neg
Non-sympatric	D2	Red deer	4.50 × 10−4	2.40 × 10−2	Neg	2.79 × 10−4	2.22 × 10−5	4.83 × 10−1
	D3	Red deer	5.05 × 10−3	6.52 × 10−2	Neg	1.21 × 10−3	1.25 × 10−4	2.04 × 10−1
	D4	Red deer	1.02 × 10−3	3.22 × 10−2	Neg	2.72 × 10−4	2.95 × 10−5	1.29 × 101
	G4	Goats	1.13 × 10−3	2.87 × 10−2	Neg	3.15 × 10−5	7.58 × 10−5	4.47 × 101
	G5	Goats	7.98 × 10−4	1.29 × 10−2	Neg	6.97 × 10−5	3.03 × 10−4	5.53 × 10−1
	G6	Goats	4.38 × 10−4	1.85 × 10−2	Neg	5.42 × 10−4	2.69 × 10−4	2.16 × 101
	B1	Cattle	2.61 × 10−6	8.13 × 10−3	1.08 × 10−2	3.65 × 10−5	3.81 × 10−1	8.07 × 10−2
	B2	Cattle	4.12 × 10−7	1.64 × 10−2	9.51 × 10−3	1.02 × 10−4	2.74 × 10−1	3.44 × 10−1
	B3	Cattle	2.94 × 10−4	7.36 × 10−3	3.72 × 10−3	3.18 × 10−5	2.62 × 10−1	8.67 × 10−2
	C3	Chamois	6.10 × 10−8	1.43 × 10−1	1.22 × 10−2	3.89 × 10−4	4.30 × 101	5.73 × 10−1
	C4	Chamois	2.03 × 10−8	1.00 × 10−1	1.01 × 10−2	2.35 × 10−4	1.44 × 101	3.39 × 10−1
	C5	Chamois	1.55 × 10−6	3.33 × 10−2	1.63 × 10−3	7.83 × 10−5	5.93 × 10–1	4.46 × 10−2

Neg: negative sample was defined as any value ≤ − 8

The highest concentration of ARG was found in a sample of goat for mcr-1 (2.16x101copies/16SRNA copies) in the domestic animal's group and for TetM (4.30 × 101 ARGs copies/16SRNA copies) in Apennine chamois in wild animals.

Between the sympatric animals, the highest concentration of gene was identified in Apennine chamois for mcr-1 (7.51 × 10−1) and in non-sympatric animals in goat for mcr-1 (4.47 × 101 ARGs copies/16SRNA copies).

The statistical analysis showed that the means of ARGs of TetA (p = 0.04), TetK (p = 0.01), TetM (p = 0.001) and mcr-1 (p = 0.01) were highly significantly (p < 0.05) different in sympatric respect to non-sympatric animals. No significant difference of ARGs mean was identified between wild and domestic animals.

Co-occurrence between ARGs and bacterial taxa

The Spearman rank sum correlation analysis allowed to achieve different results. In detail, one moderate and two strong positive correlations resulted between phylum Patescibacteria and TetA (r = 0.46 p = 0.02), TetK (r = 0.65 p = 0.0006) and TetM (r = 0.66 p = 0.0005). Additionally, a moderate positive correlation was described between Firmicutes and TetK (r = 0.42 p = 0.04) or TetM (r = 0.40 p = 0.04). Considering the gene mcr-1, a moderate correlation with Myxococcota (r = 0.42 p = 0.04) was reported and finally a moderate correlation was described for Actinobacteriota and aadA2 (r = 0.55 p = 0.005), and Bdellovibrionota and TetA (r = 0.40 p = 0.04).

Discussion

This study represents the first analysis of microbial communities and ARGs in wildlife-livestock interface in Apennines, particularly relevant for the sampling design realized in a protected area and for the rare species involved such as Apennine chamois. Based on these features, the study did not allow to collect a large number of samples, but the research constitutes a preliminary step for further and more accurate investigations.

Similar methodologies were recently carried out on wildlife or environment in north-eastern Germany [11], in Poland [23] and Italy [24] but no ecological features (i.e., geographical distribution, size population and land use) of investigated animals were evaluated.

Additionally, Lanconi et al. [25] focused on the microbial community and AMR in agricultural soils fertilized with livestock manure in Italy, without considering the potential influence of wild species co-existing in the same territories.

The results obtained in our study showed that Bacteroidota and Firmicutes, along with Proteobacteria in wild ruminants, as reported by previous studies, are the dominant colonizing bacterial groups of the digestive system in domestic ruminants [26,27]. The Firmicutes/Bacteroidota ratio in gut microbiota appears to evolve during the life of animals and humans and it results in diverse digestive ability [28]. Several factors related to individual features and geographical localization may influence the microbiota in terms of quantity and diversity [28]. For examples, in dogs a food supplementation with fiber was associated to the increase of Firmicutes and Bacteroidota, and it was described as determinant for digestion of carbohydrates and other high molecular nutrients [28]. In wildlife Proteobacteria was described as most predominant phylum in giant panda and wild sika deer, resulting useful to optimize the degradation of lignin [29]. Finally, this phylum was also related to catabolism of various components of bovine fodder [30]. This evidence suggests that the composition of dietary may influence the relative abundance of these phyla in animals. Additionally, significant differences between wild and domestic animals and sympatric and non-sympatric populations were observed in our study. The correlation analysis showed results in line with previous studies. In detail, the association between phyla and domestic animals confirmed the common bacterial phyla of fecal microbiota in domestic ruminants, as previously described by other authors [31]. In addition, the correlation between Proteobacteria and wild animals was described according to the data previously reported [27]. Finally, the correlation between non-sympatric animals with Firmicutes and sympatric animals with Proteobacteria or Patescibacteria, may be related to the different activities and diets of animals. Indeed, other authors have associated Firmicutes to the feedlot animals and Patescibacteria to grazing animals [31].

These results suggest that the diversity of microbial community may be influenced by species (wild or domestic animals) and diets as assumed by other authors [28] and that the use and level of sharing of the land play a key role to define the composition of fecal microbiota.

Regarding the ARGs distribution, a significant difference of TetA, TetK, TetM and mcr-1 genes was observed between sympatric and non-sympatric animals, highlighting one more time the potential influence of the environment to define the ARGs profile. Indeed, different habitats may accumulate different ARGs [32].

Different correlations between phyla and ARGs have already been described. In detail, the Proteobacteria phylum was previously identified as particularly abundant of ARGs [33] and this evidence may explain the correlation detected in this study between ARGs and the proteobacterial classes, recently renamed as Myxococcota and Bdellovibrionota phyla [34]. The correlation between aadA2 and Actinobacteria reported in this study, were previously justified by other authors considering that many aminoglycosides have been isolated by bacteria of this phylum [35]. The correlation of Tet genes with Firmicutes and Patescibacteria was observed and appears to be in line with previous results that include the tetracycline genes in the top 15 ARGs accumulated in different habitats (i.e., hospitals, farms, soil, water, and wastewater) [32].

Tetracycline resistance-encoding genes were the most prevalent, as described in other similar studies in Brazilian wild birds [1], in tortoises in the Galapagos Archipelago, in Ecuador [9], in guignas [10] and Andean foxes in Chile [4], and in wild boar and foxes in Poland and Germany [11,23]. The high tetracycline resistance occurrence is not surprising, considering that these antibiotics have been widely employed for treatments in human and veterinary medicine as well as in agriculture [1]. Indeed, as reported by technical staff of MNP, this antibiotic was widely used in livestock of extensive farm situated in the protected area (Angelucci S. personal communication). Noteworthy, the 75–80% of the doses of tetracyclines are excreted in urine and feces with a long-term persistence of 578 days in soil that returns in a selection of resistant bacteria [36].

Therefore, the soil resistome derived from pathogenic and non-pathogenic bacteria present in the environment may consist of ARGs relevant to both human and veterinary medicine [36].

In this view, the colistin resistance gene mcr-1 gene resulted the most abundant target observed not only in domestic animals but also in wild sympatric group, highlighting the importance of studies focused on the environmental contamination of AMR. The worldwide presence of colistin resistance gene, previously reported in human [37], animals [38,39], animal food-products [40], and the environment [41], is an alarming indication of the inevitable progression to pan-drug resistance [37]. The mcr-1 gene evidence of this study confirms the important role of wildlife as possible sentinels of this trend in the human-wildlife-environment interface [[42], [43], [44], [45]]. Similar reports carried out on wild animals showed the mcr-1 in free-ranging seabirds, such as gulls in Europe [45], North America [46,47], Oceania [48], and South America [1,49], while in wild mammals phenotypic and genetic colistin resistance was reported in isolates recovered from fallow deer in Europe [50], barbary macaques in Africa [51] and Père David's deer in Asia [52]. In Central Italy, similar studies are available only for isolates from hunted wild boar and free-ranging Apennine chamois [20,42,53].

To date, the investigations regarding animal and environment context are mainly focused on selected pathogens and ARGs of public health concerns [54]. Conversely, the results of this study suggest that new measures to understand drug resistance dynamics are necessary. In this regard, the environmental monitoring, the risk assessment, and the quantification of AMR should be realized involving wildlife sources and by means of conventional microbiological analysis improved by culture-independent strategy [54].

Conclusions

The ARGs differences found in this study were linked to species' gut microbiome composition and ecological relationship between animals, providing for the first-time new data concerning this topic in protected area in Italy and at wildlife-livestock interface.

These preliminary results represent an attempt to overcome the limitation of culture-based methods and opportunistically sampling used to study AMR, applying a multidisciplinary approach targeted on the territory.

Funding

This work was supported by the Project “Demetra” (Dipartimenti di Eccellenza 2018–2022, CUP_C46C18000530001) funded by the Italian Ministry for Education, University and Research.

Ethical approval

Not required.

CRediT authorship contribution statement

Camilla Smoglica: Conceptualization, Methodology, Software, Data curation, Formal analysis, Software, Validation, Writing – original draft, Writing – review & editing. Simone Angelucci: Visualization, Investigation, Writing – review & editing. Muhammad Farooq: Data curation, Formal analysis. Antonio Antonucci: Visualization, Investigation, Writing – review & editing. Fulvio Marsilio: Writing – review & editing, Funding acquisition, Supervision. Cristina E. Di Francesco: Conceptualization, Methodology, Software, Software, Validation, Writing – original draft, Writing – review & editing.

Declaration of Competing Interest

None declared.