A Cross-Sectional Study of Dairy Cattle Metagenomes Reveals Increased Antimicrobial Resistance in Animals Farmed in a Heavy Metal Contaminated Environment

The use of heavy metals in economic and social development can create an accumulation of toxic waste in the environment. High concentrations of heavy metals can damage human and animal health, lead to the development of antibiotic resistance, and possibly change in bovine microbiota. It is important to investigate the influence of heavy metals in food systems to determine potential harmful effects environmental heavy metal contamination on human health. Because of a mining dam rupture, 43 million cubic meters of iron ore waste flowed into the Doce river basin surrounding Mariana City, Brazil, in 2015. Following this environmental disaster, we investigated the consequences of long-term exposure to contaminated drinking water on the microbiome and resistome of dairy cattle. We identified bacterial antimicrobial resistance (AMR) genes in the feces, rumen fluid, and nasopharynx of 16 dairy cattle 4 years after the environmental disaster. Cattle had been continuously exposed to heavy metal contaminated water until sample collection (A) and compared them to analogous samples from 16 dairy cattle in an unaffected farm, 356 km away (B). The microbiome and resistome of farm A and farm B differed in many aspects. The distribution of genes present in the cattle’s nasopharynx, rumen, and feces conferring AMR was highly heterogeneous, and most genes were present in only a few samples. The relative abundance and prevalence (presence/absence) of AMR genes were higher in farm A than in farm B. Samples from farm A had a higher prevalence (presence) of genes conferring resistance to multiple drugs, metals, biocides, and multi-compound resistance. Fecal samples had a higher relative abundance of AMR genes, followed by rumen fluid samples, and the nasopharynx had the lowest relative abundance of AMR genes detected. Metagenome functional annotation suggested that selective pressures of heavy metal exposure potentially skewed pathway diversity toward fewer, more specialized functions. This is the first study that evaluates the consequences of a Brazilian environmental accident with mining ore dam failure in the microbiome of dairy cows. Our findings suggest that the long-term persistence of heavy metals in the environment may result in differences in the microbiota and enrichment of antimicrobial-resistant bacteria. Our results also suggest that AMR genes are most readily detected in fecal samples compared to rumen and nasopharyngeal samples which had relatively lower bacterial read counts. Since heavy metal contamination has an effect on the animal microbiome, environmental management is warranted to protect the food system from hazardous consequences.

The use of heavy metals in economic and social development can create an accumulation of toxic waste in the environment. High concentrations of heavy metals can damage human and animal health, lead to the development of antibiotic resistance, and possibly change in bovine microbiota. It is important to investigate the influence of heavy metals in food systems to determine potential harmful effects environmental heavy metal contamination on human health. Because of a mining dam rupture, 43 million cubic meters of iron ore waste flowed into the Doce river basin surrounding Mariana City, Brazil, in 2015. Following this environmental disaster, we investigated the consequences of long-term exposure to contaminated drinking water on the microbiome and resistome of dairy cattle. We identified bacterial antimicrobial resistance (AMR) genes in the feces, rumen fluid, and nasopharynx of 16 dairy cattle 4 years after the environmental disaster. Cattle had been continuously exposed to heavy metal contaminated water until sample collection (A) and compared them to analogous samples from 16 dairy cattle in an unaffected farm, 356 km away (B). The microbiome and resistome of farm A and farm B differed in many aspects. The distribution of genes present in the cattle's nasopharynx, rumen, and feces conferring AMR was highly heterogeneous, and most genes were present in only a few samples. The relative abundance and prevalence (presence/absence) of AMR genes were higher in farm A than in farm B. Samples from farm A had a higher prevalence (presence) of genes conferring resistance to multiple drugs, metals, biocides, and multi-compound resistance. Fecal samples had a higher relative abundance of AMR genes, followed by rumen fluid samples, and the nasopharynx had the lowest relative abundance of AMR genes detected. Metagenome functional annotation suggested that selective pressures of heavy metal exposure

