LONGITUDINAL AND CROSS-SECTIONAL STUDIES TO EVALUATE CHANGES IN COW MILK MICROBIOTA OVER THE LACTATION STAGES

To clarify the changes in milk microbiota by lactation stage, i.e., d ays in milk (DIM) for Fresh of 0 to 21, Early (DIM of 22 to 80), Middle (DIM of 81 to 200), and Late (DIM 201 or more) lactation stages in dairy cows, we performed longitudinal (12 sampling at each stage, total 48 samples) and cross-sectional sampling (Fresh [n = 7], Early [n = 12], Middle [n = 11], and Late [n = 11] stages, total of 41 samples) to comprehensively analyze the microbiota in milk samples by 16S rRNA amplicon sequencing. Although the relative abundance of bacteria of the phylum Actinobacteria increased signi ﬁ cantly in the Late lactation stage in the longitudinal study, no signi ﬁ cant changes were observed in the cross-sectional study. While no changes were observed in the milk microbiota during the lactation stages, environmental factors appeared to have a comparatively larger impact than interindividual diversity on the composition of the milk microbiota. Furthermore, the ﬁ ndings illustrated the importance of selecting appropriate study designs to clarify changes in milk microbiota throughout the lactation period. The ﬁ ndings obtained in this study not only provide useful information for interpreting previous research results, but also provide knowledge that could be helpful when designing new studies.


INTRODUCTION
Scientifi c advances have led to the development of bacterial DNA metagenomic pyrosequencing technology, which has been applied to determining the composition of microbiota in bovine milk [1,2].For example, Shinozuka et al. [3] employed such methods to describe interindividual diversity in the composition of the milk microbiota of healthy cows.The composition of milk microbiota has also been shown to be affected by seasonal [4] and environmental factors, such as bedding and airborne dust [5].In addition to annual variations in microbial composition [6], the impact of maternal nutrition on changes in the composition of milk microbiota has also been reported [7].Phase feeding, which involves nutritional design and feeding according to the lactation stage, is widely used.Differences in maternal nutritional needs at different lactation stages may also induce changes in the composition of milk microbiota.However, to the best of our knowledge, studies on the relationships between lactation stage and cow milk microbiota have been conducted only for the early stage of lactation [8], and no studies on milk microbiota of Holstein cows at any lactation stage have been conducted to date.
It is considered that elucidating the physiological changes in milk microbiota during lactation will be helpful not only for interpreting the fi ndings of previous studies, but also for the planning of further studies on milk microbiota.In the present study, we performed longitudinal and cross-sectional studies to clarify the changes in milk microbiota throughout the four lactation stages (Fresh, Early, Middle, and Late stages) in healthy cows with phase feeding.

