Insights Into the Resistome of Bovine Clinical Mastitis Microbiome, a Key Factor in Disease Complication

Bovine clinical mastitis (CM) is one of the most prevalent diseases caused by a wide range of resident microbes. The emergence of antimicrobial resistance in CM bacteria is well-known, however, the genomic resistance composition (the resistome) at the microbiome-level is not well characterized. In this study, we applied whole metagenome sequencing (WMS) to characterize the resistome of the CM microbiome, focusing on antibiotics and metals resistance, biofilm formation (BF), and quorum sensing (QS) along with in vitro resistance assays of six selected pathogens isolated from the same CM samples. The WMS generated an average of 21.13 million reads (post-processing) from 25 CM samples that mapped to 519 bacterial strains, of which 30.06% were previously unreported. We found a significant (P = 0.001) association between the resistomes and microbiome composition with no association with cattle breed, despite significant differences in microbiome diversity among breeds. The in vitro investigation determined that 76.2% of six selected pathogens considered “biofilm formers” actually formed biofilms and were also highly resistant to tetracycline, doxycycline, nalidixic acid, ampicillin, and chloramphenicol and remained sensitive to metals (Cr, Co, Ni, Cu, Zn) at varying concentrations. We also found bacterial flagellar movement and chemotaxis, regulation and cell signaling, and oxidative stress to be significantly associated with the pathophysiology of CM. Thus, identifying CM microbiomes, and analyzing their resistomes and genomic potentials will help improve the optimization of therapeutic schemes involving antibiotics and/or metals usage in the prevention and control of bovine CM.


INTRODUCTION
Mastitis is the foremost production and major economic burden confronted by the global dairy industry (Reyes-Jara et al., 2016;Hoque et al., 2019). Bovine clinical mastitis (CM) milk is now considered to host a complex microbial community with great diversity (Oikonomou et al., 2014;Falentin et al., 2016;Hoque et al., 2019). The most frequently isolated pathogens are Staphylococcus aureus, Escherichia coli, Klebsiella spp., Streptococcus spp., Mycoplasma spp., Enterobacter spp., Bacillus spp., and Corynebacterium species (Abebe et al., 2016;Gao et al., 2017;Hoque et al., 2018). Therefore, the accurate identification of pathogens that cause CM enables appropriate choices for antimicrobial treatment and preventive mastitis management (Preethirani et al., 2015;Van Boeckel et al., 2015;Cheng et al., 2019). Over the past two decades, a wide range of phenotyping and genotyping methods have been implemented to study mastitis-causing bacteria (Preethirani et al., 2015;Gao et al., 2017;Hoque et al., 2018;Cheng et al., 2019). Although culture-based techniques are traditionally used to detect CM bacteria, these methods are time-consuming and have the inherent drawback of not being applicable to non-cultivable bacteria (Baron et al., 2018). Until recently, 16S rRNA partial gene sequencing remained the most commonly used genomic survey tool to study bovine mastitis microbiomes (Oikonomou et al., 2014;Falentin et al., 2016;Cremonesi et al., 2018). However, this technique has limitations because of polymerase chain reaction (PCR) bias, lower taxonomic resolution at the species level, and limiting information on gene abundance and functional profiling (Oniciuc et al., 2018). Shotgun whole metagenome sequencing (WMS), on the other hand, better characterizes the breadth of microbial diversity in a sample and successfully provides insight into the phylogenetic composition, species and/or strain, and functional diversity for a variety of biomes (Seth et al., 2014;Oniciuc et al., 2018). This WMS typically produces high complexity datasets with millions of short reads allowing extensive characterization of the microbiome in an ecological niche, profiling its functional attributes, and gradually becoming a cost-effective metagenomic approach (Seth et al., 2014;Oniciuc et al., 2018;Hoque et al., 2019).
Currently, antimicrobial treatment is indispensable to keeping bovine udder health, animal welfare, and economic aspects in balance (Preethirani et al., 2015;Krömker and Leimbach, 2017;Cheng et al., 2019). Therefore, dependence on antimicrobials has become a widespread phenomenon on dairy farms for mastitis management, prevention, and control programs. However, efficacy of antimicrobial therapy against bovine CM pathogens is low (Cheng et al., 2019), and the use of antibiotics is confined to selected severe CM cases only (Reyes-Jara et al., 2016;Cheng et al., 2019). The prevalence of antimicrobial resistance (AMR) in bovine CM pathogens has been investigated in numerous studies (Preethirani et al., 2015;Krömker and Leimbach, 2017;Cheng et al., 2019). The secretion of antimicrobial compounds by microbes is an ancient and effective method to improve the survival of microbes competing for space and nutrients with other microorganisms (D'Costa et al., 2011). The vast diversity of bacterial species in CM milk, many with short generation times and rampant horizontal gene transfer, permit the rapid accumulation of countless resistant variants at a relatively high evolutionary pace (D'Costa et al., 2011;Weller and Wu, 2015). However, resistance in CM bacteria typically goes unnoticed until a given species becomes of clinical interest, and the associated resistome is also suspected to be a source of newly emerging resistance genes in CM pathogens (Krömker and Leimbach, 2017;Cheng et al., 2019;Hoque et al., 2019;Zaheer et al., 2019). Bacteria residing in the bovine gastrointestinal tract and udder may become resistant to these antibiotics and, once released into the milk, they may take part in horizontal transfer of antibiotic resistance genes (ARGs) to other CM bacteria of contagious and environmental origin (Cheng et al., 2019;Zaheer et al., 2019). Furthermore, AMR is a global health concern in both human and veterinary medicine (Van Boeckel et al., 2015), and thus, monitoring the emergence of AMR strains is an essential component of bovine CM prevention and control strategies (Tomazi et al., 2018;Cheng et al., 2019). Therefore, finding an effective alternative strategy for the control of bovine mastitis is a challenge for dairy producers.
The antimicrobial properties of metals have been documented throughout the history of medicine and healthcare . Metal salts such as chromium (Cr), cobalt (Co), nickel (Ni), copper (Cu), and zinc (Zn) are effective in controlling bacterial transmission and infection risks . However, their uses are limited due to their toxicity and possible detrimental environmental effects in dairy industries particularly as therapeutic agents against bovine CM pathogens. Biofilm formation (BF) is an important virulence factor for mastitis causing bacteria and contributes to resistance to different classes of antimicrobials (Singh et al., 2017). Bacterial pathogens identified in this study showed a broad spectrum of antimicrobial (antibiotics, toxic metals) resistance, and possessed biofilm forming and quorum sensing (QS) abilities, which might be potential factors in hindering CM cures; thereby leading to the persistence of the disease, and the increased risk of transmission to non-infected dairy cows. However, genetic information about resistance or in vitro assays of resistance is not enough to understand the resistome when considered in isolation rather than in combination (Van Boeckel et al., 2015;Baron et al., 2018). Here we describe the resistance potentials in the CM microbiome of four major cattle breeds (Local Zebu, LZ; Red Chattogram Cattle, RCC; Sahiwal, SW; Crossbred Holstein Friesian; XHF) of Bangladesh using both metagenomic deep sequencing (WMS) and in vitro cultural approaches. We aim to investigate the influences of metabolic genomic potentials of the microbiomes in predicting and understanding the role of resistant potentials in the pathophysiology of bovine CM. We also test if cattle breeds or host genetics influence the milk microbiota composition and susceptibility to disease and resistance to bacterial infection Curone et al., 2018;Gonzalez-Recio et al., 2018).

