We used data of monthly quarter level milk samples from all lactating cows in two Danish dairy herds throughout 1 year between January 2017 and January 2018. The herds had an average number of milking cows of, respectively, 180 and 360 cows and both herds had side-by-side milking parlors. Herd 1 was milked two times per day, and Herd 2 was milked three times per day. The bedding material in both herds was straw. Herd 1 had a mean bulk milk SCC of 294,000 (range = 261,000–324,000 cells/mL) and Herd 2 had a mean bulk milk SCC of 280,000 (range = 236,000–299,000 cells/mL) (15). We selected the length of the sampling interval to mimic the routine performance samples that take place approximately 11 times annually (Dairy Herd Improvement sampling). In total, we sampled 214, and 443 lactating cows from Herd 1 and 2, respectively (Table 1). A detailed description of the herds and the sampling procedure can be found in Kirkeby et al. (15).

In Herd 1, the cows were housed in free stalls with mattresses, and straw was used as bedding material. Lactating cows were milked two times a day in a parallel parlor of 2 × 12 cows. The procedure for milking was according to the NMC recommendations using gloves and applying pre and post-dip. In Herd 2, the lactating cows were milked three times a day also in a parallel parlor, here in 2 × 16 cows. At milking, the personnel used gloves, and teats were wiped off using clean cotton towels, but no post-milking teat dipping was implemented.