MATERIALS AND METHODS
It has been reported that postpartum milk microbiota differs between primiparous and multiparous cows [9], the cows included in the present study were in their second or third lactation, clinically healthy, and had no history of antibiotic treatment in the previous 3 months.Lactation stages were categorized based on the days in milk (DIM), which were divided into four lactation-stage groups as follows: Fresh stage, DIM of 0 to 21; Early stage, DIM of 22 to 80; Middle stage, DIM of 81 to 200; and Late stage, DIM of 201 or more.Experiment 1 was carried out using a herd maintained by Hiroshima University, and milk samples from each mammary quarter of cows (n = 3) were collected longitudinally on four occasions during 2022 (i.e., March, May, September, and December) to give a total of 48 samples.Experiment 2 was carried out using a herd maintained by the Shizuoka Prefectural Animal Industry Research Institute.Milk samples from each mammary quarter of cows were collected crosssectionally at each lactation stage (Fresh stage: n = 7; Early stage: n = 12; Middle stage: n = 11; and Late stage: n = 11) in June of 2022 to give a total of 41 samples.Herds on both farms were housed in a freestall barn with rubber mattresses and there was no history of mastitis vaccination, whi ch contains inactivated (killed) bacteria Escherichia coli and Staphylococcus aureus, the only vaccine available in Japan.Vaccination history other than mastitis was unknown.In both experiments, personnel at the dairy farm cleaned and disinfected each teat-end of the mammary gland using a cotton gauze pad moistened with 70% ethanol while wearing gloves.After four streams of milk were stripped from each quarter, milk samples were collected aseptically into sterile 15 mL conical vials by the personnel of the dairy farm.The collected milk samples were immediately stored at 4°C, transported to the laboratory, and immediately subjected to various tests.
As screening tests to exclude milk samples from cows that might have mastitis, the somatic cell count (SCC) and the activity of lysosomal N-acetyl-β-D-glucosaminidase (NAGase) were determined for each milk sample.The SCC was determined using an electronic DeLaval cell counter (DeLaval International AB, Tumba, Sweden) according to the manufacturer's instructions.The NAGase activity was determined using the β-N-Acetylglucosaminidase Assay Kit (Sigma-Aldrich Co., LLC, St. Louis, MO, USA) according to the manufacturer's instructions.The thresholds for the screening tests were a SCC of 200,000 cells/mL or more, or an NAGase activity level of 10 nmol/ min/mL or higher.After the screening tests, 36 and 35 samples remained for further analyses of milk microbiota in Experiments 1 and 2, respectively.To analyze the milk microbiota, we performed 16S rRNA gene full-length amplicon sequence analysis using a MinION™ nanopore sequencing device (Oxford Nanopore Technologies, Oxford, UK) according to the method of Hayashi et al. [10].Briefl y, genomic DNA was fi rst isolated using an ISOSPIN Fecal DNA kit (Nippon Gene, Tokyo, Japan).Milk samples were subjected to mechanical disruption by bead-beating, and DNA was isolated using silica membrane spin columns included in the kit.Next, the V1-V9 region of the 16S rRNA gene was amplifi ed from genomic DNA by polymerase chain reaction (PCR) spiked with DNA from Haloarcula japonica as an internal control using KAPA HiFi HotStart ReadyMix (Nippon Genetics, Tokyo, Japan).Finally, the PCR product was barcoded using a PCR Barcoding Kit (SQK-RBK004, Oxford Nanopore Technologies) and purifi ed using an AMPure® XP kit (Beckman Coulter, Brea, CA, USA).The amplifi ed sequences were used to create a DNA library that was loaded onto a version R9 Spot-on Flow Cell (FLO-MIN106D, Oxford Nanopore Technologies) and sequenced on a MinION™ Mk1C device (Oxford Nanopore Technologies).MINKNO W software ver.21.11.6 (Oxford Nanopore Technologies) was used for base-calling, data acquisition and the production of FASTQ fi les.The FASTQ fi les were trimmed (the fi rst 50 nucleotides of all reads were trimmed) and fi ltered (reads with a minimum average read quality score less than 10 were removed, and sequences shorter than 500 nucleotides were removed) using Nanofi lt software (Oxford Nanopore Technologies) [11].After trimming and size selection, an average of 4,123 reads per sample were retained for bacterial identifi cation.For each read, a minimap2 search was performed using 5,850 representative bacterial genome sequences stored in the GenomeSync database [12].The taxa were assigned based on the National Center for Biotechnology Information taxonomy database [13].Low-abundance taxa (fewer than 0.01% of the total reads) were discarded from the analysis.From these results, the number of taxonomic groups (Richness) and the Shannon-Wiener index (Shannon) were calculated as alpha-diversity metrics.After confi rming the inability to achieve normalization by preprocessing for both sets of data, we performed the Friedman test (a non-parametric test for 3 or more paired groups) in Experiment 1 and the Kruskal-Wallis test (a non-parametric test for 3 or more unpaired groups) in Experiment 2, respectively.These statistical analyses were performed with the R software package (ver.4.1.3,R Core Team, 2022) [14].The analysis of similarities (ANOSIM) test was used as a beta diversity metric to evaluate the differences between the lactation stages based on the Bray-Curtis distance measure [15] implemented in the PAST software package (ver.4.03) [16].A p-value of <0.05 was considered to indicate a statistically signifi cant difference.The data for specifi c genus clusters were obtained, and bacterial taxa that were signifi cantly enriched in sample groups were extracted by linear discriminant analysis effect size (LEfSe) analysis using the Galaxy online interface [17].Parameters set for the LEfSe analysis included an alpha value of <0.05 for the factorial Kruskal-Wallis test, and a logarithmic linear discriminant analysis score threshold of <2.0 for identifying discriminative features.