Ethics Statement
The protocol for milk sample collection from lactating dairy cows was approved by the Animal Experimentation Ethical Review Committee (AEERC), Faculty of Biological Sciences, University of Dhaka under reference number 79/Bio. Scs., dated: 12-12-2019.

Screening for Clinical Mastitis and Sampling
We screened 260 udder quarter milk samples collected from 260 CM affected cows belonging to 50 smallholding dairy farms in two geographical regions of Bangladesh (central region, CR = 160; southeastern region, SER = 100) (Supplementary Figure S1). The cows represented four different breeds, including local zebu (LZ), red Chattogram cattle (RCC), Sahiwal (SW), and crossbred Holstein Friesian (XHF) at their early stage of lactation (within 10-40 days post-calving). A screening test for CM was conducted using the California Mastitis Test (CMT R , Original Schalm reagent, ThechniVet, United States) (Hoque et al., 2015). Approximately 15-20 mL of milk from each cow was collected under aseptic conditions in a sterile falcon tube during morning milking (8.00-10.00 am) and kept on ice (at 4 • C) for transport to the laboratory for subsequent processing.

Metagenomic DNA Extraction and Sequencing
Genomic DNA (gDNA) from 25 randomly selected CM samples was extracted by an automated Maxwell 16 DNA extraction platform using blood DNA purification kits (Promega, United Kingdom) following previously described protocols (Hoque et al., 2019). DNA quantity and purity were determined with NanoDrop (ThermoFisher, United States) by measuring 260/280 absorbance ratios. Sequencing libraries were prepared with the Nextera XT DNA Library Preparation Kit (Head et al., 2014) and paired-end (2 × 150 bp) sequencing was performed on a NextSeq 500 machine (Illumina Inc., United States) at the George Washington University Genomics Core facility. Our metagenomic DNA yielded a total of 596.74 million reads with an average of 23.87 million (maximum = 39.75 million, minimum = 8.89 million) reads per sample (Supplementary Material) before cleaning.

Microbiome Diversity and Community Analysis
The shotgun WMS data were analyzed using both mapping-based and assembly-based hybrid methods of PathoScope 2.0 (PS) (Hong et al., 2014) and MG-RAST (MR) (Glass et al., 2010), respectively. In PS analysis, a "target" genome library was constructed containing all bacterial sequences from the NCBI Database using the PathoLib module (Hong et al., 2014). The reads were then aligned against the target libraries using the very-sensitive Bowtie 2 algorithm (Langmead and Salzberg, 2012) and filtered to remove the reads aligned with the cattle genome (bosTau8) and human genome (hg38) as implemented in PathoMap (-very-sensitive-local -k 100-score-min L,20,1.0). Finally, the PathoID (Francis et al., 2013) module was applied to obtain accurate read counts for downstream analysis. In these samples, 17.20 million reads (4.3% of total reads) mapped to the target reference genome libraries after filtering for the cow and human genome (Supplementary Material). The raw sequences were simultaneously uploaded to the MR server (release 4.0), with proper embedded metadata, and were subjected to the quality filter containing dereplication and removal of host DNA by screening for taxonomic and functional assignment. Alpha diversity (diversity within samples) was estimated using the observed species, Chao1, ACE, Shannon, Simpson and Fisher diversity indices (Koh, 2018) for both PS and MR read assignments and counts. To visualize differences in bacterial diversity, a principal coordinate analysis (PCoA) was performed based on weighted-UniFrac distances (for PS data) through the Phyloseq R package, version 3.5.1 (McMurdie and Susan, 2013) and Bray-Curtis dissimilarity matrix (Beck et al., 2013) (for MR data). We also used OmicCircos, version 3.9 (Hu et al., 2014a), an R package based on python scripts, for circular visualization of both microbiome diversity and resistance to antibiotics and toxic compounds (RATC) functional groups found in MR data for our four targeted breeds of CM cows.

In vitro Identification of Bacteria
Collected CM milk samples (n = 260) were subjected to selective isolation and identification of S. aureus, E. coli, Klebsiella, Enterobacter, Shigella, and Bacillus species according to previously described microbiological methods (Reyes-Jara et al., 2016;Gao et al., 2017;Hoque et al., 2018;Cheng et al., 2019). The pathogens were identified based on their colony morphology, hemolytic patterns on blood agar and Gram-staining (Cheng et al., 2019). Gram-positive bacteria were further confirmed based on their biochemical characteristics in indole, methyl red, Voges-Proskauer (VP), catalase, oxidase, urease and triple sugar iron (TSI) tests, and growth on mannitol salt agar. Gram-negative bacteria were confirmed based on the results of indole, methyl red, citrate (IMViC) tests and lactose fermentation on MacConkey agar (Preethirani et al., 2015;Gomes et al., 2016). Finally, all isolates were stored at −80 • C for further genomic identification.

