A longitudinal study was carried out in dairy cattle farms in the Basque Country (northern Spain) to monitor the occurrence of ESBL-/AmpC-/CP-producing E. coli in apparently healthy animals. Five commercial farms (designated F1, F2, F3, F4, and F5), representative of the style of farming in the region, were selected to be enrolled in the study. Farms were located in the three counties of the Basque Country, and the distance between farms ranged from 15–25 km for those located within the same county (i.e., F3-F4 and F1-F2, respectively) and up to 160 km (F4-F5). Before the study started, our team paid a visit to each farm and, in the presence of the farm veterinary clinicians, farmers were interviewed face to face using a questionnaire that addressed general information about farm characteristics, management practices, vaccine programs, and antimicrobial drug use. Farm size based on the combined number of lactating and dry cows, heifers, and calves, ranged between 140 and 320 animals (mean = 240), with the number of lactating cows ranging from 75 (F5) to 200 cows (F1).

Monthly visits over a 1 year-period were planned for fecal sample collection. However, one of the farms (F5) dropped out after five samplings due to operational changes; samplings in the other four farms were interrupted midway through the study due to the COVID-19 pandemic and resumed at different times after the lockdown to complete the 12 samplings scheduled. Overall, the collection of fecal samples commenced in February 2019 and ended in October 2020, and extended over 16–17 months within individual farms. Samples were collected from apparently healthy animals from different age groups defined according to the different management practices, i.e., 1–5 month-old calves, 5–22 month-old heifers, and lactating cows. At each sampling time, rectal fecal samples (minimum of 5 g) were collected with a gloved hand from five animals randomly selected within each age group, and analyzed in a single 25 g pool per age group (5 g per animal). In seven time points, heifers could not be sampled in the two farms (five sampling times in F2 and two in F4) that raised heifer replacements at a breeding center. A total of 760 rectal fecal samples were collected and analyzed in 152 pools. Additionally, environmental slurry samples were also collected from F3 and F4 (two samplings each).

Upon arrival, samples were refrigerated at 4°C and sample processing was carried out within 3 days after collection, at the latest. Pooled fecal samples (25 g) were thoroughly mixed, diluted 1:10 in buffered peptone water (BPW, bioMérieux), and incubated at 37°C for 20 ± 2 h. For the isolation of ESBL-/AmpC-producing E. coli, two loops (20 μl) of BPW were subcultured onto MacConkey agar supplemented with 1 mg/l of cefotaxime and incubated at 37°C for 20 ± 2 h. Two morphologically different colonies per plate were harvested and confirmed as E. coli by species-specific real-time PCR detection of the uidA gene (Frahm and Obst, 2003).

For the isolation of CP-producing E. coli, two loops (20 μl) of BPW were subcultured onto MacConkey agar without antibiotics. A loopful of grown colonies was then harvested for DNA extraction and subjected to a real-time PCR amplification screening targeting the CP-coding genes bla_NDM, bla_VIM, bla_KPC, and bla_OXA-48 (Ellington et al., 2016). If any of these genes tested positive, a loopful of bacterial growth from the MacConkey agar was subcultured on ChromID^® Carba Smart selective agar plates (bioMérieux), and isolated colonies were identified by uidA gene detection as above.

Between 1 and 3 isolates per plate were selected and tested to assess antimicrobial susceptibility. Minimum inhibitory concentrations (MICs) were determined by broth microdilution using two Sensititre^® MIC susceptibility plates (EUVSEC1 and EUVSEC2, Thermo Fisher Scientific) following the recommendations in Commission Implementing Decision 2013/652/EU¹ concerning antimicrobials and dilution ranges, and the results were interpreted using epidemiological cutoff values (ECOFF). For antimicrobials with no ECOFFs assigned at the time, the results were interpreted as follows: for temocillin, ECOFF was fixed at 16 mg/l based on 2020/1729/EU; for azithromycin, 16 mg/l was used as a reference based on the bibliography (Sjölund-Karlsson et al., 2011; Clinical and Laboratory Standards Institute, 2015).