INTRODUCTION
Despite the substantial advancements brought by industrialization, it is well known that incorrect management of industrial activity has the ability to cause severe damage to environmental, human, and animal health. Heavy metals are bio-accumulative elements that naturally occur in water and soil (Guilherme et al., 2005), but human activities such as mining can result in the environmental accumulation of heavy metals as toxic waste. Heavy metal contamination is a public health concern due to the presence of these toxic elements in food crops (Rai et al., 2019), milk (Zhou et al., 2019) and meat . Chronic heavy metal intake of foods derived from plants grown or animals raised in a contaminated environment can lead to bio-accumulation in the organism and multiple health risks (Jaishankar et al., 2014). Gastrointestinal cancer (Yuan et al., 2016), fetal growth restriction (Sabra et al., 2017), and neurological diseases (Chen et al., 2016;Karri et al., 2016;Bjorklund et al., 2018) are examples of heavy metal-related consequences in humans.
Increased heavy metal concentration in the environment has been linked to increased prevalence of antibiotic resistance in bacteria by a co-resistance phenomenon (Nies, 1992;Silver, 1996;Wireman et al., 1997;Choudhury and Srivastava, 2001;Baker-Austin et al., 2006). Co-resistance is characterized by the closeness between different types of resistance genes located in the same genetic element, such as plasmids and transposons, which are often transferred together. As a result of this connection, the transfer of one gene (heavy metal resistance) may occur in concert with the transfer of the closest gene (antibiotic resistance) (Baker-Austin et al., 2006). Consequently, some resistance mechanisms are shared between antibiotics and heavy metals (Silver and Phung, 1996;Levy, 2002;Mukhopadhyay and Rosen, 2002;Roberts, 2005). Efflux pumps are another important resistance system in both gram-negative and gram-positive bacteria. These transmembrane proteins are responsible for pumping antimicrobial (molecules/compounds/substances) (e.g., antibiotics, heavy metals, biocides) out of the cell, regulating their concentrations in the bacterial internal environment (Blanco et al., 2016). Multi-drug resistance efflux pumps are important examples of cross-resistance determinants, and the resistancenodulation-division (RND) superfamily is of major clinical significance in gram-negative that transports heavy metals, and hydrophobic drugs (Zgurskaya and Nikaido, 2000).
Chronic heavy metal intake may lead to health complications due to its bio-accumulative effect (Jaishankar et al., 2014). Heavy metals are absorbed in the gastrointestinal tract but some also remain in the intestinal lumen, potentially interfering in the microenvironment and its functioning (Breton et al., 2013). For example, rodent models have shown that exposure to lead and cadmium changed the gut microbiota of rats; specifically, the relative abundance of Lachnospiraceae decreased and Lactobacillaceae and Erysipelotrichaceacae increased, resulting in gut dysbiosis associated with colitis and gut inflammation (Lepage et al., 2011;Breton et al., 2013). Bacteria unique to ruminants have been shown to be sensitive to heavy metals in vitro. An in vitro study of rumen environments showed that Bacteroides succinogenes, Ruminococcus albus, Bacteroides amylophilus, and Eubacterium ruminantium were sensitive to heavy metals, such as cadmium, mercury, copper, arsenic, and chromium (Forsberg, 1978). It is well-established that the rumen resident microorganisms are responsible for feed fermentation and production of important metabolites that maintain ruminant homeostasis; thus, it is important to evaluate the impact of heavy metal exposure on ruminal microbiome.
Anthropogenic activity (e.g., mining) is responsible for the deposition of heavy metals in the environment (Nguyen et al., 2019). In 2015, the rupture of a Brazilian mining dam resulted in iron ore waste flowing into the Doce river basin, which resulted in more than 43 million cubic meters of iron ore waste flowing into the Doce river basin surrounding Mariana City in the state of Minas Gerais (Samarco, 2016;do Carmo et al., 2017). Analysis of Doce river water 2 years after the event revealed aluminum, nickel, arsenic, copper, manganese, and lead higher than the maximum value allowed by Brazilian legislation (de Carvalho et al., 2017). Little is known about the effects of heavy metal environmental contamination on food production animals and their metagenomes. This presented the opportunity to investigate the consequences of long-term exposure to heavy metal contaminated drinking water on the dairy cattle metagenome. Thus, we performed a cross-sectional study to compare the prevalence and scope of antimicrobial resistance (AMR) genes between a contaminated farm and an unaffected farm 356 km away, in addition to characterizing the impact of prolonged exposure on the nasopharynx, rumen fluid, and fecal microbiomes.

MATERIALS AND METHODS
This research was conducted from August 2018 through November 2019 (sample collection in farm A and farm B was performed in August 2018 and October 2018, respectively). Sequencing and analysis of the metagenome were performed from August to November 2019. Sixteen dairy cows from a potentially heavy metal contaminated environment (farm A) and 16 dairy cows from a non-contaminated environment (farm B) were sampled at three body sites. The research protocol was reviewed and approved by the Ethics Committee in the Use of Animals of the School of Veterinary Medicine and Animal Science of the University of São Paulo (protocol no. 6424070217).