PCR Amplification and Ribosomal (16S rRNA) Gene Sequencing
Genomic DNA of probable S. aureus, E. coli, Klebsiella, Enterobacter, Shigella, and Bacillus species was extracted from overnight cultures using the boiled method Queipo-Ortuño et al., 2008). The quantity and purity of the extracted DNA was determined as mentioned before. The 16S rRNA gene was amplified using universal primers 27F (5 -AGAGTTTGATCCTGGCTCAG-3 ) and U1492R (5 -CTACGGCTACCTTGTTACGA-3 ) (Masomian et al., 2016). Agarose gel electrophoresis (1.2% wt/vol) was used to verify the presence of PCR products. DNA sequencing was carried out at First Base Laboratories Sdn Bhd (Malaysia) using Applied Biosystems highest capacity-based genetic analyzer (ABI PRISM R 377 DNA Sequencer) platforms with the BigDye R Terminator v3.1 cycle sequencing kit chemistry.

Phylogenetic Analysis of the Microbial Communities
Taxonomic abundance of the WMS data was determined by applying the "Best Hit Classification" option in the PS pipeline using the NCBI database as a reference, with the following settings: maximum e-value of 1 × 10 −30 , minimum identity of 95% for bacteria, and a minimum alignment length of 20 as the set parameters (Masomian et al., 2016). A midpoint rooted phylogenetic tree consisting of the top 200 abundant bacterial strains, identified through PS analysis from the WMS reads of the 25 CM samples with >90% taxonomic identity, was constructed using the maximum-likelihood method in Clustal W, version 2.1 (Larkin et al., 2007) and visualized using the interactive Tree Of Life (iTOL) (Letunic and Bork, 2011). Another phylogenetic tree, which was constructed using the same approach, focused on 40 strains corresponding to the six in vitro examined CM bacteria found in 260 CM samples with >90% taxonomic identity. Using Molecular Evolutionary Genetics Analysis (MEGA) version 7.0 for the larger datasets (Kumar et al., 2016), the 16S rRNA gene sequences, amplified from all individual bacterial isolates, were aligned with each other and with relevant reference sequences obtained from the NCBI Database using MEGA 7.0, and a maximum-likelihood tree was generated using these 16S rRNA gene sequences with the Tamura-Nei evolutionary model (Kumar et al., 2016). Nodal confidence in the resulting phylogenetic relationships was assessed using the bootstrap test (1000 replicates) (Pattengale et al., 2010).

Metal Susceptibility Testing
The antibacterial effect of heavy metals was evaluated in vitro for the isolated pathogens using the agar well diffusion method (Supplementary Figure S3) (Reyes-Jara et al., 2016;Vaidya et al., 2017). Five heavy metals such as copper (Cu), zinc (Zn), chromium (Cr), nickel (Ni), and cobalt (Co) were used as salts: CuSO 4 ·5H 2 O, ZnSO 4 ·7H 2 O, K 2 Cr 2 O 7 , NiCl 2 , and CoCl 2 ·6H 2 O, respectively, to study the level of zone of inhibition (ZOI). Briefly, pure culture of the isolated pathogens from NA plates were subcultured into Mueller-Hinton agar (Oxoid TM , United Kingdom) plates, and respective agar was poured into sterile Petri dishes, which were then cooled. A total of 100 µl of cell suspension was pipetted and spread across the entire area of the agar using a sterile cotton swab. Two to five equal wells (7 mm diameter) were cut out of each agar plate using a sterile cork borer and stainless-steel needle. Varying concentrations of the metal solutions were prepared (2,4,8,16,32,48, and 64 µg/mL) and 100 µl of the prepared metal ion solution was added to each of the wells. The plates were incubated at 37 • C for 24 h to allow diffusion of the metal into the agar, and the antibacterial activity was determined by measuring the diameter of ZOI in mm . After investigating the resistance profile of the isolates at different concentrations, the minimal inhibitory concentration (MIC) of the metals was determined using the tube dilution method, by gradually increasing or decreasing the heavy metal concentrations (Reyes-Jara et al., 2016). Finally, growth of bacterial colonies was observed and the concentration that showed no growth was considered as the minimum bactericidal concentration (MBC) (Reyes-Jara et al., 2016).

Biofilm Assay and Microscopy
Microtiter plate assays were performed to screen for the BF ability of 80 randomly selected isolates using standard protocols (Schönborn et al., 2017;Singh et al., 2017;Vaidya et al., 2017). We quantified the absorbance of solubilized crystal violet (CV), in a plate reader at 600 nm using 30% acetic acid in water as the blank and TSB as the negative control. The solution was removed, and the absorbance measured at optical density-590 (OD590) (n = 3). To determine the BF ability of strains, cut-off optical density (ODc) was defined as three standard deviations above the mean OD of the negative control. Strains were classified as: non-biofilm formers, NBF (OD ≤ ODc); weak biofilm formers, WBF (ODc < OD ≤ 2 × ODc); moderate biofilm formers, MBF (2 × ODc < OD ≤ 4 × ODc), and strong biofilm formers, SBF (OD > 4 × ODc) (Tiwari et al., 2017;Vaidya et al., 2017). In this study, the ODc value was set as 0.045 and the mean OD of the negative control was 0.039 ± 0.002 .
The biofilm surfaces were then visualized using 5% TSB as nutrient rich media and FilmTracer TM LIVE/DEAD R Biofilm Viability Kit as staining materials to observe the proportion of live or active cells (fluorescent green) under Olympus BX51 upright microscope (40× objective), and finally, images were collected using an Olympus DP73 camera through cellSens entry software (Olympus Corporation, Japan) and visualized using image J software (Schönborn et al., 2017). As a negative control, we used E. coli DH5 alpha for all the in vitro resistome (antimicrobial and metal susceptibility tests and biofilm assays) analysis tests.

Microbial Functional Analysis
Metagenomic functional composition was based on the gene families from different levels of the SEED module and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2019), using the MR pipeline (Glass et al., 2010). We observed significant differences (Kruskal-Wallis test, P = 0.001) in the relative abundance of genes coding for RATC and microbial functional genomic potentials in four cattle breeds.

