Characterization and comparison of the microbiomes and resistomes of colostrum from selectively treated dry cows

Professionals in animal agriculture promote prudent use of antimicrobials to address public and animal health concerns, such as reduction of antimicrobial residues and antimicrobial resistance (AMR) in products. Few studies evaluate the effect of selective dry-cow therapy on preservation of the milk microbiome or the profile of AMR genes (the resistome) present at freshening. Our objectives were to characterize and compare the microbiomes and resistomes in the colostrum of cows with low somatic cell count that were treated or not treated with intramammary cephapirin benzathine at dry-off. From a larger parent study, cows on a New York dairy farm eligible for dry-off and with histories of somatic cell counts ≤200,000 cells/mL were enrolled to this study (n = 307). Cows were randomly assigned to receive an intramammary antimicrobial and external teat sealant (ABXTS) or sealant only (TS) at dry-off. Composite colostrum samples taken within 4 h of freshening, and quarter milk samples taken at 1 to 7 d in milk were subjected to aerobic culture. The DNA extraction was performed on colostrum from cows with culture-negative samples (ABXTS = 43; TS = 33). The DNA from cows of the same treatment group and parity were pooled (26 pools; ABXTS = 12; TS = 14) for 16S rRNA metagenomic sequencing. Sepa-rately, the resistome was captured using a custom RNA bait library for target-enriched sequencing. Sequencing reads were aligned to taxonomic and AMR databases to characterize the microbiome and resistome, respectively. The R statistical program was used to tabulate abundances and to analyze differences in diversity measures and in composition between treatment groups. In the microbiome, the most abundant phyla were Firmicutes (68%), Proteobacteria (23%), Actinobacteria (4%), and Bacteroidetes (3%). Shannon and richness diversity means were 0.93 and 14.7 for ABXTS and 0.94 and 13.1 for TS, respectively. Using analysis of similarities (ANOSIM), overall microbiome composition was found to be similar between treatment groups at the phylum (ANOSIM R = 0.005), class (ANOSIM R = 0.04), and order (ANOSIM R = −0.04) levels. In the resistome, we identified AMR gene accessions associated with 14 unique mechanisms of resistance across 9 different drug classes in 14 samples (TS = 9, ABXTS = 5). The majority of reads aligned to gene accessions that confer resistance to aminoglycoside (TS = ABXTS each 35% abundance), tetracycline (TS = 22%, ABXTS = 54%), and β-lactam classes (TS = 15%, ABXTS = 12%). Shannon diversity means for AMR class and mechanism, respectively, were 0.66 and 0.69 for TS and 0.19 and 0.19 for ABXTS. Resistome richness diversity means for class and mechanism were 3.1 and 3.4 for TS and 1.4 and 1.4 for ABXTS. Finally, resistome composition was similar between groups at the class (ANOSIM R = −0.20) and mechanism levels (ANOSIM R = 0.01). Although no critical differences were found between treatment groups regarding their microbiome or resistome composition in this study, a larger sample size, deeper sequencing, and additional methodology is needed to identify more subtle differences, such as between lower-abundance features.

In the resistome, we identified AMR gene accessions associated with 14 unique mechanisms of resistance across 9 different drug classes in 14 samples (TS = 9, ABXTS = 5). The majority of reads aligned to gene accessions that confer resistance to aminoglycoside (TS = ABXTS each 35% abundance), tetracycline (TS = 22%, ABXTS = 54%), and β-lactam classes (TS = 15%, ABXTS = 12%). Shannon diversity means for AMR class and mechanism, respectively, were 0.66 and 0.69 for TS and 0.19 and 0.19 for ABXTS. Resistome richness diversity means for class and mechanism were 3.1 and 3.4 for TS and 1.4 and 1.4 for ABXTS. Finally, resistome composition was similar between groups at the class (ANOSIM R = −0.20) and mechanism levels (ANOSIM R = 0.01). Although no critical differences were found between treatment groups regarding their microbiome or resistome composition in this study, a larger sample size, deeper sequencing, and additional