Animals and Facilities
A convenience sample of 32 dairy cows was obtained from two different farms located in Minas Gerais within 526 km of each other and were enrolled in this study.
Sixteen dairy cows were selected from a dairy farm located in Mariana City (Latitude: 20 • 22 41 S, Longitude: 43 • 25 0 W), Minas Gerais State, Brazil (farm A) in a potentially heavy metal contaminated area (de Carvalho et al., 2017). Cows were 3 years old or older to be enrolled in the study, ensuring that they were present at the farm during the environmental accident. Cows were housed in a single pen, grazed on pastures, and drank the affected river water. In addition to pasture, cows received mineral salt, wheat bran, and cornmeal. All animals were vaccinated for rabies and foot and mouth disease.
Sixteen cows were selected from a dairy farm located in Caldas City (Latitude: 21 • 55 23 S, Longitude: 46 • 23 15 W), Minas Gerais State, Brazil (farm B). The cows grazed on pastures and drank water from a non-contaminated lagoon. In addition to pasture, cows received corn and soybean meal as well as mineral salt. All animals were vaccinated for leptospirosis and foot and mouth disease.
Samples were collected from females on both farms. Both owners referred to cases of mastitis and reproductive disorders during the year before sampling time. The cows received antibiotic treatment on the farm when necessary. Unfortunately, health records were not available for every animal.

Data Collection
Two veterinarians from our research group performed physical examinations on each cow, assessing vital parameters (respiratory and heart rates, rectal temperature, and ruminal movements) and any clinical sign that indicated sickness. Animals displaying any symptoms of illness were removed from the study. Decreased milk production was also reported by the farmer for animals in the affected farm.

Sample Collection
Samples were collected from the nasopharynx, rumen, and rectum. Deep nasopharyngeal swabs (DNS) were collected using an 80 cm sterile swab (Provar, Brazil) covered by a thin sterile plastic sheath. The sample was obtained from the left nostril only. Mucus and dirt were removed from the nostrils with a paper towel, followed by inserting the swab into the nasopharyngeal cavity until met with resistance. The swab was passed 10 times over the respiratory mucosa before being removed from the nostril. The swabs were placed in cryogenic tubes and immediately stored in liquid nitrogen. Fecal samples were obtained by introducing a sterile swab into the rectum and rotating it for 15 s. Each swab was immediately stored in liquid nitrogen in a cryogenic tube. Five milliliters of ruminal fluid was obtained by introducing a plastic probe into the mouth until it reached the rumen. Using a vacuum system, the content was collected in a plastic bottle and immediately stored in liquid nitrogen. The system was cleaned with water and 70% alcohol between animals.

DNA Extraction
Before DNA extraction, all swab samples were pre-processed to obtain the maximum DNA concentration. Briefly, one milliliter of sterile PBS buffer 1× (pH: 7.4) was added to all swab samples and vortexed for 30 min. The swabs were removed, and samples were centrifuged at 3,000 × g for 30 min. The supernatant was discarded, and the pellet was resuspended in 300 µl of PBS. Rumen fluid samples were thawed and vortexed before the extraction. DNA extraction was performed using the MagMAX TM CORE Nucleic Acid Purification Kit (Thermo Fisher Scientific Inc., Waltham, MA, United States), according to manufacturer instructions using 200 µl of starting sample (concentrated swab material or rumen fluid). The DNA concentration was measured by a spectrophotometer (Nanodrop, Thermo Fisher Scientific Inc., Waltham, MA, United States).

Metagenome
Before the library preparation, the DNA concentration of each sample was also assessed by Qubit Fluorometric Quantification (Thermo Fisher Scientific, San Jose, CA, United States). Ten nanograms of each Qubit quantified genomic DNA was sheared with a Covaris E220 instrument operating SonoLab v6.2.6 generating approximately 300 bp DNA fragments according to the manufacturer's protocol. Between 10 and 100 ng of fragmented DNA was processed into Illumina compatible sequencing libraries using sparQ DNA Library Prep Kit (Quantabio, Beverly, MA, United States). Each library was barcoded with unique dual index sequences (NEXTFLEX R Unique Dual Index Barcodes, BioO Scientific). Library size and amount were confirmed with a Bioanalyzer High Sensitivity DNA chip. Polymerase chain reaction primers and reagents included in the sparQ kit were used to perform PCR, and products were purified with AMPure XP beads. Equimolar libraries were pooled and subjected to Illumina NovaSeq 6000 sequencing at 2 × 150 bp (Illumina, San Diego, CA, United States). Shotgun whole metagenome sequencing was performed at the Genome Sciences and Bioinformatics Core at the Pennsylvania State University College of Medicine, Hershey, PA, United States. Illumina bcl2fastq (released version 2.20.0.422) was used to extract de-multiplexed sequencing reads.