Statistical Analyses
The characteristics of the cow breeds with CM were compared using a non-parametric Kruskal-Wallis test for quantitative variables . The Shapiro-Wilk test was used to check normality of the data, and the non-parametric Kruskal-Wallis rank sum test was used to evaluate differences in the relative abundance of bacterial taxa at the strain level according to breed groups Curone et al., 2018;Gonzalez-Recio et al., 2018). The statistical analyses for the MR data were initially performed by embedded calls to statistical tests in the pipeline and validated further using IBM SPSS (SPSS, Version 23.0, IBM Corp., NY, United States) using the abovementioned tests. For the functional abundance profiling, the statistical (Kruskal-Wallis test and Pearson correlation) tests were applied at different KEGG and SEED subsystem levels in MR pipeline generated data (Glass et al., 2010). To evaluate the significant relationships between identified bacterial species and the study region, we used the two-sample proportions test using SPSS. Results were considered statistically significant when P < 0.05 and highly significant when P < 0.01. Mean values were used to compare the antimicrobial efficacy results of the tested antibiotics and heavy metals at varying concentrations. Standard error means were calculated to analyze the distribution of the data from the mean value, and confidence intervals of 95% were calculated for the MIC and MBC tests results to plot error bars (Schönborn et al., 2017;Vaidya et al., 2017). We also performed Pearson correlation tests to assess the relationships between the taxonomic abundance of the pathogens and AMR, both for cultural and metagenomic data. A post hoc Bonferroni test was used to compare the biofilm OD600 mean values (Schönborn et al., 2017;Vaidya et al., 2017).

RESULTS
To decipher the resistance potentials in bovine CM microbiomes, we used both in silico (WMS, 16S rRNA gene sequencing) and in vitro (culture base) approaches. The present WMS investigation leads to the direct and comprehensive evaluation of resistance to antibiotics and toxic compounds (RATC), BF and QS genes in 25 CM milk samples. Furthermore, in vitro AMR profiling of six CM causing bacteria (S. aureus, E. coli, Klebsiella, Enterobacter, Bacillus, and Shigella) isolated from 260 CM milk samples was carried out using 12 commonly used antibiotics (ampicillin, doxycycline, tetracycline, nitrofurantoin, ciprofloxacin, nalidixic acid, cefoxitin, imipenem, chloramphenicol, gentamycin, erythromycin, and vancomycin), and five toxic metals (copper, zinc, chromium, nickel, and cobalt). Moreover, we also demonstrated some functional metabolic potentials of CM microbiomes found to be associated with mammary gland pathogenesis.

Sequence Analysis
The WMS of 25 CM milk samples generated approximately 600 million reads, ranging from 8.86 to 39.75 million per sample. An average of 21.13 million reads per sample (maximum = 36.89 million, minimum = 4.71 million) passed the quality control step (Supplementary Material). We analyzed the sequencing reads simultaneously using two bioinformatics pipelines, PathoScope 2.0 (PS) and MG-RAST (MR).

Microbiome Diversity and Composition in CM
We investigated the strain-level microbial community and relative abundances in 25 CM milk samples [previously published 14 samples (Hoque et al., 2019) and 11 new samples] through WMS. The reads generated from WMS were mapped to 391 genera and 519 strains of bacteria through MR and PS analyses, respectively (Supplementary Material). The rarefaction curves based on observed species richness reached a plateau after, on average, 23.87 million reads ( Figure 1A and Supplementary Material), suggesting that the depth of coverage for most samples was sufficient to capture the entire microbial diversity within each sample. We did not, however, find any significant differences in the alpha (observed species, Chao1, ACE, Shannon, Simpson and Fisher diversity estimates) and beta (based on Bray-Curtis dissimilarity matrix) diversities among the microbial communities across the 25 CM samples (Figures 1B,C, respectively). However, significant diversity (alpha and beta) differences were observed among the CM microbiome communities across the four cattle breeds (LZ, RCC, SW, XHF) regardless of the method (i.e., either PS or MR) used to tabulate microbial abundances (PS; P = 0.005, MR; P = 0.001, Kruskal-Wallis test). In addition, this breed's specific diversity difference remained evident in the microbial ecosystem of XHF cows associated CM milk samples (Figures 1D,E, respectively). The PCoA analysis also showed significant microbial disparity (P = 0.001) among the microbiome of four dairy breeds ( Figure 1E).
The predominant bacterial phyla we found associated with CM were Proteobacteria, Bacteroidetes, Firmicutes, Actinobacteria, and Fusobacteria (contributing to >95.0% of the total sequences, Kruskal-Wallis test, P = 0.001) in the x-axis) on species richness (y-axis) in CM milk samples. The rarefaction curves representing the number of species per sample indicated that the sequencing depth was sufficient enough to fully capture the microbial diversity as existed. (B) Alpha diversity measured using the observed species, Chao 1, ACE and Shannon diversity indices through PathoScope (PS) analysis. The observed species richness (P Observed = 0.511), Chao1 (P Chao1 = 0.081), ACE (P ACE = 0.121), Shannon (P Shannon = 0.401), Simpson (P Simpson = 0.011) and Fisher (P Fisher = 0.014) diversity analyses revealed that microbiome diversity did not vary among the CM samples. (C) Beta diversity (Principal coordinate analysis; PCoA) measured on the Bray-Curtis distance method using MG-RAST tool for CM causing microbial communities (genus-level) shows that most of the CM samples clustered together (black circle), indicating no significant diversity differences. (D) Alpha diversity measured using species richness (P Observed = 0.011), Chao1 (P Chao1 = 0.001), ACE (P ACE = 0.021), Shannon (P Shannon = 0.001), Simpson (P Simpson = 0.009) and Fisher (P Fisher = 0.023) diversity matrices on PS data showed significant differences (Kruskal-Wallis test, P = 0.002) in microbial diversity across the four cow breeds (Local Zebu cows, LZ; Red Chattogram cows, RCC; Sahiwal, SW; Holstein Friesian cross, XHF). (E) PCoA plot based on weighted-UniFrac distance method at strain-level microbiome signature of four breeds of cows reveals that the CM samples appear more distantly (red circles) indicating significant group differences (P = 0.001). These differences in the microbiome signature associated with CM across the four breeds could be explained by a large percentage of variation in the first (62.6%) and second (19.7%) axes.
Simultaneously, through in vitro cultural analysis, a total of 452 isolates that belonged to six bacterial (S. aureus, E. coli, Klebsiella, Enterobacter, Bacillus, and Shigella) species were identified in 260 CM samples (including 25 WMS CM samples) collected from central (CR = 160) and southeastern (SER = 100) regions of Bangladesh (Supplementary Figure S1). The overall prevalence of S. aureus, E. coli, Klebsiella, Enterobacter, Bacillus, and Shigella species were 23.5, 18.5, 19.2, 12.3, 9.2, and 17.3% in CM samples, respectively (Supplementary Table S1). We found significant differences in the prevalence of these species (P = 0.01) when analyzing the distribution of these pathogens according to the origin of the samples FIGURE 2 | The strain-level taxonomic profile of microbiota associated with bovine clinical mastitis (CM). Taxonomic dendrogram showing the top bacterial microbiome of bovine CM milk. Color ranges identify different strains within the tree. Taxonomic dendrogram in the midpoint rooted phylogenetic tree was generated with the top 200 abundant unique strains of bacteria in CM milk metagenome based on the maximum likelihood method in Clustal W and displayed with iTOL (interactive Tree Of Life). Each node represents a single strain shared among more than 50% of the samples at a relative abundance of >0.0006% of the total bacterial community. Strains and/or species are color-coded by different order of bacteria present in >80% of samples. The strains in the phylogenetic tree are also available in Supplementary Material. Figure S4). The culture-based findings of the current study identified S. aureus as the chief etiology of bovine CM in Bangladesh, while Shigella species remained the least frequently detected CM pathogen-which corroborates the results of WMS-based taxonomic identification (Supplementary Figure S5).