Ethics
In this study, research related to animal use complied with all the relevant national regulations and institutional policies related to the care and use of animals.Experiment 1 was approved by the ethics committee of Hiroshima University (No. E19-3).In Experiment 2, the protocol for milk sample collection from lactating dairy cows was conducted in accordance with Azabu University Animal Experimentation Committee guidelines, and was approved by the ethics committee at that university (No. 200803-1).

RESULTS
After screening to exclude milk samples from cows that might have mastitis, 36 samples (n = 9 for each lactation stage) of the 48 samples in Experiment 1, and 35 samples (Fresh stage: n = 7; Early stage: n = 10; Middle stage: n = 8; and Late stage: n = 10) of the 41 samples in Experiment 2 were selected for further analyses.In Experiment 1, no signifi cant differences in Richness were observed among the lactation stages.In contrast, Shannon was signifi cantly higher in the Late stage than in the Fresh and Middle stages (Figure 1A).The mean relative abundance of milk microbiota components differed signifi cantly between the Fresh and Middle stages (p = 0.0153), the Fresh and Late stages (p = 0.0018), the Early and Middle stages (p = 0.0009), and the Middle and Late stages (p = 0.0002).Interestingly, the proportion of bacteria of phylum Actinobacteria increased signifi cantly, while those in phylum Euryarchaeota decreased signifi cantly as the lactation stages progressed (Figure 2A).In the Early stage, members of the phylum Euryarchaeota (Haloarcula) were signifi cantly more abundant, in the Middle stage, members of the phylum Actinobacteria (Corynebacteriaceae, Kocuria, and Micrococcaceae), phylum Firmicutes (Bacillaceae), and phylum Proteobacteria (Shigella) were signifi cantly more abundant, and in the Late stage, members of the phylum Actinobacteria (Microbacteriaceae and Micrococcales) and the phylum Proteobacteria (Escherichia, Enterobacteriaceae, and Pseudomonadales) were signifi cantly higher (Figure 3).In Experiment 2, no differences were observed in Richness or Shannon among any of the groups (Figure 1B).The mean relative abundance of bacteria in the phylum Firmicutes increased, while bacteria in the phylum Actinobacteria decreased as the lactation stages progressed; however, the differences were not statistically signifi cant (Figure 2B).

DISCUSSION
In Experiment 1, both the Shannon index, a measure of the evenness of bacterial species in a community, and the relative abundance of bacteria in the phylum Actinobacteria increased as the lactation stages progressed, while the relative abundance of members of the phylum Euryarchaeota decreased.These fi ndings indicate that the components of the milk microbiota change at different lactation stages.Conversely, in Experiment 2, no relationship was observed between the lactation stages and the composition of milk microbiota.In this study, DNA extracted from 2.8 × 10 3 cells of Haloarcula japonica in the phylum Euryarchaeota was added as an internal control in PCR amplifi cation.Consequently, the relative abundance of members of the phylum Euryarchaeota could be estimated by comparison with the total amount of bacteria contained in the sample, including the internal control.The reason the percentage of members of the phylum Euryarchaeota was high in the Fresh and Early stages in Experiment 1 may be because the total amount of bacteria DNA contained in the samples was small, making them appear relatively high.However, even so, the results from the two experiments were contradictory.
One possible reason for the contradictory results may be related to the timing of sample collection in the experiments.Experiment 1 was a longitudinal study that was conducted to eliminate the infl uence of interindividual diversity, and the results would have included the infl uence of the environment on the milk microbiota.In Experiment 1, samples from the Early stage were collected before summer (in May), and those from the Middle and Late stages were collected after summer (September and December, respectively), and signifi cant differences were seen in the milk microbiota compositions between the Fresh/Early and the Middle/Late stages, and the Shannon index was also signifi cantly higher in the Late stage.Celano et al. [4] reported that microbiological variability can arise due to seasonal factors.They observed that the proportions of members of the phyla Actinobacteria and Firmicutes increased in summer as a result of different microbial growth conditions.It is therefore possible that the increase observed in the proportion of members of the same two phyla in the Middle stage (September) in this study could also be attributed to seasonal factors.In addition, our results were consistent with those of a previous study [5], which found that the relative abundance of both phyla, especially the families Micrococcaceae in phylum Firmicutes and Bacillaceae in phylum Actinobacteria, increased in summer compared to winter.Although Nguyen et al. [5] reported that bacteria from several major bacterial families (e.g., Staphylococcaceae, Moraxellaceae, Ruminococcaceae, and Bacteroidaceae) that are commonly found in bedding and airborne dust can potentially infl uence the components of milk microbiota, our study did not examine these taxa, and so no such relationships could be clarifi ed.It is thus possible that the results of Experiment 1 may have been affected by season and/or environmental factors.Further, the meta-genomic analysis conducted in this study was performed using 16S rRNA gene sequences, which cannot be used to determine the absolute number of each bacterial species.Consequently, to more accurately evaluate the origin of bacteria that may have affected the composition of milk microbiota, it will be necessary to investigate the absolute numbers of bacterial taxa using other methods, such as quantitative PCR [18].
In Experiment 2, which was designed to negate the potential infl uence of environmental factors, lactation stage was observed to have no effect on the composition of milk microbiota.By analyzing milk collected from cows within two months after parturition, Zhu et al. [8] reported that the components of milk microbiota change during the early lactation stage.In this study, the whole milk stage was divided into four stages, so there may not have been any signifi cant differences among stages.In addition, although Experiment 2 was a cross-sectional study that was not affected by environmental factors, it has been reported that marked differences can exist in the milk microbiota of milk samples from different cows [3].Thus, given the possibility of such interindividual differences, we could not conclude whether there was a relationship between the lactation stages and the composition of milk microbiota.
To clarify the infl uence of lactation stage on milk microbiota, we conducted a longitudinal study to negate the infl uence of interindividual diversity and a crosssectional study to negate the effect of environmental infl uences on the components of milk microbiota.The results showed signifi cant differences in the composition of microbiota in the longitudinal study, but not in the cross-sectional study.Although it was not possible to conclude whether the lactation stage affects the composition of milk microbiota because the results of the two experiments were different, the fi ndings showed that environmental factors play a relatively greater role than interindividual diversity in affecting the composition of milk microbiota.Although multiple factors are thought to affect the composition of milk microbiota, this study specifi cally elucidated the extent of infl uence attributable to two key factors.The fi ndings therefore not only provide useful information for interpreting previous research results, but they also provide knowledge that could be helpful when designing future studies.manuscript.TK and KK supervised the entire study.All of the authors provided critical feedback, assisted with the research, and prepared the fi nal manuscript.