Data Analysis
Microbiome Fastq files were uploaded to the MG-RAST server (Keegan et al., 2016), concatenated and submitted to their standard pipeline to determine the relative abundance of the microbiota. In the MG-RAST pipeline, sequences were subject to quality control, including dereplication, removal of host-specific species sequences (Bos taurus, UMD v3.0), ambiguous base filtering, and read length filtering. Fragments were mapped against a nonredundant database (N5nr). Organism abundance was analyzed with a maximum e-value of 1 × 10 −5 , minimum identity cutoff of 60%, and a minimum alignment length cut-off of 15 (default values).
The normality of the data was assessed by the Shapiro-Wilk test (Shapiro and Wilk, 1965). Variances were analyzed by Bartlett's test (Bartlett, 1937). Relative abundance was calculated by dividing the number of reads of each taxa by the total number of reads of all taxa. Relative abundances of each phylogenetic level (phylum and genus) were analyzed by the Wilcoxon ranksum test (Wilcoxon, 1945) according to herd (farm A and farm B). Results were described as mean (M) and standard error of the mean (SEM). Differential abundance analysis was performed at the genus level to compare farm A and farm B within each anatomical site with the DESeq2 package in R. The analysis was performed with Wald significance test and mean dispersal type, and P values were adjusted with False Discovery Rate. All statistical analyses were completed in R software version 3.6.0 (RStudio Team, 2020). Significance was accepted at P < 0.05.

Resistome
Read quality was visualized on a random subset of 10 samples with FastQC (Andrews, 2010). Reads were trimmed with Trimmomatic for a minimum length of 50 bp (Bolger et al., 2014). Paired-end reads were merged in FLASH (Magoc and Salzberg, 2011) and converted to FASTA format in Seqkit (Shen et al., 2016). Reads were aligned to a custom AMR gene database built in BLAST (Altschul et al., 1990). The BLAST database was created from the publicly available FIGURE 1 | Hierarchy of antimicrobial resistance genes. Each gene or operon ("Group" level) belongs to a Mechanism, a Class, and a Type. Pairwise comparisons between farms and anatomical sites were made between Types, Classes, Mechanisms, and Groups to determine resistome differences.
MEGARes 2.0 AMR database 1 , which incorporates genes from CARD, ARGANNOT, AMRFinder, NCBI, ResFinder, BacMet, Lahey, and ResFinder for a total of 7,868 genes or operons that confer metal, biocide, drug, or multi-level resistance (Doster et al., 2019). The BLAST alignment resulted in a tabular list of reads that aligned to resistance genes or operons, hereon referred to as "gene counts" (Figure 1). Contaminant genes were detected and removed with the decontam R package based on negative control samples from DNA extraction (n = 2) and PCR (n = 1) (Davis et al., 2017). A mock community sample was used as a positive control at both DNA extraction and PCR, and both underwent sequencing for quality control purposes (Zymo Research, United States).
To account for bias, gene counts were normalized both to gene length and pseudo-sequencing depth. Gene lengths varied for the different mutations found in each gene, so gene counts were divided by the minimum gene length and then the number of merged reads per sample. The normalized data was scaled to relative abundance where the sample with the highest normalized gene count was set to 1. After scaling to relative abundance, we realized the data was strongly zero skewered and the distribution was not corrected by variable transformation. Therefore, we also analyzed the data in terms of gene presence, where data was converted to binary presence per sample (1 if the gene was present in the sample, 0 if not) as this represented the heterogeneity of the data. This is called "gene prevalence" throughout. Higher gene prevalence is defined as higher per-sample average gene presence; total gene presence per farm divided by the number of samples in each farm.
Non-parametric tests were used for pairwise comparisons of relative gene abundance. Significance thresholds were set at a False Discovery Rate corrected P < 0.05 (Benjamini and Hochberg, 1995). Separate comparisons were made for resistance type, class, mechanism, and group (Figure 1). Comparisons of relative gene abundance between farms were made with the Wilcoxon Rank-Sum test with False Discovery Rate adjustment for multiple comparisons. Comparisons of relative gene abundance among anatomical sites at each farm were made using the Kruskal-Wallis test with False Discovery Rate adjustment for multiple comparisons.

Metagenome Functional Annotation
To determine whether farm A was enriched in functions that promote resistance to both antimicrobials and heavy metals, we utilized mi-faser, a high-precision metagenome analysis pipeline which maps sequencing reads to molecular functions encoded by the read-corresponding genes, to annotate our data as a set of molecular functions (Zhu et al., 2018). Pairedend fastq files were uploaded to the mi-faser submission site 2 whose pipeline performed quality control with fastp (phredquality = 20, readlength = 40) before reads were mapped to the mi-faser extended reference database. Job output files including enzyme catalog (E.C.) numbers, functional annotation, and read counts per annotation were exported as .csv for all 90 samples. A Wilcoxon Rank-Sum test was performed to compare the number of unique reads identified between farms by sampling sites. Reads were then normalized according to the total read count in each individual metagenome and these data together with their metadata (Sample ID, farm, sampling site) were imported into the open-source machine learning software Waikato Environment for Knowledge Analysis (WEKA) (Frank et al., 2016). Within WEKA, the "InfoGainAttributeEval" attribute evaluator algorithm and search method "Ranker" were used to calculate the relative worth of each E.C. by measuring the information gain with respect to three nominal classes (1) "farm" (two levels), (2) "sampling site" (three levels), and (3) a composite class of "farm by site" (six levels). The data set was reduced by selecting the top 10 ranked E.C.s from each approach for pathway mapping in the Kyoto Encyclopedia of Genes and Genomes (KEGG) Mapper (Kanehisa and Sato, 2020).