Resistome Composition of CM Associated Microbiome
For analyses of the resistome composition in CM microbiomes, the SEED module of the MR pipeline provided a comprehensive picture. Using SEED, 147,040 reads aligned to 30 resistance to antibiotics and toxic compounds (RATC), and 10 BF and quorum sensing (BF-QS) functional groups across the CM samples, with different abundances (Supplementary Material). The RATC genes were classified into two unique groups, 19 antibiotic resistance and 11 toxic metal resistance groups (Figure 4 and Supplementary Material). This WMS analysis showed a significant association (Pearson correlation, P = 0.001; Non-parametric Spearman's Correlation, P = 0.003) between the number of reads aligned to bacterial genomes and the number of reads mapped to RATC genes (Supplementary Material). Among the RATC functional groups, multidrug resistance to efflux pumps (MREP, 28.6%), CmeABC operon (8.9%), resistance to fluoroquinolones (RFL, 6.2%), mdtABCD cluster (5.5%), methicillin resistance in Staphylococci (MRS, 3.8%), BlaR1 regulatory family (BlaR1, 3.4%), MexE-MexF-OprN (2.4%), and beta-lactamase resistance (BLAC, 2.2%) were the dominating antibiotic resistance genes (ARGs) found in CM milk microbiomes (Figure 4A and Supplementary Material). In addition to ARGs, the WMS analysis also detected a number of metal and toxic metals resistant genes in CM microbiomes. Among them, cobalt-zinc-cadmium resistance (CZCR, 19.3%), copper homeostasis (CH, 9.6%), arsenic resistance (AR, 2.9%), copper homeostasis: copper tolerance (CHCT, 2.3%), and resistance to chromium compounds (RCHC, 1.4%) were the predominate resistant genes ( Figure 4A and Supplementary Material). Although the relative abundance of these RATC genes varied among the microbiomes of the four breeds (LZ, RCC, SW, and XHF), their resistome composition did not vary significantly (P = 0.692) by taxonomic diversity of respective breeds (Figure 4B and Supplementary Material).
The in vitro antibiogram profiling of 221 individual isolates of the six bacteria, revealed that S. aureus isolates had the highest resistance to doxycycline, ampicillin, tetracycline, and erythromycin (73.0-88.0%) and moderate resistance to chloramphenicol, ciprofloxacin, and nitrofurantoin (50.0-58.0%) (Figure 6 and Table 1). The isolates of another Gram-positive bacterium (Bacillus) demonstrated the highest resistance against doxycycline, ampicillin, nalidixic acid, and erythromycin (60.0-84.0%). However, E. coli isolates exhibited the highest resistance against tetracycline, doxycycline, nalidixic acid, and ampicillin (77.0-93.0%) and moderate resistance to chloramphenicol, nitrofurantoin, gentamicin, and ciprofloxacin (40.0-63.0%). The isolates of Klebsiella, Enterobacter, and Shigella displayed the highest resistance to doxycycline, nalidixic acid, tetracycline, and ampicillin (70.0-100.0%) and moderate resistance to ciprofloxacin, gentamicin, nitrofurantoin, and chloramphenicol (30.0-70.0%). In this study, imipenem and cefoxitin remained as the most sensitive antibiotics against four Gram-negative bacterial (E. coli, Klebsiella, Enterobacter, and Shigella) species, while the two Gram-positive (S. aureus and Bacillus) species were mostly sensitive to imipenem, cefoxitin, and vancomycin (Figure 6 and Table 1). Taken together, the antibiogram profile revealed that all of the selected CM pathogens are becoming multidrug resistant (MDR, resistant to ≥3 antibiotics) and the highest resistance was found against tetracyclines (tetracycline and doxycycline), followed by quinolones (nalidixic acid), and penicillin (ampicillin) groups of antibiotics (Figure 6 and Table 1).

Microbiome Functional Analysis
We also investigated the possible links between chemotaxis and pathogenicity through the identification of putative genes or proteins associated with both flagellar motility and bacterial chemotaxis. The KEGG pathway analysis identified 48 protein families associated with flagellar motility in prokaryotes, and among them, flagellar hook-length control protein, FliK ( Another regulatory and cell signaling gene, endoplasmic reticulum chaperon grp78 (BiP), was also found as the single most abundant (93.8%) gene in the proteolytic pathways of the CM associated bacterial strains (Supplementary Figure S7 and Supplementary Material). A deeper look at microbial genes associated with phages-prophages, transposable elements, and plasmids revealed that pathogenicity islands related proteins such as methionine-ABC transporter substrate-binding protein (33.8%), GMP synthase (27.7%), tmRNA-binding protein; SmpB (16.0%), heat shock protein 60; GroEL (16.0%), and SSU ribosomal protein; S18p (6.1%) were predominantly abundant among the CM pathogens (Supplementary Material). We also found significant associations between the number of reads assigned to genes coding for AMR and biofilm-formationand-quorum-sensing (Pearson correlation, P = 0.0001), and in the relative abundance of genes coding for biofilm-formation and flagellar activities (Pearson correlation, P = 0.004). The SEED module analysis also enabled us to identify 28 different protein functions associated with oxidative stress responses that were mostly represented by catalase related proteins (26.7%), Cu-Zn-Fe-Mn mediated superoxide dismutase (12.7%), H 2 O 2inducible genes activator (7.8%), and paraquat-inducible protein B (7.3%) (Figure 9 and Supplementary Material).

DISCUSSION
Previously, we reported that the bovine CM milk microbiome harbors genes that may promote microbial pathogenicity, including AMR compared to the healthy state (Hoque et al., 2019). In this study, we employed a combination of both in silico (WMS) and in vitro (culture-based) approaches to elucidate the resistance potentials in CM associated microbiomes. Currently, it is feasible to use metagenomic sequencing as a diagnostic tool for the detection of microbiome and its associated resistome in environmental samples (Pérez-Cobas et al., 2013;Schaik, 2015;Lanza et al., 2018). Shotgun metagenomics (WMS) investigations are progressively being used to analyze the ensemble of genes that may encode antibiotic resistance in various microbial ecosystems, which are defined as the resistome (Forsberg et al., 2012;Schaik, 2015;Su et al., 2017;Lanza et al., 2018;Zaheer et al., 2018). The increasing prevalence of AMR in bacteria in bovine CM milk is one of the most important challenges that public health faces. It has been proposed that livestock production systems may contribute to an increased prevalence of antimicrobial resistant bacteria, and/or associated genes in the environment and therefore pose a risk to human health. Globally, more than FIGURE 9 | Projection of the clinical mastitis (CM) milk metagenome onto KEGG pathways. The whole metagenome sequencing (WMS) reveals significant differences (Kruskal-Wallis test, P = 0.001) in functional microbial pathways. A total of 28 genes associated with oxidative stress were found in CM microbiomes. Black lines with green circles delineate the distribution of the stress related genes according to their class across the CM metagenome. The diameter of the circles indicates the relative abundance of the respective genes. More details about these genes can be found in the text and Supplementary Material. 57 million kilograms of antibiotics are used annually in food animal production (Van Boeckel et al., 2015), which may select for antibiotic-resistant bacteria that persist throughout the meat and milk production chain. Investigation of the microbiome, and/or associated resistome of dairy animals especially in milk from cows with mastitis, and their environment may provide valuable data and models to estimate the public health risk of antibiotic-resistant human infections associated with antibiotic use in dairy animals (Zaheer et al., 2019). Our present findings are sufficiently enriched in taxonomic resolution and predicted protein functions, and corroborates the findings of several previous studies (Oikonomou et al., 2014;Falentin et al., 2016;Patel et al., 2017;Hoque et al., 2019). The occurrence of bovine mastitis could be affected by cattle breeds Curone et al., 2018;Gonzalez-Recio et al., 2018), and the diversity of CM-causing pathogens is associated with a broad range of host-defense mechanisms as part of its immunological arsenal (Thompson-Crispi et al., 2014;Li et al., 2019). We found significant differences in taxonomic diversity and abundances among the CM microbiomes of four dairy breeds. The XHF cows suffering from CM had higher microbial diversity at strain-level, and a significant proportion of the microbiota found to be shared with that of the other three breeds (LZ, SW, and RCC). Consistent with the results of earlier studies Curone et al., 2018;Gonzalez-Recio et al., 2018;Li et al., 2019), the taxonomic profile of the CM microbiomes found in four breeds of cows were dominated by phyla Proteobacteria, Bacteroidetes, Firmicutes, Actinobacteria, and Fusobacteria. This breed specific variation in taxonomic richness and diversity of the microbiome, especially in XHF and LZ cows, could be associated with their increased disease resistance or immune response Curone et al., 2018;Gonzalez-Recio et al., 2018), and rumen microbial features (e.g., taxa, diversity indices, functional categories, and genes) (Li et al., 2019). However, this study demonstrated that the resistome is not significantly correlated with different dairy breeds suggesting that resistance potentials of CM microbiome might be common among dairy breeds. In our study, lack of congruence between resistome and microbial communities in four dairy breeds suggested that other factors, in addition to antimicrobial use, might be associated with changes in resistome, and microbiome composition. It may also reflect the possibility that a small fraction of the total bacterial population was resistant to antibiotics (MacLean and Vogwill, 2015;Rovira et al., 2019), and that horizontal gene transfer disrupted the link between microbiome and resistome composition (Johnson et al., 2015). However, further investigations are necessary to evaluate the real effect of breed specific bacteria on cow mammary gland diseases. Recent understandings regarding evolutionary relationships of major CM causing bacteria are primarily based on 16S rRNA gene phylogenetic identification along with a few individual gene or protein sequences (Naushad et al., 2016), which often produces conflicting phylogenies. This study explored the possibility that the prevalence of CM milk pathogens could vary according to geographical locations, and farming (semi-intensive to intensive grazing system in SER, semi-intensive to free-range grazing systems in CR) systems (Reyes-Jara et al., 2016). These differences may imply that the etiology of bovine CM in Bangladesh could be related to the breed/host genetic factors Curone et al., 2018;Gonzalez-Recio et al., 2018;Li et al., 2019), types of feeding, and farm locations and types (Reyes-Jara et al., 2016), types of antibiotics and/or metals used for treatment resistomes, and other predisposing factors, as described in other countries (Preethirani et al., 2015;Reyes-Jara et al., 2016;Cheng et al., 2019).
Data presented here, coupled with the data reported in our earlier study (Hoque et al., 2019) provide important insights into the resistance potentials in CM microbiomes. Our results are concordant with MDR bacteria reported elsewhere from the milk of clinically infected cows Tomazi et al., 2018;Cheng et al., 2019), buffalo cows (Preethirani et al., 2015), and humans (Penders et al., 2013;Patel et al., 2017;Baron et al., 2018). Our findings linked multidrug resistance to efflux pumps (MREP), CmeABC operon, mdtABCD cluster, BlaR1 family, methicillin resistance in Staphylococcus (MRS), resistance to fluoroquinolones (RFL), and multiple metals resistance to CZCR and AR as the predominantly abundant antibiotics and toxic compounds resistance (RATC) functional groups in CM microbiomes, suggesting that bovine CM milk microbiome constitutes a good reservoir for AMR (Kumar and Varela, 2012;Penders et al., 2013;Baron et al., 2018;Hoque et al., 2019). It has been reported that efflux pumps regulated by two-component systems in several pathogens, including A. baumannii and K. pneumonia, provide multidrug resistance, which may limit the treatment options against bacterial infections of the mammary glands (Li and Nikaido, 2009;Tiwari et al., 2017). Relative over-expression of efflux pumps enhances resistance to antimicrobials by reducing the accumulation of antibiotics inside the bacterial cells, providing sufficient time for the bacteria to adapt to the antibiotics (slow phase antibiotic efflux), and through mutations or alteration of antibiotic targets (Yao et al., 2016;Tiwari et al., 2017). The CmeABC operon is highly potent against multiple antibiotics, promotes the emergence of ARGs, and confers exceedingly high-level resistance to fluoroquinolones (Singh et al., 2017). Therefore, multidrug resistance to efflux pumps and multiple heavy metals resistance represented ubiquitous resistance mechanisms among CM microbiomes, which might be associated with unethical overuse of antibiotics in dairy animals (Preethirani et al., 2015;Curone et al., 2018;Tomazi et al., 2018;Cheng et al., 2019;Zaheer et al., 2019), and extensive application of toxic chemicals and metals in agricultural use (Reyes-Jara et al., 2016;Vaidya et al., 2017), or might have a function in the gut microbiome that is still unknown (Penders et al., 2013;Hu et al., 2014b;Ciesinski et al., 2018). The RATC genes detected in this study are of particular interest because there is concern that the use of this class of antibiotics or metals in veterinary medicine, particularly for food animals, may contribute to the development of resistance to this class of antimicrobial options in humans (Penders et al., 2013;Hu et al., 2014b).
The in vitro antibiogram assays revealed that a high proportion of multidrug resistant (MDR) bacteria was frequently observed, where tetracyclines (tetracycline and doxycycline), quinolones (nalidixic acid), penicillins (ampicillin), and phenols (chloramphenicol) were the most common resistant antibiotics against the six selected CM pathogens. This finding of high MDR patterns for CM pathogens is in line with many previous studies on bovine mastitis (Preethirani et al., 2015;Cheng et al., 2019). The AMR profile of bovine CM pathogens for different antimicrobials could vary according to the type and origin of bacteria (Preethirani et al., 2015;Van Boeckel et al., 2015;Cheng et al., 2019) and host-population such as bovine (Tomazi et al., 2018;Cheng et al., 2019) and bubaline (Preethirani et al., 2015) cows. Consistent with bacterial needs, heavy metals can be transformed (e.g., oxidized, reduced, methylated, or complexed) and used as a source of energy, terminal electron acceptors, and/or enzyme structural elements (Drewniak et al., 2016). The highest abundance of CZCR genes among CM pathogens is mainly due to the presence of Co, Zn, and Cd detoxification systems (Drewniak et al., 2016). Although understanding of the uncontrolled spread of ARGs in bovine mastitis pathogens (Cheng et al., 2019) is growing, information on toxic compounds or heavy metal resistance remains unavailable. In this study, heavy metals (Cr, Co, Ni, and Cu) tested for antibacterial sensitivity showed good efficacy, but knowledge on their mode of action is limited. Thus, with the increase of MDR bacteria in CM, there is an imperative need for new biocidal and antimicrobial formulations. The MIC and MBC tested metals revealed effective antimicrobial efficacies against a wide range of AMR pathogens (Reyes-Jara et al., 2016;Vaidya et al., 2017;Ciesinski et al., 2018). We found that Cr and Co compounds had the highest antimicrobial efficacy (MIC) against all of the tested bacteria as supported by several previous studies Murcia et al., 2018). This finding is likely related to the better lipophilicity of Cr and Co, resulting in increased antimicrobial activities or chelating effects (Murcia et al., 2018). However, both Cr and Co, should be used cautiously either in medication or in feed supplementation after calculating appropriate concentrations, to avoid cytotoxicity. These in vitro resistance patterns of CM bacteria corroborate the resistome profile of the WMS, since multidrug resistance to efflux pumps, resistance to fluoroquinolones, beta lactamases, and methicillin were the predominant antibiotics resistance functional groups, whereas cobalt, zinc, cadmium, copper, arsenic, and chromium resistance were the predominating toxic metals resistant genes in CM causing microbiomes. Molecular methods complemented with culture-based diagnostics have been historically implemented to document AMR in bacteria, however, the advent of shotgun metagenomics (WMS) has revolutionized the method of addressing relevant problems like diagnosis and surveillance of infectious diseases, and the issue of AMR (De, 2019). Biofilm formation is an important virulence factor that may result in recurrent or persistent udder infections (Melchior et al., 2006), and treatment failure through increased resistance to antibiotics and protection against host defenses (Schönborn et al., 2017). The relative overexpression of genes encoding lsrACDBFGE operon, biofilm adhesion biosynthesis, protein YjgK cluster, and QS: autoinducer-2 synthesis in CM microbiomes is in accordance with several earlier reports (Gomes et al., 2016;Schönborn et al., 2017;Hoque et al., 2019). In this study, the relative abundance of the predicted proteins for BF and QS varied significantly among the six selected bacterial taxa. Bacterial BF and QS abilities can be the strain specific or genetically linked traits, representing a selective advantage in pathogenesis of bovine CM (Schönborn et al., 2017). Biofilms can enhance proliferation of reactive oxygen and nitrogen species that can survive antibiotic treatment leading to the transfer of ARGs (Gomes et al., 2016). In this study, overall, 76.2% of the isolates of the six selected CM pathogens were found as biofilm producers, and their ability to produce biofilm varied significantly (Gomes et al., 2016;Schönborn et al., 2017). A large number of food spoilage and/or pathogenic bacteria, including Enterococcus faecalis, Enterobacter spp., Pseudomonas spp., Klebsiella spp., S. aureus, E. coli, B. cereus, and others, have already been reported to be associated with biofilm-formation (BF) from diary niches (Melchior et al., 2006;Gomes et al., 2016;Schönborn et al., 2017;Singh et al., 2017;Vaidya et al., 2017) which supports our current findings.
Bacterial chemotaxis mediated by flagellar activities (Duan et al., 2013) and the flagella mediated virulence factors are found in many pathogenic species of bovine CM microbiomes, making them a potential target for new antibacterial therapeutics (Duan et al., 2013). The intra-and interspecies cell-to-cell communication in bovine CM microbiomes were associated with 26 different genes, which might have vital roles in the early phase of mastitis for attachment to, or entry into the udder tissues and virulence regulation (Zatakia et al., 2018), and bacterial colonization in mammary tissues like other suitable sites (Matilla and Krell, 2017). The cheA-cheY twocomponent system mediated bacterial chemotaxis also facilitates the initial contact of bacteria with mammary gland epithelial cells, and contribute to effective invasion (Dons et al., 2004). The two-component signal transduction system BarA-UvrY regulates metabolism, motility, BF, stress resistance, virulence, and QS in CM pathogens by activating the transcription of genes for regulatory small RNAs (Zere et al., 2015). The up-regulation of genes coding for proteolytic activity, grp78 during hostpathogen interactions in CM, is associated with endoplasmic reticulum (ER) stress which further triggers proteolytic activities to initiate the mechanism of pathogenesis and cell death (Hirai et al., 2018). Catalase activity is a marker of bovine mastitis, which plays a central role in milk redox control and markedly increases during the pathophysiology of bovine CM (Andrei et al., 2016). Our present findings are in line with previous reports (Andrei et al., 2016;Darbaz et al., 2019) that an elevated oxidative stress mediated by catalase activity might have originated either from the mammary gland and/or bacterial cells. During the pathogenesis of bovine mammary gland, bacteria are not rapidly killed by the phagocytic activity of bovine macrophages; rather, they survive within macrophages during prolonged infection due to secretion of catalase and superoxide dismutase, which by degrading H 2 O 2 , inhibit the ROS mediated killing mechanism of the host (Andrei et al., 2016;Darbaz et al., 2019). The majority of published studies on bovine mastitis describe 16S rRNA genebased community structure evaluations, whereas published reports on shotgun deep sequencing metagenomics of CM microbiome studies remain scarce. We completed an in-depth report simultaneously describing microbiome diversity along with previously unreported opportunistic strains, and their associated resistome in bovine mastitis.