Based on their phenotypic AMR profile, sampling time, and age group isolation source, 41 isolates (27 from F4, 11 from F1, and one each from F2, F3, and F5) were selected for WGS. For in-depth genome characterization, genomic DNA was extracted from pure cultures using NZY Microbial gDNA Isolation kit (NZYtech) and subjected to long-reads (Oxford Nanopore Technologies, ONT) WGS. For ONT sequencing, a library was prepared using the Ligation Sequencing Kit (SQK-LSK109). Native barcoding genomic DNA kits (EXP-NBD104 and EXP-NBD114) were used for sample multiplexing except for three isolates that were sequenced in singleplex. Libraries were run in FLO-MIN106 (R9.4.1) or FLO-MIN111 (R10.3) flow cells on a MinION Mk1C device (ONT). For validation purposes, five isolates also underwent short-reads (Illumina) WGS; genomic DNA was submitted to Eurofins Genomics, where libraries were prepared based on the NEBNext Ultra II FS DNA library prep kit (Illumina) and sequenced with Illumina NovaSeq 6,000 (150-bp paired-end reads). The output files generated by ONT sequencing were basecalled in high-accuracy mode (HAC) and quality-filtered using Guppy (Qscore >7 in v4.2 and v4.3, and Qscore >8 in v5.0). Then, reads were adapter-trimmed and filtered by length and quality, as described before (Tello et al., 2022) and the resulting fastq reads were de novo assembled using Unicycler (Wick et al., 2017b). For one particular sample, Flye assembler (Kolmogorov et al., 2019) was used after retrieving inconsistent results in the draft genome generated with Unicycler, and the resultant assembly was the one further used in this study. For isolates sequenced by both technologies, Illumina reads were pre-processed for assembly as described elsewhere (Tello et al., 2020) and the outputs were further used to generate hybrid Nanopore-Illumina assemblies with Unicycler (Wick et al., 2017b). As previously described (Tello et al., 2022), isolates were subjected to in silico typing to determine their serogroup and phylogroup. MLST profiles were determined from unassembled long-reads using Krocus (Page and Keane, 2018). New sequence type (ST) assignations were obtained after submitting WGS reads to the Enterobase database (Zhou et al., 2020). Draft genomes were processed to predict plasmid- and chromosome-derived contigs using PlasFlow (v.1.1; Krawczyk et al., 2018). Molecular characterization of the isolates, including screening of ARGs, chromosomal point mutations associated with AMR, virulence factors detection, and plasmid replicon identification were performed, as previously described (Tello et al., 2020). Databases used for molecular characterization (ResFinder, PointFinder, PlasmidFinder, and ecoli_vf) were all updated on 20 October 2021. ResFinder hits were filtered at 90% coverage and identity and those with values below 100% were individually revised for frameshifts and amino acid changes, removing those considered not potentially functional. Virulence genes were filtered at 75% identity and 95% coverage, and the pattern of presence/absence of these genes was used as a typing scheme for genetic diversity. Genome annotations were carried out with Prokka (Seemann, 2014) and RAST (Aziz et al., 2008), and were graphically represented using SnapGene v.5.2.4.² Genome alignments were performed using MAUVE (Darling et al., 2010) in Geneious Prime v. 2020.2.4 software.³ Blast Ring Image Generator (BRIG) v.0.95 was used for plasmid structural comparison (Alikhan et al., 2011).

Phenotypic resistance profiles and the genetic determinants of resistance (GDR) in each sequenced sample (chromosome and plasmids) were represented in heatmaps. The plasmid heatmap was graphed along with a dendrogram illustrating the similarity among plasmids based on their AMR pattern. The hierarchical clustering analysis for the dendrogram was performed with the unweighted pair-group method with arithmetic mean (UPGMA) based on the Jaccard distance matrix, using the function hclust (v.3.6.1) of the R statistical package v.3.6.3. To identify the shared and unique phenotypic antimicrobial resistance profiles among the different age groups within each farm, Venn diagrams were constructed with the online tool InteractiVenn (Heberle et al., 2015).