RESULTS
Thirty-two dairy cows were selected from two dairy farms, one of which had heavy-metal contaminated drinking water. At the time of sampling, no clinical manifestations that indicated any disease was observed in any animal. Farm A's owner reported to decreased in milk production after the dam failure (250 L milk/per to 200 L milk/per day). No differences were noted by farm B's owner.
Rumen fluid, fecal swabs, and DNS were collected from each cow, totaling 96 samples. Six samples from farm B (one fecal swab, three rumen fluids, and two DNS) did not yield enough DNA concentration and were not considered for sequencing, so a total of 90 samples were sequenced and analyzed, in addition to appropriate negative and positive controls. Detailed sequencing data on total number of sequences and base pairs, qualified reads per sample and average length of sequences are presented in Table 1.

Resistome
A total of 549 AMR genes (called "groups" by the MEGARes database as it can include either genes or operons; called "gene" throughout) were analyzed after removing contaminant gene alignments (n = 74) and genes found only in one sample (n = 117). The distribution of AMR gene presence was highly heterogeneous between samples, as most of the 549 genes were found in only a few samples (Supplementary Figure S3).
Farm A had a higher relative abundance and higher resistance prevalence in all four AMR classes of biocides, drugs, metals, and multi-compound than farm B (Figure 4). Higher prevalence of biocide resistance in farm A was driven by multi-biocide resistance ABC (normalized per-sample average gene presence 1.375 farm A and 0.628 farm B) and RND (normalized persample average gene presence 3.5 farm A and 1.581 farm B) efflux pumps. Metal resistance prevalence was driven by tellurium resistance and multi-metal resistance, notably a multi-metal resistance RND efflux pump (per-sample average gene presence 0.438 farm A and 0.047 farm B). Significant multi-compound resistance genes that were denoted as multi-compound resistance for drugs and biocides included two RND efflux pumps and one Major Facilitator (MFS) efflux pump, which were almost absent from farm B samples; RND efflux pump "MEXY" was present in only 1 sample in farm B, while RND efflux pump "MTRD" and the MFS efflux pump "KDEA" were completely absent from farm B. Several types of drug resistance were more abundant in farm A than farm B (Supplementary Figure S4). Mechanisms conferring resistance to aminoglycosides, betalactams, glycopeptides, macrolide-lincosamide-streptogramin B (MLS), mupirocin, mycobacterium, nucleosides, tetracyclines, and multi-drug resistance were all more highly abundant and prevalent in farm A than farm B.
Within each farm, there were strong differences in AMR abundance and prevalence among anatomical body sites and the differences were consistent between farms. Generally, in both farms fecal samples had higher gene presence than rumen fluid samples, and DNS had the least amount of gene presence (Figure 5). Resistance types of biocides, drugs, metals, and multi-compound resistance differed between body sites on both farms. Interestingly, out of the 549 AMR genes measured via Kruskal-Wallis, 289 genes had significant differences in relative abundance among body sites on farm A, and 271 genes differed on farm B. Of those, only 75 on farm A and 20 on farm B had higher gene prevalence in either rumen fluid or DNS samples than the fecal samples, confirming that AMR prevalence is generally higher in fecal samples.

Metagenome Functional Annotation
The median number of unique functions identified by the mifaser for each sampling site by farm is reported in Table 2. A greater number of unique functions were identified in fecal samples from farm B (P = 0.05). Still, no difference in functional annotation was detected in reads from rumen fluid or DNS samples between farms.
Normalized read counts for the 10 highest ranked E.C.s for each nominal class were stratified by farm and sampling site and visualized on a heatmap (Figure 6). The normalized read counts for the 10 E.C.s most informative in predicting farm are given in Figure 6A, and those E.C.s were mapped to a number of KEGG pathways and biological functions, including microbial metabolism in diverse environments, lipopolysaccharide synthesis, and T cell receptor signaling ( Table 3). The normalized read counts for the 10 E.C.s most informative in distinguishing body site are visualized in Figure 6B. Those E.C.s were mapped to a number of KEGG pathways and biological functions, including microbial metabolism in diverse environments, methane metabolism, and plant cell wall degradation ( Table 4). The normalized read counts for the 10 E.C.s most informative in differentiating the composite farm by site nominal class are illustrated in Figure 6C, and those ECs were mapped to numerous KEGG pathways and biological functions, including microbial metabolism in diverse environments, amino acid metabolism, DNA repair, and gene expression checkpoints, as well as plant cell wall degradation (Table 5).