CONCLUSION
Bovine CM milk is a potential reservoir of diverse groups of microbes harboring a diverse resistome and other virulence factors. Our findings reveal that genes coding for multidrug and multiple metal resistance are ubiquitously present in the CM microbiome, and protect bacteria from the antibacterial effects of antimicrobials, extruding them out of cells. The efflux pumps mediated by multidrug resistance and multiple metals (e.g., cobalt, zinc, cadmium, arsenic, chromium) resistance were the predominating genomic functional groups to be considered as the potential key factors for the persistence of bovine CM, and the failure of conventional therapies against CM-related pathogens. Additionally, BF, QS, bacterial flagellar movement and chemotaxis, regulation and cell signaling, and oxidative stress associated genes can be a great benefit to bacteria against the host's immune system. Microbial resistance and genomic functional potentials represent important mechanisms for antibiotic resistance, leading to treatment failure and persistence of the disease. An accurate and timely identification of CM-associated pathogens and analyses of their resistance potentials are necessary for the selection of proper therapeutics (antibiotics and/or metals), and the development of effective, safe, and economical treatment regimens for bovine mastitis and sustainable dairying. Although the baseline data presented here are promising, further studies are recommended using a larger sample size, and with the inclusion of gut microbiome sampling in addition to the milk samples to elucidate the specific resistance potential of CM pathogens as well as for direct testing of microbiome and resistome transfer across this axis.