To evaluate differences between age groups and farms in the shedding prevalence of cefotaxime-resistant E. coli and in the occurrence of phenotypic antimicrobial resistance for each antimicrobial, multivariate logistic regressions were performed including age group and farms as the explanatory variables. Adjusted odds ratios (OR_adj) were used as the measure of association between positivity and the explanatory variables and were expressed with their confidence interval at 95% (95% CI). Differences were considered statistically significant if p < 0.05. Simpson indices were estimated to calculate the diversity of phenotypic antimicrobial resistance profiles for each farm.

Following common practice in dairy farms in the Basque Country, all farms were closed production systems where replacement heifers originated from the same farm. Two of the farms (F2 and F4) raised their heifer replacements off-site in two different breeding centers. In both cases, animals leave the farm at 3–4 months of age and return already pregnant a few months before calving. A blanket antimicrobial treatment program was routinely used at dry-off that included the intramammary application of antimicrobials and teat sealant. The antimicrobials used for intramammary dry-cow therapy (DCT) were benzylpenicillin-benethamine/framycetin sulfate (Mamyzin) in F1 and F2, and cephapirin benzathine (Cefa-safe) in F3, F4, and F5. Farms participating in the study also used antimicrobials belonging to 12 antimicrobial drug classes for the treatment of disease in calves and cows. The antimicrobials most commonly used were third- and fourth-generation cephalosporins, followed by fluoroquinolones, and tetracyclines. Other antimicrobials used included penicillins, aminoglycosides, macrolides, and sulfonamides. Parenteral administration of fluoroquinolones was the most common treatment for mastitis during lactation in all except farm F2 where mastitis was not treated with antimicrobials. Third- and fourth-generation cephalosporins were the most common drugs used to treat reproductive diseases, diarrhea, and lameness that warranted systemic antimicrobial treatment. Reproductive diseases for which the producer opted to use antimicrobials included metritis, retained placenta, or other diseases related to reproduction.

Vaccination programs were quite different among farms. For example, vaccination against mastitis was only performed in F1. The vaccination program in F1 included vaccines against Infectious Bovine Rhinotracheitis (IBR), clostridia, and mastitis; IBR, Leptospira, and respiratory pathogens (Parainfluenza, bovine respiratory syncytial virus - BRS, Mannheimia) in F2; no vaccines at all were used in F3; IBR, Bovine viral diarrhea virus, and clostridia in F4; and clostridia, diarrhea in calves, and respiratory pathogens (Parainfluenza, BRS, Mannheimia) in F5.

Escherichia coli was isolated in cefotaxime-containing media in 92 of the 152 pooled fecal samples analyzed (60.5%) and in the 4 slurry samples collected from F3 and F4. These included samples collected from all farms and age groups, but differences in frequencies among age groups and farms were observed (Figure 1). Overall, isolation frequency of cefotaxime-resistant E. coli was higher in lactating cows [OR_adj = 4.71 (1.76–12.64), p = 0.002] and calves [OR_adj = 4.21 (1.59–11.18), p = 0.004] compared with heifers, and lower in F1 and F2 compared with the other three farms (LR χ² = 21.55, p < 0.001).

When available, between 1 and 3 cefotaxime-resistant E. coli isolates per age group and sampling date were selected in each farm for antimicrobial susceptibility testing. Thus, 187 isolates from fecal samples (73 isolated from lactating cows, 40 from heifers, and 74 from calves) and 10 isolates from slurry were analyzed. Since isolates had been obtained by selective isolation in a medium containing cefotaxime, they were all resistant to cefotaxime and ampicillin. Most isolates were also resistant to cefepime (99.0%) and ceftazidime (98.0%). Resistance to cefoxitin was detected in 36 isolates (18.3%), but 19 of them displayed a MIC value just one dilution step above the ECOFF. All 197 isolates were susceptible to tigecycline and colistin. Two isolates obtained from the same pool of feces collected from calves in F4 were resistant to all β-lactams tested, including temocillin and carbapenems (ertapenem, imipenem, and meropenem).