DISCUSSION
To better understand the consequences of long-term exposure to contaminated drinking water in the metagenome of dairy cattle, we sequenced the metagenome of fecal, rumen fluid, and deep nasopharyngeal samples, and characterized the differences between animals in farms from heavy metal contaminated and non-contaminated areas. Our data revealed differences in bacterial taxa between cows from both areas, and the abundance and prevalence of AMR genes was higher in the affected farm (farm A) than the unaffected farm (farm B). This was the first study to assess the effect of an environmental in the metagenome of animals exposed to drinking water contaminated with high levels of heavy metals.

Microbiome
Long-term exposure to heavy metals seemed to interfere in the fecal microbiota at phylum-level, since Firmicutes is described as the most abundant phylum in fecal samples of dairy cows in studies unrelated to heavy metal contamination (Liu et al., 2016;Xu et al., 2017;Hagey et al., 2019) and in the present research, Firmicutes was more abundant in farm B. Rodent models showed differences in gut microbiota after heavy metal exposure, possibly due to the direct effect of non-absorbed elements (Breton et al., 2013). In addition, the long-term effect of heavy metal pollution in the gut microbiota of Bufo raddei (Mongolian toad) revealed decreases in species diversity and in population of probiotic bacteria . Our findings suggest that long-term exposure to heavy metal contaminated drink water may also interfere in the microbiota of dairy cattle feces at genus-level, probably due to the presence of more sensitive and more tolerant groups. Analysis of rumen fluid at phylum-level revealed Bacteroidetes and Firmicutes as the major phyla. The rumen fluid microbiota is mainly composed of gram-negative bacteria in animals fed with high forage diets (Hungate, 1966). In our present research, dairy cattle from both herds had a forage-based diet and our results support this finding. Fibrobacter succinogenes (genus Fibrobacter) are butyrate-producing bacteria present in the rumen fluid, and this genus was higher in farm A compared to farm B. Prevotella and Bacteroides were the most abundant genera in almost all samples of this study. Prevotella has been reported to be the most abundant genus in the rumen (Chaucheyras-Durand and Ossa, 2014;Jami et al., 2014;Lima et al., 2015;Tong et al., 2018), whose species are related to degradation of starch, hemicellulose, and protein (Henderson et al., 2015). Prevotella has also been positively correlated with milk production (Lima et al., 2015;Indugu et al., 2017). Here, Ruminococcus and Eubacterium were numerically higher in farm A and farm B, respectively. The study conducted by Tong et al. (2018) showed higher abundances of Ruminococcus 2 and Eubacterium coprostanoligenes in high yield cows. Eubacterium has also been associated with low residual feed intake (Elolimy et al., 2018), and this genus is involved in cellulose degradation in the rumen and some species    are involved in the synthesis of butyrate (Flint et al., 2007), an energy source for cattle. We hypothesize that the lower abundance of Eubacterium and Prevotella in cows of farm A may be contributing to lower reported milk production for this farm. In fact, cows in farm A produced 250 L milk/per day before the dam failure, and, recently, the milk production is 200 L milk/per day. Research on nasopharyngeal samples is often conducted to understand the bovine respiratory disease or changes in the microbiota in feedlot cattle. For the first time, we compared the deep nasopharyngeal microbiota of dairy cows from a potential heavy metal contaminated environment and a non-contaminated environment using animals with no signs of respiratory disease. Analysis at phylum-level revealed Proteobacteria, Firmicutes, Actinobacteria, and Bacteroidetes as the most abundant groups. Similar results were described by Lima et al. (2016) and Gaeta et al. (2016), in which studies Proteobacteria was also the most abundant phylum in healthy and pneumonic groups. At the genuslevel, Eubacterium, Streptococcus, Acinetobacter, and Bacillus were the major groups observed in this study. These findings are in accordance with what has been previously reported, with Acinetobacter being commonly detected in the nasopharyngeal microbiota of cattle (Gaeta et al., 2016;Holman et al., 2017;Zeineldin et al., 2017). Streptococcus is an important genus related to pneumonia in humans, particularly Streptococcus pneumoniae (Jain et al., 2015), and has been detected in the respiratory tract of cattle (Klima et al., 2019). Bacillus are ubiquitously present in the environment and could be inhaled by cattle and detected in the upper respiratory tract. Some Bacillus species have the ability to produce antimicrobial peptides (Shelburne et al., 2006), that have been proven to inhibit the in vitro growth of Pasteurella multocida, Mannheimia haemolytica, and Histophilus somni (Xie et al., 2009), important pathogens involved in bovine respiratory disease complex. The presence of Bacillus as one of the major genus in the nasopharyngeal samples may be one reason for the lower abundance of Pasteurella, Mannheimia, and Histophilus observed in animals enrolled in this study, and, consequently, might have contributed to the absence of signs of respiratory disease.