DATA AVAILABILITY STATEMENT
The sequence data reported in this article have been deposited in the NCBI database (BioProject PRJNA529353 for metagenome sequences, NCBI accession numbers: MN 620423-MN 620430 for 16S rRNA gene sequences).

ETHICS STATEMENT
The protocol for milk sample collection from lactating dairy cows was approved by the Animal Experimentation Ethical Review Committee (AEERC), Faculty of Biological Sciences, University of Dhaka.

AUTHOR CONTRIBUTIONS
MNH, MS, AI, and MAH conceived and designed the overall study. MNH surveyed and collected all the field samples. MNH, RC, KG, OS, and OI carried out laboratory work including DNA extractions and sequencing, microbiological (cultural, biochemical) examinations, antimicrobial (antibiotics, metals) sensitivity tests, and biofilm assays. MAH and KC contributed chemicals and reagents. MNH, RA, and AI conceived, designed and executed the bioinformatics analysis. MNH interpreted the results and drafted the manuscript. MS, AS, KC, and MAH contributed intellectually to the interpretation and presentation of the results. Finally, all authors have approved the manuscript for submission.

FUNDING
The Bangladesh Bureau of Educational Information and Statistics (BANBEIS), Ministry of Education, Government of the People's Republic of Bangladesh (Grant No. LS2017313) partially supported this work.