Co-resistance to other antimicrobial classes was also observed in most isolates (161/197, 81.7%) and 72.1% (142/197) showed multidrug resistance (MDR, resistance to 3 or more antimicrobial classes). Overall, resistance to tetracycline (53.8%), nalidixic acid (45.7%), ciprofloxacin (66.5%), sulfamethoxazole (69.0%), trimethoprim (48.7%), and chloramphenicol (47.7%) was very frequent, while resistance to gentamicin (29.9%) and azithromycin (14.2%) was lower and mainly associated to F5. The prevalence of resistance to each antimicrobial tested did not differ between age groups. However, statistically significant differences between farms were observed in the occurrence of resistance to several antimicrobials. Compared to other farms, F1 and F5 presented a significantly higher prevalence of tetracycline, chloramphenicol, and trimethoprim (all with p < 0.001). Resistance to gentamicin (p < 0.001), azithromycin (p < 0.001), ciprofloxacin (p = 0.002), and nalidixic acid (p = 0.009) was higher in F5 than in other farms, while resistance to cefoxitin was significantly higher in F1 and F2 (p = 0.003).

A total of 45 different profiles of microbiological resistance (Supplementary Table S1) resulting from the combination of antimicrobial agents that showed MICs above the ECOFF were observed in the study. Each phenotypic resistance profile was designated a letter of the Latin alphabet, and their distribution within each farm is represented in Figure 2. Within each farm, the number of different profiles ranged between 5 and 16 along the 12 samplings, the lowest diversity being found in F1 (Simpson index = 0.609) and the highest in F4 (Simpson index = 0.905). In F1, resistance to tetracycline, chloramphenicol, sulfamethoxazole, and trimethoprim remained stable during the entire study, whereas resistance to gentamicin, ciprofloxacin, and nalidixic acid was only observed in the second half of the study. This observation might reflect a shift in the circulating resistance profiles, where profile B, which dominated at the beginning of the study in all the age groups, was displaced by profile G in the second half of the study. On the contrary, the highest diversity in resistance profiles was observed in F4, where the three predominant profiles (A, C, and H) coexisted with 12 other profiles, with A and H dominating in the first half of the study, and profile C in the second half. Profile A only included resistance to ESBLs (penicillins and cephalosporins), whereas C and H included resistance to additional antimicrobials (Figure 2; Supplementary Table S1).

In F2, the prevalence of cefotaxime-resistant E. coli was the lowest, and fewer isolates were recovered and typed, particularly in heifers. Still, 11 different profiles were identified among 24 isolates, but profile D was the only one recovered in more than one sampling, in S7 in lactating cows and in S8 both in calves and lactating cows. In F3 high diversity in resistance profiles was observed, with a total of 16 different profiles, and a shift in the predominant resistance profiles occurred with time as happened in F4. Finally, the most outstanding feature of isolates recovered in F5 was the MDR pattern of all of them, with resistance to gentamicin and azithromycin being common in all age groups. On the other hand, isolates recovered from slurry samples shared their resistance profiles with isolates from fecal samples collected within the corresponding farms (Figure 2; Supplementary Figure S1).

ONT sequencing provided a median of 60,729 reads per sample (IQR = 22,781–405,021) in a median of 631 Mb per sample (IQR = 501–1,011 Mb) corresponding to a median coverage of 114X (IQR = 84X–182X; Supplementary Table S2). Upon assembly, the 5 isolates sequenced by both Illumina and ONT technologies, and 24 of the 36 ONT sequenced isolates resulted in circularized chromosomes. In all cases, the chromosome size of the assembled draft genome corresponded to the expected size of E. coli (median = 4,999,307 bp; IQR = 4,871,651 bp - 5,059,042 bp). Plasmid replicons were identified in a total of 125 contigs that in most cases (120/125, 96.0%) were assembled into complete circular plasmids. At least one plasmid replicon was identified in each isolate. IncF type plasmids were the most common (38/125, 30.4%), followed by IncB/O/K/Z (15/125, 12.0%), IncX1 (13/125, 10.4%), and IncY (13/125, 10.4%) along with 13 other replicon types. Screening for ARGs and SNPs associated with AMR identified 41 acquired ARGs and point mutations (9) in 4 other genes, coding for resistance to antimicrobials representing 9 different classes (Figure 3). The combination of GDRs detected in each isolate resulted in 22 different genotypic profiles of resistance (Supplementary Table S3). Sixty-two plasmids contained at least one ARG (Figure 4). None of the IncL, IncP, IncX4, or Col plasmids carried ARG genes.