Declaration of confl icting interests
The author(s) declared no potential confl icts of interest with respect to the research, authorship, and/or publication of this article.

Informed consent
Informed consent was obtained from clients for client-owned animals included in this study.The owners understood the experimental procedure and agreed that the results related to this investigation could be published in Scientifi c Journal Acta Veterinaria-Beograd.

Figure 1 .
Figure 1.Changes in the alpha-diversity index (Richness, Shannon) by lactation stage in the longitudinal study (A) and the cross-sectional study (B).The lactation stages were categorized into four lactation-stage groups according to the days in milk (DIM); Fresh stage (DIM of 0 to 21), Early stage (DIM of 22 to 80), Middle stage (DIM of 81 to 200), and Late stage (DIM of 201 or more).The index was compared among groups by the Friedman test in the longitudinal study (A), and by the Kruskal-Wallis test in the crosssectional study (B).A p-value of <0.05 was considered to indicate statistical signifi cance.

Figure 2 .
Figure 2. Changes in the mean relative abundance of the milk microbiota components by lactation stage in the longitudinal study (A) and the cross-sectional study (B).The graphs show the proportion of each bacterial phylum inferred from polymerase chain reaction amplifi cation and pyrosequencing of 16S rRNA amplicons.A and B: Bacterial composition of milk microbiota over time according to the days in milk (DIM).The lactation stages were categorized into four lactation-stage groups based on DIM; Fresh stage (DIM of 0 to 21), Early stage (DIM of 22 to 80), Middle stage (DIM of 81 to 200), and Late stage (DIM of 201 or more).Milk microbiota were compared among groups by an analysis of similarities with Bonferroni corrections.A p-value of <0.05 was considered to indicate statistical signifi cance, and is indicated by asterisks.

Figure 3 .
Figure 3. Taxonomic cladogram generated from a linear discriminant analysis effect size (LEfSe) analysis in the longitudinal study.Cladogram showing the relative proportions of the bacterial species of the milk microbiota that increased signifi cantly among lactation stages estimated by the Kruskal-Wallis test.Red symbols indicate Early-stage samples, blue symbols indicate Middle-stage samples, and green symbols indicate Late-stage samples.