Resistome
Long-term persistence and continuous deposition of heavy metals in the environment (e.g., water source and soil) enable microbial interactions and transference of resistance genes, potentiating the emergence of AMR bacteria in food animals and increase the risk of human exposure through the food chain (primary production, food industry, and household) (Oniciuc et al., 2019;Pérez-Rodríguez and Mercanoglu Taban, 2019). Efflux pumps are protein systems related to cross-resistance (Baker-Austin et al., 2006;Martinez et al., 2009), as several antimicrobial molecules (e.g., heavy metals and antibiotics) are pumped out of the cell by the same system. Clinically important systems are members of RND superfamily in gram-negative bacteria (e.g., Escherichia coli, Campylobacter jejuni, and Pseudomonas aeruginosa), MFS and ABC in gram-positive bacteria and mycobacteria (e.g., Staphylococcus aureus, Bacillus spp., and Mycobacterium tuberculosis) (Piddock, 2006;Routh et al., 2011). Our results suggest that both cross-resistance and co-resistance mechanisms may be occurring more often in farm A due to the higher prevalence of efflux pumps (including multi-metal resistance RND and multibiocide resistance ABC), in addition to protein types conferring resistance to several antibiotic classes, possibly enhanced by environmental contamination. Furthermore, the detection of several resistance mechanisms warn for possible presence of multidrug resistant bacteria in feces and rumen that may contaminate raw milk (Tóth et al., 2020), carcasses, abattoir surfaces, workers, and consumers, which have been extensively documented in the literature (Schlegelová et al., 2004;Lavilla Lerma et al., 2014;Savin et al., 2020). AMR development is a natural event that can be exacerbated by therapy administration and/or growth promoters (Pérez-Rodríguez and Mercanoglu Taban, 2019). In the present research, we evaluate the long-term exposure to heavymetal contaminated drinking water in dairy cows as a potential new pathway to AMR development in farm animals. Our data showed that fecal samples exhibited the highest prevalence of AMR genes compared to rumen fluid and DNS, which implicates future sampling designs using shotgun sequencing.

Functional Annotation
The symbiotic relationships among rumen fluid microbial communities are responsible for producing microbial protein by ammonia nitrogen utilization. The ammonia assimilation occurs by two main pathways in which the participation of glutamate synthase is fundamental. This enzyme is classified according to the electron donor; NADH-dependent, NADPHdependent, which is commonly found in several bacteria, and the ferredoxin-dependent glutamate synthase, which is found in plants and cyanobacteria (Pengpeng and Tan, 2013). Here, the results showed that ferredoxin-dependent glutamate synthase (E.C. 1.4.7.1) was overrepresented in rumen fluid samples of both farms, and in DNS of farm B. Cyanobacteria were overrepresented in two DNS samples of farm B, which may explain the results for this site. Glutamine-dependent asparagine synthase (E.C. 6.3.5.4), also identified in ammonia assimilation in rumen bacteria (Pengpeng and Tan, 2013), was overrepresented on farm B ruminal samples. This in addition to the presence of ferredoxindependent glutamate synthase suggests that the production of microbial protein is higher in the rumen environment of farm B dairy cows.
The 2,3-diphosphoglycerate dependent phosphoglycerate mutase (E.C. 5.4.2.11) was overrepresented in rumen fluid compared to feces but highly overrepresented in DNS, and probably independent of farm. A metal independent enzyme is present in vertebrates, yeast, and several bacteria (e.g., Haemophilus influenzae). This enzyme participates in glycolysis and gluconeogenesis processes and methane metabolism (Jedrzejas, 2000). In addition, DNA ligase (E.C. 6.5.1.2) is also represented in feces and rumen fluid, but it is overrepresented in DNS samples of farm B. This enzyme, commonly found in bacteria, participates in the DNA replication by forming phosphoric-ester bonds and repairs DNA breaks. The NADP-dependent glucose 6-phosphate dehydrogenase (E.C. 1.1.1.49), overrepresented in farm B, participates in the pentose phosphate pathway, which mainly produces NADPH and ribose-5-phosphate, a precursor of nucleic acids (Stanton, 2012). Finally, the DNA-directed DNA polymerase (E.C. 2.7.7.7), also overrepresented in farm B, is responsible for the synthesis of a new strand of DNA by adding nucleotides (Roettger et al., 2010), and it is also implicated in repair of DNA breaks. Taken together, these findings lead us to hypothesize that cows of farm B have more active bacterial replication compared to farm A, and bacteria in farm A are more susceptible to mutation and death, as they might have decreased ability to maintain DNA integrity. Thus, more research is required to understand differences in bacterial metabolism and behavior in cows of heavy metal contaminated and noncontaminated farms. It is worth mentioning that DNS samples were characterized by high variation in functional annotation, which may be due to low quantity DNA and the difficulty of DNA extraction from swabs. To better understand what is occurring in the nasopharyngeal environment, sampling and extraction protocols need to be optimized with more representative samples.
This broad, exploratory study is limited in the ability to draw mechanistic conclusions. One notable limitation is the lack of baseline pre-environmental disaster data to identify if the metagenome and resistome differed between the two farms preexposure. Because the AMR gene count data was normalized and scaled to account for gene length and sequencing depth bias, count-based statistical analyses (i.e., alpha and beta diversity, differential abundance analysis) were not possible. Also, diet can influence microbiota of dairy cattle (Loor et al., 2016;Dill-McFarland et al., 2018) but we did not directly measure this effect. In addition, cows within each farm were given the same diet, and the confounding effects of diet are already accounted for by a farm comparison. All samples were processed at once in the same batch, eliminating concerns regarding a true batch effect in terms of contamination or sequencing errors. Cows were the same breed, though we do not have information regarding their relatedness.