The collection of quarter milk samples was performed in line with NMC standards (http://www.nmconline.org/sampling.htm). As described in Kirkeby et al. (15), the teats in Herd 1 were sprayed with a foaming teat wipe-off product (Viri Foam,Novadan ApS, Kolding, Denmark), whereas the teats in Herd 2 were predipped using milk wash from Trinol (Hobro, Denmark). The teats were then cleaned with cotton towels soaked in water and a minimum of 4 squirts of milk were discarded from each quarter. After this teats were sanitized using single-service wet wipes (MS Lavettes, MS Schipper, Bladel, the Netherlands) soaked in 90% ethanol, with one towel per teat. The teats were further cleaned before milk sampling using 90% ethanol spray followed by air-drying for a minimum of 30 s. For sampling, pre-labeled sterile sample vials were applied. After filling of the vials with 15 ml milk, they were recapped immediately. All samples were stored in thermally insulated boxes, with cooling elements shipped to the laboratory. The samples were then preserved using 0.5% boric acid, and bacteriological analysis took place within 48 h of collection. Bacteriological analysis was conducted as explained in Kirkeby et al. (16). Briefly, identification was performed using bacteriological culture (BC) where a sterile glass loop was used to streak 10 μLof each sample to a quarter of a blood agar plate with 0.1% esculine, corresponding to one plate per cow. Then followed aerobic incubation at 37°C for up to 48 h. A positive sample was identified if there was minimum 1 colony. Pathogens were identified by morphology under a microscope. Final confirmation was carried out using MALDI-TOF MS using Biotyper Version V3.3.1.0 (DB 5989, Bruker Daltonik, 2012). Contaminated samples were determined if there were more than three different pathogens present. The farmers did not receive information on the results and did therefore not change their management procedures in response to any findings.

The statistical analysis was conducted in R 3.6.3 (17), and figures were generated using the package ggplot2 (18). We obtained data on the dates for dry-off per cow during the study period from the Danish cattle database. We estimated quarter-level prevalence and incidence per sampling interval for both herds. We corrected missing data points to increase data size by applying the following rules (15). A missing value between two positive samples was corrected to positive; a missing value between two negative samples was corrected to negative; and missing values between any negative and positive samples were corrected to negative. Additionally, if a missing value followed two positive samples, we corrected it to be positive (15). We considered infections to be recovered once a negative sample occurred after a positive sample and no positive sample followed that. The recovered animals were either spontaneously recovered or treated with antimicrobials. Data from cows after treatment was excluded from this study.

Poisson regression is a commonly applied method for estimating transmission rates in both epidemic and endemic scenarios (8, 19). We here wanted to assess and compare the transmission rates of Corynebacterium spp. using three different methods, namely Poisson regression and two other methods (20). One advantage of these methods is flexibility because longitudinal field studies can include many introductions and exits of animals between samplings. These methods are less sensitive to these effects than the Poisson regression. Furthermore, Poisson regression and method 1 are sensitive to the time between sampling intervals because they do not consider reinfections between sampling. Method 2 takes this into account, correcting for the possibility of multiple infections between samplings (20).

To estimate the transmission rate of Corynebacterium spp. with Poisson regression, we used a model with a log-link, assuming that each newly infected quarter is infected halfway between two sampling points, as described in Zadoks et al. (19):

where log(IN)^ Indicates the expected log number of newly infected quarters at each sampling interval, β is the estimated transmission rate of Corynebacterium spp. S_int and I_int are the number susceptible and infectious quarter days at risk, respectively. N_int is the total number of susceptible and infectious quarter days at risk. We used the term log(S_int · I_int/N_int) as an offset in the regression analysis, meaning that this model reflects contagious transmission directly between cows. Like in Dalen et al. (8), we also constructed a similar model with log(S_int) as the offset, they mimic a fixed rate of infection, independent of the proportion of infected animals present at the sampling points. This model reflects environmental transmission where pathogens come from a reservoir in the environment. We then compared these two models using the Akaike Information Criterion (AIC), to investigate which model fits the data best. We carried out these steps for both herds.

Additionally, we estimated the corresponding recovery rate, α, using the following equation [as described in Dalen et al. (8)]:

where log(Cm)^ accounts for the expected log number of recovered quarters per monthly sampling interval, α is the estimated recovery rate, and I_int is the number of infectious quarter days at risk per monthly interval. We calculated Wald 95% confidence limits of the transmission and recovery rate using the confint function in R.

Kirkeby et al. (20) proposed two methods for estimating the transmission rate. Both ways assume that the system is in equilibrium. These two methods do not require calculating the number of quarter days at risk for each sampling interval, but only the numbers of animals present in each state (I_N, I, S, and N) at each sampling time are required to estimate the transmission rate. The new methods yield an estimate of the transmission rate per sampling interval compared to the Poisson regression, which provides one combined transmission rate (20). However, to inspect the robustness of the Poisson regression, we performed this analysis for each sampling interval, obtaining one estimate of the transmission and recovery rates per interval.

For method 1, we estimated the transmission rate per sampling interval using the following equation:

Where β is the transmission parameter, I_N is the number of new infections per sampling interval, I is the number of infected quarters, T is the length of the sampling interval, S is the number of susceptible quarters and N is the total number of quarters after the sampling interval T. The corresponding recovery rate (α) was calculated in a similar way:

with the same parameters as described in (3). Method 2 is given by the following equation:

Where the abbreviations are the same as described for method 1. However, method 2 allows for multiple infections and spontaneous recovery between sampling intervals Kirkeby et al. (20). The corresponding recovery rate (α) is defined by:

with the same parameters as described above. Confidence limits for the transmission and recovery rates per sampling interval were derived by bootstrapping for these two methods. We sampled from all quarters 1,000 times with replacement. When sampling quarters, all data points over time for each selected quarter where used. By this, we obtained 1,000 estimates of the transmission and recovery rates for each sampling interval. We then calculated the 95% confidence intervals from the bootstrapped estimates of the transmission and recovery rates.

We further evaluated the epidemiological fitness of Corynebacterium spp. in a completely susceptible population by the basic reproduction number, R₀, given by [adapted from May (21)]:

where β is the transmission rate, and τ is the duration of infection and the inverse recovery rate (α). Here, we used the 1,000 estimates of the transmission and recovery rates from the bootstrapping and took the mean of these estimates per method and estimated the R₀. The 95% confidence limits were also calculated from these estimates.

In herd 1, 4,665 quarter milk samples were collected from 214 lactating cows (Table 1). Detailed description can be found in Supplementary Tables 1, 2. We corrected 126 and 141 missing records to positive in herd 1 and 2, respectively, as described above. Likewise, we corrected 202 and 567 missing records to negative in herd 1 and 2, respectively. The mean quarter-level prevalence was 24% (range = 12.1–38.6) (Figure 1A). During the 12-mo study period, 320 new infections with Corynebacterium spp. were detected. The mean quarter-level incidence was estimated to be 8% (range = 0.8-20.5) (Figure 1B). The mean number of recovered quarters was estimated to be 34.2 per sampling interval. In herd 2, 10,381 milk samples were collected from 443 lactating cows. The mean quarter level prevalence was estimated to be 11.7% (range 5.9–18) (Figure 1C). In total, 508 new infections with Corynebacterium spp. were found during the study period. This corresponds to a mean quarter level incidence of 6.5% (range 2.1–17.4) (Figure 1D). The mean number of recovered (cured) quarter infections was 51.6 per sampling interval.

For both herds, we compared the fit of a model where transmission was dependent on the proportion of infected animals in the herd at the sampling points (Equation 1) with the model describing a fixed transmission rate, as described above. In both herds, the former model had lower AIC (215.4) than the latter (261.8), suggesting a contagious component in the transmission of Corynebacterium spp., e.g., via the cluster during milking or other places in the herd.

The estimated rates are presented in Table 2. In herd 1, the mean estimated transmission rate was 0.016 cases/quarter day, and the standard estimated recovery rate was 0.012 cases/quarter day. In general, there was reasonably slight variation between the estimates using the different methods. In herd 2, the mean transmission rate was estimated to 0.018 cases/quarter day, and the mean recovery rate was calculated to be 0.016 cases/quarter day (Table 2).

When visually inspecting the estimated transmission and recovery rates per sampling interval for herd 1, all three methods followed similar transmission patterns (Figure 2). Generally, both transmission and recovery rates peaked in April, August, and November. However, the temporal way in recovery rate using Poisson regression appeared different from method one and method 2 with a shift in the temporal peaks (Figure 3). In herd 2, we observed that all three methods followed similar transmission and recovery patterns, except for the recovery rate estimated with Poisson regression that differed in the temporal peaks, identical to the findings for Herd 1. For all three methods, we observed that the transmission rate varied during the 12-mo sampling period, with an initial decrease, then a plateau from July to October when it decreased again before increasing again in November onwards before a final dip in January. The recovery rate followed the same course for method one and method 2 (Figures 3E,F). The estimates using Poisson regression were different from the two other methods (Figure 3).

The estimated basic reproduction number, R₀, for herd 1 and 2 for the different methods are presented in Table 2. In herd 1, it was estimated to be 1.13 (95% CI: 1.02–1.23) using the Poisson regression method (Table 2). This indicates that each infectious quarter would infect 1.13 other quarters on average during its infectious period in a completely susceptible population. For method 1 and 2, the estimated R₀ was 1.34 (95% CI: 1.30–1. 38) and 1.34 (95% CI: 1.30–1.38), respectively. The estimated mean R₀ of the three methods was 1.27. In herd 2, the R0 was 1.03 (95% CI: 0.94–1.11) using the Poisson regression method (Table 2). The R₀ for methods 1 and 2 were estimated to 1.14 (95% CI: 1.130–1.16), and the mean R₀ of the three methods was 1.10 (Table 2).