A selection of isolates, mainly from farms F1 and F4, were analyzed by WGS to confirm whether the distribution of the different AMR profiles within the farms was due to different strains coexisting throughout the sampling period or to successive colonization by different strains. The 41 isolates were assigned to 18 MLST types, including two novel STs, i.e., ST-11626 in F4 and ST-12870 in F3. In F1, the 11 isolates tested belonged to 4 ST types, and 2 of them (ST-69 and ST-2930) included more than one isolate (Figure 3). Thus, ST-69 was represented by 5 isolates from the 3 age groups recovered at samplings S2, S4, and S7, which were identical in all other features, i.e., phylogroup (D), serotype (O15:H18), phenotypic resistance profile (B), and genotypic profile (a). ST-2930 included 4 isolates recovered in samplings S6 and S12 from calves and lactating cows that also shared all their genetic features, i.e., they were all assigned to phylogroup A, serotype O100:H25, and genotypic profile c. However, they split into two phenotypic resistance profiles differing only in susceptibility to FOX (profile E, 1 susceptible isolate, MIC_FOX = 2; and profile G, 3 resistant isolates, MIC_FOX = 16). The remaining two isolates sequenced (E0858 and E1072) were recovered during the second and last samplings, and had unique features.

In F4, 27 isolates were sequenced and assigned to 4 phylogroups (A, B1, C, and D) and 13 different ST types. Eight of the ST types were identified in more than one isolate (n = 2–4) and 5 were represented by a single isolate (Figure 3). As in F1, the most prevalent MLST type was ST-69. It included isolates recovered from lactating cows (S7), heifers (S9), calves (S11), and slurry. They all shared the same phylogroup (D), serotype (O15:H18), and genotypic (b) and phenotypic AMR profiles (C), suggesting that the same clone spread after sampling 7 in all animal groups and was also detected in slurry. As described in F1, differences associated with FOX among otherwise similar isolates were also observed in F4 within ST-23 (profiles AP and K) and ST-109 (profiles AB and Q). Other identical clones infecting several animals in F4 were those with MLST types ST-88 (n = 2), and ST-448 (n = 3).

On the other hand, several isolates with identical ST differed in other features (Figure 3). These included differences in resistance due to the carriage of ARG-harboring plasmids (i.e., ST-4981 and ST-69) and occasionally also in chromosomally encoded features (i.e., ST-58 isolates, which differed in ARG-harboring plasmids as well as serotype and chromosomally encoded ARGs). Finally, the isolates from F2, F3, and F5 were unique in all their features. When virulence genes were examined, patterns of presence/absence were highly conserved within ST types, with ST-58 and ST-69 being the only exceptions (Supplementary Table S4). Thus, for these epidemiologically related isolates, the typing scheme based on the presence/absence of virulence genes confirmed the diverse genetic profile inferred from the combination of all other features (phylogroup, serotype, and GDR profile).

ARG-harboring plasmids were present in all but two of the isolates (E1072 and E1027); 16 isolates carried a single plasmid and 23 carried 2 types of plasmids with ARGs. Overall, ARGs were present in 9 different types of plasmids, and 7 of them harbored ESBL-encoding genes, alone or in combination with several other ARGs (Figure 4). GDRs associated with ESBL production were only sporadically located in the chromosome. These included bla_CTX-M-14 (n = 1), bla_CTX-M-15 (n = 5), and the point mutation (nt 42 C → T) in the ampC promoter (n = 1). However, differences between farms were found concerning the diversity and location of these genes.

In F1, the most prevalent ESBL-encoding gene was bla_CTX-M-1 (9/11 isolates), present in isolates recovered from all animal groups and at different sampling times along the study (Figure 3). This gene was located in IncX1 type plasmids which always carried the same repertoire of ARGs. These plasmids were structurally compared using MAUVE and showed a high degree of similarity demonstrating the presence of largely conserved collinear coding blocks (data not shown). Hence, in addition to bla_CTX-M-1, IncX1 plasmids harbored the aminoglycoside resistance genes aadA2, ant(3″)-Ia, and aph(3′)-Ia, a trimethoprim resistance gene (dfrA12), a sulfamethoxazole resistance gene (sul3), and a chloramphenicol resistance gene (cmlA1; Figure 4). This plasmid was present in all isolates assigned to ST-69 and ST-2930.