CONCLUSION
Long-term persistence of heavy metals in the environment may interfere in the microbiome of dairy cows. These data suggest that exposure to heavy metal contamination results in selection for bacteria that confer resistance to biocides, drugs, and metals and that AMR genes are most readily detected in fecal samples compared to the rumen fluid and nasopharyngeal samples. In addition, metagenome functional annotation data suggested that selective pressures of heavy metal exposure could be skewing pathway diversity toward fewer, more specialized, functions. Differences in microbiota indicated that exposure to heavy metal contaminated water may interfere in dairy cattle microbiota and, consequently, in animal productivity. Further research is warranted to determine if AMR is transferred via the food chain by milk or meat consumption, which would have substantial implications for human health. Since heavy metal contamination has an effect on animal microbiomes, environmental management, such as correct dam construction and maintenance and AMR surveys are warranted to protect the food system from hazardous consequences.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The name of the repository/repositories can be found below: https://www.mg-rast.org/linkin.cgi?project= mgp92419, https://github.com/gandalab/amr-heavymetal.

ETHICS STATEMENT
The animal study was reviewed and approved by Ethics Committee in the Use of Animals of the School of Veterinary Medicine and Animal Science of the University of São Paulo (protocol no. 6424070217). Written informed consent was obtained from the owners for the participation of their animals in this study.

AUTHOR CONTRIBUTIONS
NG, LG, and EG contributed to conception and design of the study. NG, DC, JC, and MA collected the data. NG, EB, and AM organized the database. NG, EB, and AM performed the statistical analysis. NG wrote the first draft of the manuscript. EB and AM wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

ACKNOWLEDGMENTS
Authors are thankful to Yuka Imamura Kawasawa and the Genome Sciences and Bioinformatics Core at the Penn State College of Medicine for performing the sequencing protocol, farmers for access, Istvan Albert and the Penn State Bioinformatics Consulting Center for guidance on metagenomic data analysis and Virology Laboratory at the Department of Microbiology, Institute of Biomedical Sciences, University of São Paulo for processing samples (extraction step).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2020.590325/full#supplementary-material Supplementary Figure S3 | AMR gene distribution across samples is highly heterogeneous and over 50% of AMR genes are present in fewer than 10% of samples. Antimicrobial resistance gene presence is heterogeneous across samples. The number of genes in each range is shown at the top of the bar. Over 50% of AMR genes (317 of 549) are present in fewer than 10% of samples. One gene (PH23S for Phenicol resistance) is present in all samples. The one gene present in all samples: Drugs| Phenicol| Phenicol-resistant_23S_rRNA_mutation| PH23S| RequiresSNPConfirmation.
Supplementary Figure S4 | Drug resistance prevalence and sample heterogeneity across farms. Drug resistance gene presence across farms. The top (blue) bar plot and dot matrix are the number of samples that share common gene presence. The left (gray) histogram is the frequency of each drug resistance type across samples in each.
Supplementary Table S1 | DNS genera that had the most differential abundance (largest log2 fold change) between farm A and farm B. The 25 most increased genera in farm A and the 25 most increased genera in farm B are shown. P values were adjusted with False Discovery Rate.
Supplementary Table S2 | Ruminal fluid genera that had the most differential abundance (largest log2 fold change) between farm A and farm B. The 25 most increased genera in farm A and the 25 most increased genera in farm B are shown. P values were adjusted with False Discovery Rate.
Supplementary Table S3 | DNS genera that had the most differential abundance (largest log2 fold change) between farm A and farm B. The 25 most increased genera in farm A and the 25 most increased genera in farm B are shown. P values were adjusted with False Discovery Rate.