INTRODUCTION
The risk of an IMI is highest immediately postcalving relative to the remainder of a cow's lactation. This risk is associated with chronic subclinical infections present at dry-off or acquisition of a pathogen over the dry period (Todhunter et al., 1991;Bradley and Green, 2000;Green et al., 2002). In addition to the presence of an IMI at dry-off, cow and quarter-level characteristics, such as high milk production and poor teat end condition can also be associated with an increased risk of detecting an IMI within a week after calving (Dingwell et al., 2004;Robert et al., 2008). For these reasons, mitigation strategies developed and implemented since the 1960s have focused on blanket, or prophylactic treatment of all quarters of all cows with a long-acting antimicrobial at dry-off. In the United States, this was recently supported by a survey by the National Animal Health Monitoring System in 2014 that indicated that over 90% of cows are treated at dry-off and 90% of operations use antimicrobial drugs for this purpose (USDA APHIS, 2016). The overall use of antimicrobials on US dairies, including the large contribution of use at dryoff, has led to increased scientific interest and public awareness regarding the administration of preventative doses of antimicrobials to animals. This includes drug classes critical to human medicine such as third and fourth generation cephalosporins, which are currently on the World Health Organization's list of critically important antimicrobials for human medicine (WHO, 2020). Nearly all US cows treated at dry-off and approximately 75% of cows treated for mastitis receive a β-lactam antimicrobial. Further, of the cows treated for mastitis 50% receive a third-generation cephalosporin; for treated dry cows, the proportion is 28% (USDA APHIS, 2016). Concerns include the promotion of antimicrobial resistance (AMR) and the public health sequelae: potential of zoonotic infection with resistant pathogens and loss of effectiveness of antimicrobials for their intended use(s) (Threlfall et al., 2000;Aarestrup, 2005).
To explore the promotion of AMR in production animal systems, several studies have investigated the effects of antimicrobial use on tissue or manure microbiomes (microbial communities), as well as the presence of AMR genes in those sample types Rovira et al., 2019;Vikram et al., 2019). One system and product that has not been thoroughly evaluated is the conventional dairy farm and milk. Clinical trials on dairies have evaluated effects of antimicrobial use on the milk microbiome, whereas other studies have explored the association between antimicrobial use on dairies and AMR of clinical mastitis isolates in vitro (Pol and Ruegg, 2007;Lindeman et al., 2013). Metagenomic sequencing approaches, including of the hypervariable regions of the 16S rRNA gene, have been used to evaluate the effect of intramammary (IMM) antimicrobial use on the milk microbiome when cows were treated for clinical mastitis (Ganda et al., 2016 or at dry-off (Bonsaglia et al., 2017;Derakhshani et al., 2018). As compared with their untreated coun-terparts, infected or healthy quarters that were treated with IMM antimicrobial products were not different in their microbial profiles within days of treatment or within days of freshening (Ganda et al., 2016Bonsaglia et al., 2017). Only one study has performed deep sequencing of bovine milk with the objectives of characterizing the microbial organisms present and annotating the functionality or transcriptome of the sequences, including mechanisms of resistance (Hoque et al., 2019). Although the study provides differential assessment between a limited number of healthy and mastitic samples, it does not evaluate the effect of antimicrobial use on the microbiome or resistome. For these reasons, we chose to focus our sequencing strategy on one area of high IMM antimicrobial use on dairies, at dry-off.
To this end, we characterized the microbiome and resistome of colostrum from healthy cows, as within a selective treatment protocol, these animals would not receive IMM treatment. Within this study, pooling of colostrum microbiome samples before DNA amplification allowed us to evaluate community-level diversities and AMR; pooling is currently considered a viable measure in population-level association research studies (Ray et al., 2019). Our first objective was to determine if the microbiome in pooled colostrum of cows treated with antimicrobials at dry-off differs from cows that are not treated with antimicrobials at dry-off. Our second objective was to explore and compare the presence of AMR genes between the pooled colostrum groups. Differences in the microbiomes detected and the AMR genes that they possess might further elucidate the effect of prophylactic use of antimicrobials at dry-off.

Study Design
A subset of 76 cattle colostrum samples from a splitherd randomized clinical trial were subjected to DNA extraction, 16S rRNA metagenomic sequencing, and enriched shotgun sequencing for characterization of the microbiome and presence of antimicrobial gene accessions, respectively.

Study Animals and Project Farm
A randomized clinical trial was performed from June 2016 to March 2017, at a New York State commercial dairy milking 1,800 Holstein cows (Vasquez et al., 2018). The study, which explored a selective dry-cow treatment protocol, was approved by Cornell University Institutional Animal Care and Use Committee (approval #2013-064). Cows with the following crite-ria were eligible for dry-off, and therefore, enrollment into the project: pregnant >220 d or pregnant >180 d and producing <11.4 kg of milk per day. This farm used DHIA testing, which included once per month SCC. DairyComp 305 (Valley Ag Software) was used for recording mastitis events. All cows included in the current study subset (n = 307) had a SCC ≤200,000 cells/mL on the most recent test, an average SCC over the last 3 mo of ≤200,000 cells/mL, and no more than 1 clinical mastitis event in the current lactation. Any clinical event must not have occurred in the last 14 d and the projected dry period for each cow had to be less than 100 d. Within the larger trial, these low-risk cows were then randomly assigned to receive IMM antimicrobials (ToMorrow; Boehringer Ingelheim) and an external teat sealant (ABXTS; T-hexx Dry; Hydromer Inc.) on all 4 quarters or application of the teat sealant only (TS) to all 4 quarters at dry-off. Only cows with negative milk aerobic culture results at 1 to 7 DIM along with a negative composite colostrum sample were eligible for inclusion in the subset (n = 76; Figure 1).

Sample Collection
Quarter-level milk samples and composite colostrum samples were taken from cows at 2 time points dur-  (Vasquez et al., 2018) into the current study. ABXTS = antibiotics; TS = teat sealant. ing the trial: at freshening (composite colostrum), and at 1 to 7 DIM before milking (quarter level). In this investigation, quarter-level samples were used for screening, whereas the colostrum samples were used for sequencing. As colostrum samples would be submitted for sequencing, extra measures were taken to ensure colostrum samples were not contaminated. Cows were placed into restraining stocks and visible debris was removed from the udder. Teats were dipped with an iodine-based premilking teat disinfectant and wiped with individual paper towels. For both quarter-level and colostrum samples, teat barrels and teat ends were scrubbed with individual 70% isopropyl alcohol-soaked gauze and allowed to dry. This was repeated until no visible dirt was present on the teat end or gauze. Each quarter was then stripped twice, and the milk discarded before sample collection. Then the desired volume of colostrum was retrieved into a hinged-top sterile container (5 mL for quarter-level samples, 40 mL for composite colostrum samples). Quarter-level samples were retrieved by study personnel, and composite colostrum samples were retrieved within 4 h of calving by 2 regularly trained on-farm staff. All sample tubes were marked with the date and the cow number, as well as with the initials of the person sampling. Samples were immediately placed in a −20°C freezer for storage and transported in an insulated cooler to the university laboratory freezer (−20°C) on a weekly basis.

Culture and Pooling Strategy
Colostrum samples were thawed at room temperature, inverted several times, and aliquoted into 2, 15-mL conical tubes using sterile technique. One tube was refrozen in a −20°C freezer, and the second was used to perform aerobic culture according to National Mastitis Council guidelines for identification of aerobic organisms (National Mastitis Council, 2017). A sterile cotton swab was used to streak 0.03 mL of milk on trypticase soy agar containing 5% sheep blood and 1% esculin (PML Microbiologicals). Culture plates were incubated aerobically at 37°C for up to 48 h and any samples with colony growth were documented. Additionally, quarter-level samples were subjected to aerobic culture at Quality Milk Production Services (Ithaca, NY) using the same techniques. Identification of specific mastitiscausing organisms was performed using the morphological (colostrum and quarter-level samples) and analytical methods (quarter-level samples) outlined in Vasquez et al., 2018. Only colostrum from cows with negative results for all samples were retained for further analysis. A cow with an IMI at dry-off or freshening is not indicative of a healthy microbiome and a dominant pathogen would influence the analysis. Thirty-three ABXTS cows and 43 TS cows qualified for inclusion into the study (total = 76).
To maximize the information obtained via sequencing for the financial means available, a pooling strategy was developed. Individual samples were aggregated into 26 pools based on treatment group, parity, antimicrobial and mastitis history, dry period length, and DIM at dry-off (TS = 14, ABXTS = 12; Supplemental Table S1; https: / / github .com/ avasquez321/ supplementalmaterials; Vasquez et al., 2021). Primary comparisons included differences in the microbiomes and resistomes between the treatment groups.

Extraction, Amplification, and Sequencing Methods
Extraction was performed on a 50-µL aliquot of each of the 76 composite colostrum samples using DNeasy PowerFood Microbial Kit (Qiagen), according to manufacturer's instructions with one modification. To maximize yield, purified DNA was eluted in 50 µL of elution buffer (10 mM HCL) and passed through the spin filter twice. The DNA concentration and purity were determined with the Qubit dsDNA HS Assay Kit using the Qubit 2.0 Fluorometer and NanoDrop One/ One Microvolume UV-Vis Spectrophotometer (Thermo Fisher Scientific). Each volume of extraction (50 µL) was pooled using the strategy outlined in Supplemental Table S1 (n = 26), with pools totaling 50 to 250 µL of amplicons each. Extraction of DNA from blank (negative control) samples was not performed.
The 16S quantitative PCR was used to quantify the bacterial concentration of each pool. A 20-µL reaction was used, containing 10 µL of 2xSYBR master mix (Thermo Fisher Scientific), 1 µL of 4 µM forward and reverse primers, 8 µL of molecular grade water, and 1 µL of DNA template. Femto Bacterial DNA standards (Zymo) served as a 4-point, 1:10, standard curve for quantification of samples. Two-step PCR was used with the following parameters: 95°C for 10 min, 40 cycles of 95°C for 15s, and 60°C for 45s (Witzke et al., 2020). A melt curve was generated to evaluate for nonspecific products. Primer oligos were taken from the NEBNext DNA Microbial Enrichment kit (New England Bio-Labs Inc.), and the sequences were as follows: forward (5′-CCATGAAGTCGGAATCGCTAG-3′) and reverse (5′-GCTTGACGGGCGGTGT-3′).
Extremophiles were then spiked into each sample at 1% of quantitative PCR-determined bacterial content. Two strains were used in a 2:1 ratio, Deinococcus radiodurans: Salinibacter ruber. Isolated gDNA (25 µL) containing the extremophile spike-in was sent to the Novogene Corporation for 16S rRNA gene sequencing. The hypervariable V4 region of the 16S rRNA gene was amplified using universal 16S primers 515F (5′-GTGC- Vasquez et al.: COLOSTRAL MICROBIOMES AND RESISTOMES AFTER SDCT CAGCMGCCGCGGTAA-3′) and 806R (5′-GGAC-TACHVGGGTWTCTAAT-3′) for characterization of the microbiome. The library of amplicons was sequenced on the HiSeq 2500 Sequencing System using on the Illumina platform (Thermo Fisher Scientific) with a targeted depth of 100,000 reads per sample following Novogene's protocol and quality control standards. As a negative control, equal volumes of nuclease-free sterile water were included with each batch of samples processed. These controls were included in procedures for PCR amplification and in preparation of sequencing libraries. Sequencing of these negative control samples did not yield product and therefore were not included in downstream analysis.
The remaining volumes of extracted pooled amplicons were subjected to target enrichment for sequencing of the resistome. The SureSelect XT-HS Target Enrichment Kit (Agilent) was used for library preparation, hybridization, and capture of resistant genes (resistome) via a custom AMR bait library, MEGaRICH, as developed by Noyes et al. (2017). Modifications to the library preparation protocol were performed as follows: the shearing methods were altered to produce larger fragments of 250 to 350 bp. This was done by reducing the Covaris M220 (Thermo Fisher Scientific) peak incident power to 50 and shearing 500 ng gDNA in 130 µL of 10 mM Tris-EDTA (TE) buffer (10 mM Tris, 1 mM EDTA) for 240 s. Samples were then purified using AMPure Beads (Beckman Coulter Life Sciences) in a 0.9:1 ratio (117:130 µL) and resuspended in 52 µL of nuclease-free water. End-repair, a-tailing, and ligation were all performed according to manufacturer's instruction. Pretarget enrichment amplification was done with 7 cycles, hybridization was performed at 70°C and post-target enrichment amplification was then done with 13 cycles. The PCR cleanup was performed using AMPure XP beads in a 1:1 ratio. One microliter of amplicons was used to confirm the quality of library via Tapestation 2200 (Agilent). The 26 pools were then pooled at equal nanomolar amounts (7 nM) and size selected via Blue Pippin (Sage Science) before submission to the University of Colorado, Denver, Anschutz. Sequencing was performed using the Novaseq 6000 with a paired-end 150 cycle (2 × 150) on the Illumina platform (Thermo Fisher Scientific) following the standard protocols and quality control standards. Reagent-blank negative control samples were not included in preparation of target-enriched sequencing libraries. Sequencing depth was targeted at 50 million reads per sample.

Quality Filtering and Alignment
For the microbiome, 16S reads were processed, aligned to taxonomic databases, and analyzed using the Quantitative Insights Into Microbial Ecology version 2 pipeline (Bolyen et al., 2019). Briefly, DADA2 was used for quality filtering and denoizing (Callahan et al., 2016) and taxonomic classification was performed using a naïve Bayes classifier trained on the GreenGenes database (McDonald et al., 2012). Results for each pool were tabulated into count tables and imported into the R statistical program version 3.5.1 (https: / / r -project .org) for downstream analysis.
For the characterization of the resistome, the sequencing reads from target enrichment were analyzed with the AmrPlusPlus 1.0 workflow . Briefly, sequence quality control was performed using Trimmomatic to remove low quality reads and adaptor contamination (Bolger et al., 2014). The Burrows-Wheeler-Aligner (BWA) software and SAMtools were then used to map sample reads to a reference Bos taurus genome and to filter out host-associated contaminant reads (Li et al., 2009;Merchant et al., 2014). Next, the MEGARes v1.02 database, which contains sequence data for approximately 4,000 hand-curated AMR genes , was used with BWA for alignment and classification of these reads to antimicrobial resistance gene accessions. To reduce false-positive classifications of AMR gene accessions, only gene accessions for which >80% of the reference sequence was aligned to by at least one read was included in all downstream analysis.

Statistical Analysis
Summary statistics including tabulations and comparisons of the numbers of reads for each sample were compared using Wilcoxon Rank-sum tests with the wilcox.test function in R. For study design metadata, primary comparisons of interest were between ABXTS and TS, but any effect of lactation group was also investigated. No statistical differences were found between pools of different lactations (not reported). The counts for the extremophiles were identified and extracted from the samples before statistical analysis.
The MetagenomeSeq R package (Paulson et al., 2013) was used for count normalization with cumulative sum scaling (CSS) and differential abundance testing to compare treatment groups was done using zero-inflated Gaussian models (Paulson et al., 2013;Doster et al., 2018). Relative abundances of microbiome and resistome features in each sample were tabulated and analyzed at different taxonomic levels , but data in this manuscript is reported at the class and mechanism levels of the resistome and to the phylum, class and order levels of the microbiome. This reduces bias associated with inconsistent naming criteria, unclassified reads, and unreliability of results at lower classification levels (Caro-Quintero and Konstantinidis, 2012;Hall and Schwarz, 2016). For plotting of the microbiome, the phyla and classes amounting to greater than 2% relative abundance were compiled; those with lower abundances were summed into a category of "other." For the resistome, heatmaps were constructed using the log 2 transformed CSS counts and plotted using a gradient scale. Alpha diversity comparisons for both the microbiome and resistome (richness, the number of unique bacteria or mechanisms of resistance and Shannon index, accounting for number and proportion of the same) were calculated using the vegan package in R (Oksanen et al., 2016). Statistical comparisons of diversity indices were performed using nonparametric Wilcoxon Rank-sum tests with the wilcox.test function. The vegan R package was also used to compare congruence of the microbiome and resistome ordinations using nonmetric multidimensional scaling (NMDS) on Euclidian distances that were calculated on Hellingertransformed counts (using the metaMDS function in R). Differences in microbiome and resistome composition clustering between ABXTS and TS were tested with analysis of similarities (ANOSIM). P-values are presented for all comparisons.

RESULTS
Raw sequencing data for all pools can be accessed in the National Center for Biotechnology Information repository, BioProject number PRJNA724672. Before pooling, the range of DNA concentration per sample was 3.4 to 119.8 ng/µL (mean = 29.2 ng/µL, SD = 30.3 ng/µL).
Resistome composition was similar between groups at the class (ANOSIM R = −0.20, Stress = 0.11, P > 0.99) and mechanism levels (ANOSIM R = −0.01, Stress = 0.06, P = 0.5; Table 2). No individual AMR classes were differentially abundant between drug classes (All adjusted P > 0.1). Of the counts conferring to aminoglycoside resistance (35%), nearly all were classified as encoding aminoglycoside O-phosphotransferases (98% of counts within this class). These were not differentially abundant between treatment groups (P = 0.22), however, all other mechanisms that were present, such as class A β-lactamases and tetracycline resistance mechanisms were differentially abundant between groups (P < 0.02). The TS group pools had greater relative abundances compared with the ABXTS samples for all but the tetracycline resistance ribosomal protection proteins and tetracycline transcriptional repressors. Within the tetracycline class, most were classified as encoding tetracycline resistance ribosomal protection proteins (72%). Complete counts, including those for specific genes can be seen in Supplemental Table S3 (https: / / github .com/ avasquez321/ supplementalmaterials; Vasquez et al., 2021). The NMDS ordination depicted overlap for both class and mechanism ( Figure 4); treatment groups did not cluster apart. Stress values for these ordinations are provided in Table 2.

DISCUSSION
The results of this trial suggest that the effect of antimicrobial use at dry-off on colostrum microbiomes as well as on the presence of AMR genes in colostrum is minor compared with other unmeasured factors. We predicted that β-lactam resistance genes in the colostrum resistome would be enriched in the presence of direct selection pressure from cephaparin benzathine use. Consequences of such findings might have resulted in scrutiny of current antimicrobial use practices and subsequent motivation for prudent use to prevent the increase and transfer of critical AMR genes from livestock farms. In a previous study, enrichment of AMR genes in manure was found to be associated with tetracycline use in feedlot systems, supporting the implementation of new veterinary feed-directive regulations, which eliminate use of antimicrobials in feed without direction of a veterinarian . Alternatively, recent studies on retail ground beef products found that the resistomes and microbiomes of products with "raised without antibiotics" labeling were similar to products without these claims; products differed more by production system and by supplier than by label type (Vikram et al., 2018;Doster et al., 2020b). The current study, which investigated colostrum, did not find any statistical differences in the overall composition of microbes or resistomes between pools from dairy cows treated with an IMM cephalosporin product versus pools from cows not treated with an antimicrobial at dry-off. Though the overall composition of the microbiomes and resistomes did not differ between groups, our analysis found nonappreciable but statistical differences in abundances when low abundance taxonomic classes, orders, and microbial AMR mechanisms were evaluated. Effect of IMM antimicrobial use on the microbiota of other organ systems and the subsequent contribution of AMR genes via environmental excretion was not evaluated.
Though not ascertained in the current study, differential findings in the microbiomes of colostrum between treatment groups could suggest that administration of antimicrobials at dry-off facilitate the creation and maintenance of a dysbiosis, or disruption of the healthy colostrum microbiome. The shortcomings of using 16S sequencing for milk include a low yield of biomass and the inability to confirm that the microbial sequences amplified represent an active population within the mammary gland; the idea of a milk microbiome is challenged in a recent opinion article (Rainard, 2017). Notwithstanding, studies using a variety of sampling and sequencing methods are supportive of a standard, yet small population of commensal bacteria in bovine milk that is impacted by a health event such as mastitis (Oikonomou et al., 2014;Addis et al., 2016;. Additionally, findings from human, mouse, and bovine studies have each suggested cause and effect relationships between microbiomes and health status (Ley et al., 2006;Krajmalnik-Brown et al., 2015;Lima et al., 2017). These studies showed that, independent of antimicrobial use, a disruption or alteration of the normal gastrointestinal or colostral flora is associated Vasquez et al.: COLOSTRAL MICROBIOMES AND RESISTOMES AFTER SDCT Table 2. Analysis of similarity comparing the microbiomes and resistomes of pools of colostrum from cows receiving cephapirin benzathine and teat sealant at dry-off (ABXTS) to pools of colostrum from cows receiving teat sealant only at dry-off (TS) with subsequent disease; specifically, autism, diabetes, and mastitis. In regards to the milk microbiome, groups continue to investigate whether these populations are functional in combating pathogenic inhabitants, as was shown to be the case for the human breastmilk microbiome against Staphylococcus aureus (Heikkilä and Saris, 2003). This functionality would support the finding that shifts in microbiota resulting from antimicrobial use can increase susceptibility to disease (Willing et al., 2011). In the case of the milk microbiome and IMM antimicrobial administration in dairy cattle, no shortterm shifts were documented for in healthy nor mastitic quarters treated with IMM antimicrobials nor were any longer-term health outcomes observed (Ganda et al., 2016;Ganda et al., 2017;Vasquez et al., 2019). Though our cross-sectional study did not investigate longitudinal changes, it is congruent with the previous studies in that no differences in the overall microbiomes were , and drug mechanism (D) levels using nonmetric multidimensional scaling (NMDS), between pools of colostrum from cows treated with cephapirin benzathine and external teat sealant at dry-off (blue; n = 12 pools for microbiome, 5 pools for resistome) and those only treated with external teat sealant (pink; n = 14 pools for microbiome, n = 9 pools for resistome). Individual dots represent pools. AMS = antimicrobial resistance; MDS1 and MDS2 = multidimensional scales 1 and 2.
found between groups when samples were taken from treated versus not treated cows at later postadministration time points.
Caution should be taken when comparing results of our study to others. Our samples are quarter-level, pooled composites containing multiple animals, which as a method has increased risk for contamination and can preferentially increase abundance of phyla such as Clostridiales (Andrews et al., 2019). A few previous studies have investigated the effect of using long-acting Vasquez et al.: COLOSTRAL MICROBIOMES AND RESISTOMES AFTER SDCT Figure 5. Heat map of log-normalized counts of reads aligning to antimicrobial resistance classes (A) and mechanisms (B) in pooled colostrum from cows treated with antibiotics and external teat sealant at dry-off (ABXTS; n = 5 pools) and those only treated with external teat sealant (TS; n = 9 pools). Abundance calculated from normalized counts (cumulative sum scaling) of taxonomically classified reads. "Other" is a composite of all phyla individually representing <2% of counts. MLS = macrolides, lincosamides, and streptogramines; MFS = major facilitator superfamily.
antimicrobial formulations at dry-off on the milk microbiome. Derakhshani et al. (2018) sequenced DNA extracted from late-lactation milk and colostrum from 9 treated cows to characterize the microbiome over the dry period. Similar to our study, more than 75% of reads in colostrum were classified to the phyla Firmicutes and Proteobacteria. Bacteroidetes and Actinobacteria were also the third and fourth most commonly represented phyla. In regards to longitudinal trends, the authors observed a high degree of commonality between the microbiota of both time periods, again supporting resiliency of the natural populations. However, they also concluded from the few alterations found, that there was an effect of IMM penicillin G/novobiocin administration on the microbiota profiles. Differences were found between the dry-off milk samples and the colostrum samples for relative abundances of Firmicutes and Proteobacteria, Chao diversity indices, and proportional differences of specific genera. Conversely, in the absence of a control group, the association between antimicrobial use and changes in microbiome is not clear; differences could be attributed to physiological changes over time, which are well known to occur during the dry period (Sordillo and Nickerson, 1988;Hurley, 1989;Capuco et al., 1997). In the current study, differences (2 and 3%) were found between treatment groups for the Clostridia and Deltaproteobacteria classes. Clostridiales, within the Clostridia class, is part of the core milk microbiome (Taponen et al., 2019); therefore, one would expect it to be in greater abundance in the nontreated group rather than the treated group. These results should be interpreted lightly, however, as we did not find any clustering in the overall microbial composition between groups. Influential pools or the samples within can also drive these differences between treatment groups. One TS pool, for example, had a very low abundance of Clostridiales (Supplemental Table S2).
The overall relative abundance of Firmicutes (34%) was appreciably lower in Derakhshani et al. (2018) than in our study (>60%), likely due to differences in farm location, management practices, or laboratory methods. Use of a narrow spectrum antimicrobial versus use of a broad-spectrum cephalosporin may also have led to preferential treatment and elimination of primarily gram-positive microflora in the previous trial. However, a broad-spectrum dry-off product was administered in Lima et al. (2017) and colostrum samples also had appreciably lower abundances of Firmicutes (40.8%) than in our study . Still, the 4 most commonly represented phyla again were Firmicutes, Bacteroides, Proteobacteria, and Actinobacteria. Finally, no effects on the milk microbiome were described by Bonsaglia et al. (2017) when analyses included numbers of reads, Shannon and Chao diversity indices, and ANOSIM for cows treated with an IMM cephalosporin. The authors also account for Firmicutes, Bacteroides, Proteobacteria, and Actinobacteria phyla amounting to >90% of reads (Bonsaglia et al., 2017). Perhaps no differences were found between groups due to the resiliency of milk microbiomes as documented in the aforementioned studies, or perhaps because of the short duration of antimicrobial use for IMM treatments. Relevantly, when feedlot cattle were treated with an injectable cephalosporin product that maintains effective tissue concentrations for up to 2 wk (FDA, 2003), differences were observed between diversities of microbiomes in manure from high versus low antimicrobial exposure feedlot pens . As the IMM product used in the current study is intended to have prolonged activity after dry-off, lack of antimicrobial effect is more likely due to a combination of microbiome resiliency and study design. Average dry periods in the current study were approximately 55 d, allowing for return of natural populations following treatment and before sample collection. The timing of sample retrieval is certainly a limitation of our study as serial samples immediately after dry-off might have allowed us to characterize microbial changes. However, producers must follow withdrawal periods postadministration for these products; evaluating the effect in early lactation has more benefit for both fresh-cow and public health concerns.
Though no important alterations were found in the microbiomes of pools of treated cows, the contribution of antimicrobial residues, the potential enrichment of AMR genes, and the subsequent transfer of mechanisms between microbiota might represent a more significant concern (Sommer et al., 2010). We did not observe any important differences between resistome composition in the colostrum of treated versus nontreated cows. Despite targeted enrichment and a targeted sequencing depth of 50 million reads per sample, we found zero AMR accessions in many of our pools (12/26), and our overall alignments to gene accessions per sample was small (368 counts/sample) compared with studies evaluating manure at a depth of ~25 million reads (61,751 counts/ sample) and ground beef at a depth of 100 million reads (8,372 counts/sample; Zaheer et al., 2018;Doster et al., 2020b). In the current analysis, AMR gene accessions represented 0.6% of bacterial metagenomic sequences on average, which is in line with other studies (Noyes et al., 2016;Rovira et al., 2019). These findings suggest that AMR genes are in low abundance in colostrum compared with other types of samples, indicating that the risk of transmission of these genes via the consumption of raw colostrum, as with calves, would be low.
Deep whole-genome sequencing has been performed on milk to elucidate the "natural presence" of AMR mechanisms. Hoque et al. (2019) thoroughly evaluated microbial, viral, and archaeal composition as well as proteomic characterization of healthy and mastitic milk samples (n = 21). The study found a variety of resistance genes including multidrug resistance operons, fluoroquinolone resistance mechanisms, aminoglycoside adenylyl transferases, β-lactamase resistance, and others. Considering that cows were not treated with antimicrobials before sampling, the study suggests that the presence of AMR genes is not dependent on recent treatment events. Pathogenic and nonpathogenic bacteria within samples can harbor these genes, albeit a rare proportion of the total reads (Kazimierczak et al., 2009). Natural presence could certainly explain the counts mapping to the nonadministered, antimicrobial classes not approved by the Food and Drug Administration, in our samples. Fluoroquinolones and aminoglycosides are not approved for use in lactating dairy cows, yet they are present in one and 9 pools, respectively. Other reasons for presence of these drug classes are co-selection. In this study only 1 of 9 pools had hits to the aminoglycoside class alone, whereas the remainder had counts to multiple classes, supporting the idea of co-selection. Our findings are similar to other studies that have found AMR accessions in animal manure from facilities where animals were raised without antimicrobials (Santamaria et al., 2011;Vikram et al., 2017;Vikram et al., 2018). Finally, antimicrobial use in other production groups at the same facility can also contribute to the population of AMR genes in the farm environment and within the samples.
Pooling in our study limits the ability to determine if the colostrum from one or more cows was influential to our findings. Though treatment group assignment dictated whether a cow received antimicrobials at dry-off or not, some of the cows in this study received antimicrobials at another point during the lactation. This was the case for 5 TS pools and 4 ABXTS pools (Supplemental Table S1). Unfortunately, given the small sample size it is also difficult to determine if there is an association between previous antimicrobial use and reads aligning to the AMR database. Of the 5 TS pools containing cows with histories of antimicrobial use, 4 had counts for AMR accessions and 2 of 4 ABXTS pools with histories of use did. However, we found AMR accessions in 5 pools containing cows with no documented exposure to antimicrobials. This again supports the phenomenon that acquisition of AMR genes does not depend on antimicrobial usage in individual animals alone, and that many classes of resistance are ubiquitous (D'Costa et al., 2011;Stokes and Gillings, 2011;Cadena et al., 2018).
One should note that the presence of AMR gene accessions does not necessarily equate to expression or ability to confer resistance (Forslund et al., 2013). Though the extra step does not confirm expression of the gene, a previous study by common authors have added SNP confirmation within their methods to screen for the presence of SNPs within the generated sequences that map to the database. The understanding is that a SNP at the wrong location might limit the ability to provide resistance (Doster et al., 2020a). This was not performed within the current study, but counts were retrospectively compiled and evaluated. There were a total of 1,373 counts mapping to genes within AMR classes with suggested SNP confirmations, 50% of which were contributed by 1 TS sample (Supplemental Table S4; https: / / github .com/ avasquez321/ supplementalmaterials; Vasquez et al., 2021). Overall 25 and 15% of the suggested genes mapped to fluoroquinolone and rifampin drug classes, respectively. As these drug classes were not the primary focus of our study, we did not further investigate these counts.
One statistical difference we described was for diversity indices between treatment groups, specifically, for Shannon indices in AMR class and mechanism. This demonstrates the influential nature of one or 2 pools and the limitations of a small sample size. Two of 9 pools in the TS group had counts mapping to 6 and 9 different drug mechanisms (4 and 8 classes) whereas none of the 5 pools represented by the ABXTS group had counts mapping to more than 3 mechanisms or classes. Influential pools could also explain why differences were found when comparing specific drug classes or mechanisms between groups, particularly for drug classes such as class A β-lactamases, which had higher abundances in the TS group pools. A larger sample size would be beneficial as substantial variation in the microbiome or resistome can exist between pools. Within a study, such variation has been described among clinically healthy quarters, between different cows from the same farm, and between different product processing locations (Oikonomou et al., 2014;Falentin et al., 2016;Ganda et al., 2017;Vikram et al., 2018). The effect of sample size on statistical power is poorly understood in metagenomic and AMR sequencing studies and to date, there are no recommended guidelines for calculating sample sizes, other than increasing the number of samples as natural variation increases within a study population. In this study, sample size was limited by cost and pooling was performed to maximize the information retrieved.

CONCLUSIONS
There are limited number of approved IMM antimicrobials available in the United States for dry-cow administration, and all are β-lactam antimicrobials. We Vasquez et al.: COLOSTRAL MICROBIOMES AND RESISTOMES AFTER SDCT evaluated the effect of administration of a first-generation cephalosporin, cephapirin benzathine, on the microbiome and resistome of colostrum samples from low SCC cows on one farm. The 4 dominate phyla making up the microbiome were Firmicutes, Proteobacteria, Actinobacteria, and Bacteroidetes. There was a low prevalence of AMR accessions with aminoglycosides, tetracyclines and β-lactams being the most frequent classes. Exploration of the subtle differences between treatment groups for specific taxonomic levels is needed to determine the effect of sample size and to investigate biological relevance. While overall differences in the microbiome and resistome between treated and untreated cows were not found in this study, a larger sample size, deeper sequencing, and additional methodology would help to confirm our findings.