In F4, a higher diversity of ESBL-encoding genes was observed (Figure 3). These included bla_CTX-M-14 (n = 11), bla_CTX-M-15 (n = 8), bla_CTX-M-27 (n = 3), bla_CTX-M-32 (n = 3), and bla_SHV-12 (n = 1). In addition, bla_NDM-1 was detected in one isolate (E1110). All were located in plasmids except for 4 chromosomally encoded bla_CTX-M-15. The bla_CTX-M-14 gene was always located in IncB/O/K/Z type plasmids that did not carry any additional ARGs (Figure 4). This plasmid was found in E. coli of different ST, genotypic and phenotypic profiles, isolated from slurry and animals of all age groups throughout the study. The bla_CTX-M-15 gene, which was also detected in all animal groups and environmental samples, was the predominant ESBL-encoding gene in isolates recovered in the second half of the study. This gene was located in IncY plasmids (all 4 ST-69 isolates) or in the chromosome (ST-4981 and ST-58). Besides bla_CTX-M-15, IncY plasmids harbored 6 other ARGs (Figure 4). The bla_CTX-M-27 gene was detected in the IncF plasmid of 3 ST-533 isolates, along with 5 other identical ARGs. Three isolates (ST-23) carried the bla_CTX-M-32 gene in an Incl1 plasmid, and the bla_SHV gene was present in an IncX3 plasmid in one isolate recovered from lactating cows in the last sampling. The isolation in F4 of a carbapenem-resistant E. coli harboring the bla_NDM-1 gene in an IncC plasmid was a significant finding extensively reported elsewhere (Tello et al., 2022).

In F1, ST-2930 isolates, besides IncX1, also carried an IncI1 plasmid that harbored another 5 ARGs [aac(3)-IVa, aph(3″)-Ib, aph(6)-Id, aph(4)-Ia, and tet(C)], while ST-69 isolates carried an IncF plasmid that harbored another four different ARGs [bla_TEM-1B, ant(3″)-Ia, sul1 and tet(A)]. A different IncF plasmid with a different collection of ARGs was present in one isolate (E0858). This isolate also carried a chromosomally encoded mutation in the ampC promoter. No plasmids were detected in the remaining isolate of F1 sequenced (ST-925), which harbored several GDRs in its chromosome, including bla_CTX-M-14 (Figure 3). Resistance to (fluoro)quinolones was always associated with point mutations in the gyrase and topoisomerase genes (gyrA, parC, and parE) and only observed in E. coli strains assigned to ST-2930 and ST-925 (Figure 3).

In F4, the 3 ST-23 isolates, which carried an Incl1 plasmid harboring the bla_CTX-M-32 gene, also carried an IncF plasmid, resulting in an identical genotypic profile. Instead, ST-533 isolates (E0888, E0892, and E0896) carried a second plasmid (IncHI2) with a different repertoire of ARGs (Figure 4). These IncHI2 plasmids were structurally compared (Supplementary Figure S2) and showed extensive sequence similarity, but E0896 lacked an 11.000 bp fragment that included aadA2, cmlA1, ant(3″)-Ia, and sul3 genes which was present in E0888 and E0892. Other genes coding only for resistance to narrow-spectrum β-lactamases like bla_OXA-1 (n = 3), bla_TEM-1A (n = 4), and bla_TEM-1B (n = 14) were mostly located in IncF plasmids, and less frequently in the chromosome or other type of plasmids such as IncX2 and IncY (Figures 3, 4). Resistance to (fluoro)quinolones in F4 was associated with point mutations in the gyrase and topoisomerase genes (gyrA, parC, and parE) in 10 isolates, and with the qnrS1 gene in another 10 (along with qnrB19 in one of them). Interestingly, the gene that codes for resistance to lincosamides, lnuF, was present in 5 isolates recovered from all age groups in the second half of the study. lnuF was always located in IncF and IncX3 plasmids (Figure 4).