Molecular Epidemiology of Rotavirus A in Calves: Evolutionary Analysis of a Bovine G8P[11] Strain and Spatio-Temporal Dynamics of G6 Lineages in the Americas

Rotavirus A (RVA) causes diarrhea in calves and frequently possesses the G6 and P[5]/P[11] genotypes, whereas G8 is less common. We aimed to compare RVA infections and G/P genotypes in beef and dairy calves from major livestock regions of Argentina, elucidate the evolutionary origin of a G8 strain and analyze the G8 lineages, infer the phylogenetic relationship of RVA field strains, and investigate the evolution and spatio-temporal dynamics of the main G6 lineages in American countries. Fecal samples (n = 422) from diarrheic (beef, 104; dairy, 137) and non-diarrheic (beef, 78; dairy, 103) calves were analyzed by ELISA and semi-nested multiplex RT-PCR. Sequencing, phylogenetic, phylodynamic, and phylogeographic analyses were performed. RVA infections were more frequent in beef (22.0%) than in dairy (14.2%) calves. Prevalent genotypes and G6 lineages were G6(IV)P[5] in beef (90.9%) and G6(III)P[11] (41.2%) or mixed genotypes (23.5%) in dairy calves. The only G8 strain was phylogenetically related to bovine and artiodactyl bovine-like strains. Re-analyses inside the G8 genotype identified G8(I) to G8(VIII) lineages. Of all G6 strains characterized, the G6(IV)P[5](I) strains from “Cuenca del Salado” (Argentina) and Uruguay clustered together. According to farm location, a clustering pattern for G6(IV)P[5] strains of beef farms was observed. Both G6 lineage strains together revealed an evolutionary rate of 1.24 × 10−3 substitutions/site/year, and the time to the most recent common ancestor was dated in 1853. The most probable ancestral locations were Argentina in 1981 for G6(III) strains and the USA in 1940 for G6(IV) strains. The highest migration rates for both G6 lineages together were from Argentina to Brazil and Uruguay. Altogether, the epidemiology, genetic diversity, and phylogeny of RVA in calves can differ according to the production system and farm location. We provide novel knowledge about the evolutionary origin of a bovine G8P[11] strain. Finally, bovine G6 strains from American countries would have originated in the USA nearly a century before its first description.


Introduction
Neonatal calf diarrhea (NCD) is a leading cause of disease and death in dairy and beef calves [1].It generates important economic losses to the livestock industry due to poor growth, death, and expenses in treatments [1,2].Also, in cow-calf beef farms, treating diarrheic calves is labor-intensive and physically risky [1].The etiology of NCD is complex and involves enteric pathogens (virus, bacteria, and protozoa) along with predisposing factors such as the environment, immunity, and nutrition [2][3][4].
Rotaviruses (RVs) belong to the family Sedoreoviridae and genus Rotavirus.Currently, the International Committee on Taxonomy of Viruses (ICTV) recognizes nine species (groups) that are designated RVA to RVD and RVF to RVJ [5].The serological crossreactivity and amino acid (aa) sequence identity of VP6 as well as the host ranges allow for us to distinguish these species [6,7].RVA is not only the main pathogen associated with NCD [2] but also the most clinically relevant species infecting other ruminants [3].Previous studies in different countries report a variable frequency of RVA infections in diarrheic calves, ranging from 27% to 79.9% [8][9][10][11][12][13][14].Aside from other factors, virus detection methods and management practices adopted in cattle farms may explain this variability.Although subclinical infections contribute to environmental contamination [3], few studies have compared RVA infections in diarrheic and non-diarrheic (subclinical) calves from beef and dairy farms [8,13,14].
The whole genome classification system of RVA strains is based on nt percentage identity cutoff values, which allow for us to define the genotypes for all the 11 gene segments.Thus, the schematic nomenclature Gx-P[x]-Ix-Rx-Cx-Mx-Ax-Nx-Tx-Ex-Hx (x = Arabic numbers) denotes the genotypes for the VP7-VP4-VP6-VP1-VP2-VP3-NSP1-NSP2-NSP3-NSP4-NSP5/6 coding genes, respectively [27].This classification enables us to elucidate the most likely origin and trace the evolutionary patterns of animal and human strains [28].Also, this approach is very important because interspecies transmission and reassortment events are major ways to increase the genome diversity of RVA strains [3,29].
Continuous knowledge of the epidemiology and genetic diversity of RVA is essential to implement preventive measures [16].The aims of this study were as follows: (a) to compare the frequency of RVA infections and the distribution of G/P genotypes in diarrheic and non-diarrheic calves from beef and dairy farms in major livestock regions of Argentina; (b) to elucidate the evolutionary origin of an unusual G8 strain and analyze the G8 lineages; (c) to infer the phylogenetic relationship of RVA study strains and assess the influence of the production system and the location of farms on their relationship; (d) investigate the evolution and spatio-temporal distribution of the major G6 lineages in American countries.The resulting data of this study show differential trends when comparing the infection frequency, G/P genotypes, genetic diversity, and phylogeny of RVA strains between calves from beef and dairy farms.Our research provides novel knowledge about the evolutionary origin of bovine G8 strains.Finally, in American countries, bovine G6 strains would have originated in the USA nearly a century before its first description in diarrheic calves from this country.

Farms, Calves, and Sample Collection
From September 2007 to December 2010, a total of 32 farms were selected by convenience (non-probability method) [46], because veterinary practitioners requested to detect enteric pathogens during NCD outbreaks.The farms were located in Buenos Aires province, where the largest cattle population of Argentina is reared (14,883,528 out of 40,023,083; ~37% of all cattle) [47].Beef farms (n = 14) were located in eight districts inside (7/8) or just outside (1/8) the beef production region of "Cuenca del Salado and Depresion de Laprida" (Supplementary Table S1).This region accounts for 48% of the cattle population in Buenos Aires province and is the most important production region in Argentina [48].Geographically, this region can be split into "Cuenca del Salado" and "Depresion de Laprida" [49].Dairy farms (n = 18) were located in 10 districts and three dairy regions, namely "Cuenca Mar y Sierras" (n = 14), "Cuenca Oeste" (n = 3), and "Cuenca Abasto Norte" (n = 1) [50].Three districts included both beef and dairy farms (Supplementary Table S1).

RNA Extraction and G/P Genotyping
To conduct the multiple based genotyping of RVA field strains in Ag-ELISA positive samples, dsRNA was extracted from fecal suspensions with TRIzol (Invitrogen, Carlsbad, One dairy farm, designated DR, was visited in 2009, 2010, and 2011 due to severe NCD outbreaks.On this farm, the fecal samples collected in 2009 were included for statistical, genotype, and phylogenetic analyses, whereas the samples collected in 2010 and 2011 were included for genotype comparisons and phylogenetic analysis.

RNA Extraction and G/P Genotyping
To conduct the multiple based genotyping of RVA field strains in Ag-ELISA positive samples, dsRNA was extracted from fecal suspensions with TRIzol (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions.The dsRNA pellet was air dried and resuspended in 50 µL of nuclease-free water (Biodynamics, Buenos Aires, Argentina).

Sequencing of VP7 and VP8*
Ag-ELISA positive samples from diarrheic calves were selected according to the following criteria: single G and/or P genotype in the sample, production system (beef or dairy), and farm location in production regions.The RT-PCR products of VP7 (1062 bp) and VP4 (VP8* subunit, 877 bp) genes were purified using illustra ExoProStar (GE Healthcare, Little Chalfont, UK).Then, sequencing by the dideoxy chain termination method was conducted bi-directionally with generic primers [11].This analysis was performed on an automated DNA sequencer ABI 3500 (Applied Biosystems, Foster city, CA, USA) at the Genomic Unit of the Biotechnology Institute (CICVyA, INTA, Hurlingham, Argentina).

Sequencing of RVA Gene Segments
For the full genome sequencing of a selected RVA strain, dsRNA was extracted from the fecal suspensions of an Ag-ELISA positive sample (4385) using the QIAamp Viral RNA Mini Kit (Qiagen/Westburg, Leusden, The Netherlands).The dsRNA gene segments were reverse transcribed and amplified using the Qiagen OneStep RT-PCR kit (Qiagen/Westburg, Leusden, The Netherlands), with primers and cycling conditions described previously [53,54].Briefly, the stages of RT (50  C/3 min for VP7, VP6, and NSP1-NPS5 genes.Both cycle conditions were followed by a final extension stage (72 • C/10 min).The RT-PCR products were purified with the MSB Spin PCRapace kit (Invitek, Berlin, Germany) and sequenced by using the BigDye Terminator Cycle Sequencing Reaction kit (Applied Biosystems, Foster city, CA, USA) on an automated sequencer ABI PRISM 3130 (Applied Biosystems).Sequencing reactions were conducted bi-directionally with the primers used for each RT-PCR (primers available upon request).Primer walking sequencing with internal primers was performed to cover the ORF sequences of the longest genes (VP1-VP4).This study was conducted at Dr. Jelle Matthijnssens laboratory at the Rega Institute for Medical Research, University of Leuven, Belgium.

Sequence Analysis, Genotype Assignment, and Percentage Identity
Sequencing electropherograms were inspected and corrected with the ChromasPro V1.5 software (Technelysium Pvt. Ltd., 2009).To obtain the ORF sequences of all the 11 gene segments, the nt sequences generated with forward, reverse, and internal primers were aligned and assembled using BioEdit v7.0.2.5 [55].These sequences were evaluated to accomplish the guidelines proposed by the Rotavirus Classification Working Group (RCWG) [56].The genotype for each gene segment was assigned using the Rotavirus A Genotype Determination Tool, which was available in the Virus Pathogen Database and Analysis Resource (ViPR) [57].Currently, this tool has migrated to the Bacterial and Viral Bioinformatics Resource Center.
The Nucleotide BLAST tool (https://blast.ncbi.nlm.nih.gov/Blast.cgi;accessed on 20 March 2023) was run to retrieve similar sequences deposited in the GenBank database.Particularly, the nt sequences of G8 and/or P [11] strains were collected if their backbone genotype constellations were reminiscent to that of the study strain.Multiple sequence alignment files were prepared using ClustalW with default settings, as integrated in the MEGA 7.0 software [58].These sequence datasets comprised complete and partial ORF sequences of each gene segment.Identity matrices, prepared with the aligned nt and deduced aa sequences, were calculated using the p-distance algorithm [30,59] and the pairwise distance analysis [28] implemented in MEGA 7.0.
To infer the phylogeny of G6, P [5], and P [11] RVA strains according to the production system and the location of farms in Buenos Aires province, additional ORF sequences corresponding to VP7/VP8* (nine strains), VP7 (11 strains), and VP8* (12 strains) were included from previous studies in Argentina [41,44].Subsequently, phylogenetic reconstruction was carried out using the Kimura 2-parameter substitution model and the Neighbor-Joining (NJ) method as implemented in MEGA 7.0 [28,41].

Determination of G8 Lineages
The VP7-G8 genotype dataset included 183 complete and partial (≥73% completeness) ORF sequences.This dataset, which was designated RVA-G8(lineages) (see Supplementary Materials), did not include sequences of RVA strains with 100% identity and obtained from the same host species and/or country.The phylogenetic tree was inferred using the Kimura 2-parameter substitution model and the NJ method [28,41].Thereafter, a cutoff value to best discriminate G8 lineages was estimated.Briefly, a matrix of genetic distance values was prepared using the pairwise distance algorithm and the p-distance model available in MEGA 7.0 [58].The pairwise distance frequency graph was created by plotting the pairwise distance values, expressed as percentages, in the abscissa (x-axis) and their frequencies in the ordinate (y-axis).The genetic distance values between strains of the same lineage were designated "intra-lineage distances", whereas the genetic distance values between strains of different lineages were designated "inter-lineage distances".The cutoff value enables to differentiate between the intra-and inter-lineage genetic distances [28,41].The lineage numbers were assigned in accordance with previous papers [38,40].

Phylodynamic and Phylogeographic Analyses of G6 RVA Strains
Three VP7 sequence datasets [RVA-G6, RVA-G6(III), and RVA-G6(IV)] (see Supplementary Materials) were prepared for evolutionary analyses of bovine G6 RVA strains from American countries.All nt sequences were confirmed as non-recombinant using the Recombination Detection Program (RDP) v4.100 [63].The presence of substitution saturation was not detected by the Xia's test implemented in the DAMBE software v7.2.102 [64].The phylogenetic signal was evaluated by the likelihood mapping method available in IQ-TREE v1.6.12.The phylogenetic noise was low to rather low for each sequence dataset (15.9%, 32.2%, and 7.5%, respectively).A positive correlation between the genetic divergence and sampling time was observed using the Root-to-Tip analysis with TempEst v1.5.3 [65] (see Supplementary Materials).
The spatio-temporal process was evaluated on time-scaled genealogies over discrete sampling locations (countries) using a Bayesian stochastic search variable selection procedure [68] and an asymmetric model.This analysis was performed in the BEAST package v1.8.4.A Bayes Factor (BF) test was used to weigh the significance of the linkages between countries, using the SpreaD3 software v0.9.6 [69].BF > 3 was considered a well-supported link.

Statistical Analyses
The statistical analysis was conducted in R v4.2.2 (R Core Team, Vienna, Austria).The data were analyzed using binary and multinomial logistic regression models that were fitted with the R packages 'survey' v4.1-1 [71] and 'svyVGAM' v1.2 [72], respectively.In binary logistic regression models, bias reduction methods were applied (r package 'brglm2 v0.9) [73].In all the analyses, the year of sampling was included as a stratifying factor and the farm was considered as a clustering factor.
The frequency of RVA infections was compared between production systems (beef or dairy) using binary logistic regression.The analysis was adjusted by the clinical condition of the calves (diarrheic or non-diarrheic).Also, the effect measure modification by clinical condition was assessed in the multiplicative scale by including in the model the product term between production system and clinical condition.The frequency of RVA infections was compared between clinical conditions using binary logistic regression and adjusting the analysis by the production system.The effect measure modification by production system was also assessed in the multiplicative scale following the aforementioned methodology.For this analysis, only the 27 farms with RVA-infected calves were included in order to avoid positivity violations [74].Also, animals infected with bovine coronavirus (BCoV) (three dairy calves from two dairy farms) were excluded from this analysis.The latter was conducted considering that the interaction between RVA and BCoV might increase the risk of diarrhea and that this interaction could not be addressed properly.
The overall distribution of G and P genotypes or G/P combinations was compared between production systems using multinomial logistic regression.Also, separated binary logistic regression models were used to compare the proportion of each genotype between production systems.The False Discovery Rate correction [75] was applied to account for multiplicity.
A comparative analysis of genotype constellations and sequence identities was performed across the 11 gene segments of strain 4385VT_D_BA and those of other reference G8 and/or P [11] strains distributed worldwide.In this analysis, the nt and deduced aa sequences of each homologous genotype were compared using the ORF sequences (Table 3).Of note, strain 4385VT_D_BA shared 10 genotypes with the bovine strain B383 (G15P [11]) and nine genotypes with the artiodactyl bovine-like strains Rio_Negro (G8P [1]) and Chubut (G8P [14]), both from guanacos [76], and the caprine bovine-like strain 0040 (G8P [1]) [54].All of them were reported in Argentina.Taken together, the backbone genotype constellations I2-R2/R5-C2-M2-A3/A11/A13-N2-T6-E2/E12-H3 were conserved among G8 strains obtained from various animal species of the order Artiodactyla, as well as humans, believed to be examples of interspecies transmissions.Particularly, the 11 ORF sequences of strain 4385VT_D_BA were found to share high genetic similarity not only with the aforementioned strains but also with two artiodactyl bovine-like G8 strains, namely C75 (G8P [14]), collected from a vicugna in Argentina [53], and 562 (G8P [14]) and 1115 (G8P [1]), both obtained from alpacas in Peru [78] (Table 3).Phylogenetic studies were conducted to assess the most likely evolutionary origin of the G8 strain under study.

Phylogenetic Analysis of the G8P[11] RVA Strain
The ML trees constructed with the 11 ORF sequences of strain 4385VT_D_BA (details in Table 2) are shown in Figure 2a-f.For the VP7 gene, strain 4385VT_D_BA was most closely related to the caprine bovine-like strain 0040, which both segregated in a subcluster that only included Argentinean artiodactyl bovine-like strains (Figure 2a; % identity in Table 3).Regarding the VP4 gene, strain 4385VT_D_BA was most closely related to the Uruguayan bovine G10P [11] strains LVMS1837 and LVMS2625 (99% nt identity) [14], which grouped inside a large cluster of other Argentinean strains from calves [41] (Figure 2a).
For the VP6-I2 genotype, strain 4385VT_D_BA was segregated in a subcluster that contained four Uruguayan bovine P [11] strains (93% nt identity) (Figure 2b).For the VP1-R5 genotype, strain 4385VT_D_BA was most closely related with the caprine bovine-like strain 0040.Both strains also grouped together with the bovine strain B383 and the artiodactyl bovine-like strain Rio_Negro (Figure 2b; % identity in Table 3).All these strains formed a subcluster of Argentinean R5-carrying strains.
Remarkably, in the phylogeny of the VP2-C2 genotype, strain 4385VT_D_BA clustered most closely with the Peruvian alpaca strains 562 and 1115, as well as with the bovine strain B383 (Figure 2c; % identity in Table 3), rather distantly related to other known C2carrying RVA strains.For the VP3-M2 genotype, strain 4385VT_D_BA clustered between the Argentinean artiodactyl bovine-like stains Chubut and C75, in one branch, and the bovine strains B383 and Y136 and the caprine bovine-like strain 0040, in another branch (Figure 2c; % identity in Table 3).The strain Y136 (G8P [11]) was previously detected in a Brazilian calf [39].
closely related to the caprine bovine-like strain 0040, which both segregated in a subcl ter that only included Argentinean artiodactyl bovine-like strains (Figure 2a; % identity Table 3).Regarding the VP4 gene, strain 4385VT_D_BA was most closely related to Uruguayan bovine G10P [11] strains LVMS1837 and LVMS2625 (99% nt identity) [1 which grouped inside a large cluster of other Argentinean strains from calves [41]  For the VP6-I2 genotype, strain 4385VT_D_BA was segregated in a subcluste contained four Uruguayan bovine P [11] strains (93% nt identity) (Figure 2b).For the R5 genotype, strain 4385VT_D_BA was most closely related with the caprine bovin strain 0040.Both strains also grouped together with the bovine strain B383 and the dactyl bovine-like strain Rio_Negro (Figure 2b; % identity in Table 3).All these s formed a subcluster of Argentinean R5-carrying strains.
The mean intra-lineage distance values ranged from 0.4 to 7.3%, whereas the mean inter-lineage distance values ranged from 12.4 to 17.8% (Figure 3B).According to the NJ tree topology (Figure 3A) and the pairwise distance frequency graph (Figure 3C), the cutoff value that best discriminates G8 lineages was determined at 10% genetic distance (90% nt identity) (see dotted line (Figure 3C)).The strains showing distance values lower than 10% frequently belonged to the same lineage, and strains with higher distance values frequently belonged to different lineages (Figure 3C).

Phylogenetic Analyses of VP7 and VP8* among G6 and G10 Strains
Twenty-two RVA strains from diarrheic beef and dairy calves of our study were selected for phylogenetic analyses.The G/P genotypes and G6 lineages [G6(III) or G6(IV)] determined by both SnM RT-PCR assays were in complete agreement with the results obtained by the Rotavirus A Genotype Determination Tool (ViPR) and the phylogenetic analysis of VP7, respectively.
The ML trees of VP7-G6 and VP8*-P [5] sequences showed that 10 strains from diarrheic beef calves segregated in two subclusters either within the G6(IV) (Supplementary Figure S2) and P [5](I) (nine strains) lineages (Supplementary Figure S3).One subcluster mainly grouped G6(IV)P [5](I) strains from Uruguay (Castells et al., 2020) together with five strains from the subregion "Cuenca del Salado".The other subcluster grouped G6(IV)P [5](I) strains reported previously in Argentina [41] together with the five strains from the subregion "Depresion de Laprida"; these strains were phylogenetically very closely related (Supplementary Figures S2 and S3).Of note, within the G6(IV) lineage, the strain 5595DR_D_BA obtained from a diarrheic dairy calf was closely related to the North American bovine strain B641 (Supplementary Figure S2).The ML tree of VP7-G6 sequences also showed that six strains from diarrheic dairy calves were segregated in different branches of the G6(III) lineage.These branches included bovine P [11] strains from Argentina [41] and Uruguay [14] (Supplementary Figure S2).Likewise, in the ML tree of VP8*-P [11] sequences, most of these strains clustered within the P [11](VI) lineage in different branches (Supplementary Figure S4).
The ML analysis of VP7-G10 sequences revealed that three strains from diarrheic dairy calves clustered in the G10(VI) lineage.Two strains showed a very close phylogenetic relationship with bovine G10P [5] strains from Argentina [41], whereas strain 3DR_D_BA displayed a rather close relationship with bovine G10P [11] strains from Uruguay [14] (Supplementary Figure S5).Regarding the phylogenetic analysis of P [11] sequences, two strains segregated in the P [11](VI) lineage, whereas strain 3DR_D_BA showed a rather close relationship with Uruguayan strains, all of which may represent a new lineage.The only VP8*-P [11] sequence obtained from a diarrheic beef calf (strain 11LC_B _BA; G6(IV)/G10P [11]) clustered in the P [11](VI) lineage and was closely related to strain 27412VA_ D_BA (G6(III)P [11]), which was collected from a diarrheic dairy calf (Supplementary Figure S4).

Phylogeny of G6 RVA Strains According to Production Systems and Farm Location
To evaluate the phylogenetic relationship of G6(IV)P [5] and G6(III)P [11] strains according to the production system and the location of the farm within production regions, a selection of RVA strains from this and previous studies [41,44] was analyzed using the VP7 and VP8* sequences available in GenBank.Only strains from Buenos Aires province were included in the analysis.The clustering patterns are shown in NJ trees in Figure 4.
"Cuenca del Salado" grouped within subcluster-2 (99% nt identity) and subcluster-3 (98% nt identity).Also, for the VP8*-P [5] sequences, most of these strains grouped within subcluster-1 (99% nt identity).Regarding the VP7 and VP8* sequences of RVA strains belonging to the G6(III) lineage and P [11] genotype, respectively, the NJ trees showed no clustering pattern according to the farm location (Figure 4A,B).Notably, according to our records, almost all G6(IV)P [5] strains were derived from beef calves, whereas almost all G6(III)P [11] strains were derived from dairy calves; reflecting that beef and dairy farms are located in different production regions of Buenos Aires province (Figure 4C).[5] and P [11] genotypes.Phylograms were reconstructed using the Kimura 2-parameter substitution model and the neighbor-joining method (MEGA 7.0 software).The phylogram of VP7-G6 lineages was rooted using strain RVA/Cow-wt/ARG/B383/1998/G15P [11].The analysis only included similar sequences and with available information.Each strain is represented with a number and symbol; white symbols correspond to beef farms, whereas black symbols correspond to dairy farms, as follows: □ "Cuenca del Salado"; △ "Depresion de Laprida"; ■ "Mar y Sierras"; ▲ "Oeste"; • "Abasto Norte"; ▼ "Abasto Sur".Bootstrap values (1000 replicates) >70% are shown as  [5] and P [11] genotypes.Phylograms were reconstructed using the Kimura 2-parameter substitution model and the neighbor-joining method (MEGA 7.0 software).The phylogram of VP7-G6 lineages was rooted using strain RVA/Cow-wt/ARG/B383/1998/G15P [11].The analysis only included similar sequences and with available information.Each strain is represented with a number and symbol; white symbols correspond to beef farms, whereas black symbols correspond to dairy farms, as follows: "Cuenca del Salado"; "Depresion de Laprida"; "Mar y Sierras"; "Oeste"; • "Abasto Norte"; "Abasto Sur".Bootstrap values (1000 replicates) >70% are shown as branch nodes.The name of the strain, type of production system, and province of collection are indicated.Abbreviations: D/B, Dairy or Beef; BA, Buenos Aires.(C) The map of Buenos Aires province shows RVA strains according to the production system (dairy or beef) and the location of farms in districts and production regions; each strain is represented with numbers and symbols as indicated previously.The location of the farms does not represent GPS coordinates.

Phylodynamic and Phylogeographic Analyses of G6(III) and G6(IV) Lineages
In American countries, and particularly in Argentina, the evolution and possible origin of bovine RVA strains belonging to G6(III) and G6(IV) lineages remain unknown.
According to the spatio-temporal reconstruction of all G6 strains together, the most probable ancestral location for the G6(III) lineage strains was ARG in the 1980s (95% HDP, 1963-1989), while for the G6(IV) lineage strains this was the USA in the 1940s (95% HDP, 1919-1959) (Table 4; Figure 5A).Similar dates were obtained for the RVA-G6(III) and RVA-G6(IV) datasets analyzed separately (Supplementary Figure S6A,B).The spread processes of bovine G6 RVA strains were inferred.The highest migration rates (BF > 1000) for all G6 strains together were from ARG to BRA and ARG to URY (Figure 5B).This last route of spread was also decisive (BF >1000) either for G6(III) or G6(IV) lineage strains (Supplementary Figure S7).Another route of spread with strong support (100 < BF < 1000) for all G6 strains was from the USA to ARG (Figure 5B).Routes of diffusion with low support (3 < BF < 100) for the G6(IV) lineage strains were from CAN to ARG and USA to BRA.Other low-supported pathways were observed between BRA and VEN or CAN and MEX (Supplementary Figure S7).

Discussions
NCD generates important economic losses due to increased mortality, treatment costs, and poor growth [16].Epidemiological studies of enteric pathogens in many countries have demonstrated that RVA and Cryptosporidium spp.are the most frequent causes of NCD [8][9][10]12,13,80].Hence, ongoing epidemiological studies are of great importance to assess the impact of RVA infections on livestock health.In this study, RVA was detected in at least one calf in ˃80% of beef and dairy farms during NCD outbreaks.Similarly, in the UK and Spain, calves shedding RVA were found in 58% to 93% of beef and dairy farms [8,13,81].The broad distribution of this virus is due to its environmental resistance, which allows for long persistence in calf feces (up to 9 months) [82], the short period of

Discussions
NCD generates important economic losses due to increased mortality, treatment costs, and poor growth [16].Epidemiological studies of enteric pathogens in many countries have demonstrated that RVA and Cryptosporidium spp.are the most frequent causes of NCD [8][9][10]12,13,80].Hence, ongoing epidemiological studies are of great importance to assess the impact of RVA infections on livestock health.In this study, RVA was detected in at least one calf in >80% of beef and dairy farms during NCD outbreaks.Similarly, in the UK and Spain, calves shedding RVA were found in 58% to 93% of beef and dairy farms [8,13,81].The broad distribution of this virus is due to its environmental resistance, which allows for long persistence in calf feces (up to 9 months) [82], the short period of incubation, and the large amount of virus shed by diarrheic calves for several days, even after clinical recovery [83].
In this study, the epidemiology of RVA showed a different trend in beef and dairy farms.The frequency of RVA infections in both diarrheic and non-diarrheic calves was higher in beef farms.This finding supports previous studies conducted in Argentina and Brazil (beef calves, 23% and 46% vs. dairy calves, 16% and 24%) [80,86].According to these data, NCD associated with RVA infections has a greater impact on beef cattle [82].In beef herds, bull breading usually lasts 90-120 days, which results in many calvings during a short period.Thus, virus transmission between calves of similar age is enhanced by their gregarious behavior (crowding) and direct contact with feces and contaminated udders [82,86].Also, at the beginning of the calving season, calves are usually infected subclinically and become biological amplifiers [87].Conversely, in the dairy farms of our region, calves are frequently reared in individual stakes or hutches until weaning at two months of age.This practice may reduce virus transmission by direct contact.Subclinical RVA infections have been seldomly investigated [14,88].In our study, RVA was detected in non-diarrheic beef (10%) and dairy (<5%) calves.Virus shedding during the incubation period or after clinical recovery can occur in infected calves [89], but infections can remain subclinical when high titers of specific passive maternal antibodies are present in the intestinal lumen [83].
The surveillance of G/P types among bovine RVA field strains is important because it allows for the discovery of new or emerging strains [20,90], provides insights into the epidemiology of animal infections [18] which may represent interspecies transmission events [77,90,91], and encourages the development and evaluation of vaccines [16,19,25,92].Our results are in line with previous data available in five continents where G6 (56.7%),P [5] (25.9%), and P [11] (21.5%) are the most frequent genotypes infecting cattle [16].However, we observed consistent differences in the frequency and diversity of G/P genotypes according to the production system.The G6(IV)P [5] combination was the most frequent in beef calves (≈90%), while the G6(III)P [11] combination and mixed genotypes prevailed in dairy calves (≈60%).These findings agree with a previous study in Argentina [11].However, we observed that mixed genotypes and mixed RVA outbreaks were very uncommon in beef calves and farms, respectively.In contrast, mixed RVA outbreaks were very frequent in dairy farms, as reported previously [44].Remarkably, only three studies report a different distribution of G/P types according to the production system.In beef herds from Sweden and the USA, G6 (67% and 100%) was more frequent than G10 (47.5%).Conversely, in dairy herds, G10 (17.5% and 54%) and G6 + G10 (10%) prevailed over G6 (32%) [17,93].In Brazil, G6P [5] (65.5%) was most frequent in beef herds, while G10P [11] (38.4%) and G6P [11] (30.8%) prevailed in dairy herds [25].Future studies should consider that production systems and cattle breeds, among other factors, may bias the fitness of different RVA genotypes.
So far, in Argentina and other countries, it is unknown which factors govern the differential distribution of G/P genotypes in beef and dairy farms.Probably, livestock management practices, breeds, and the unusual exchange of cattle between beef and dairy farms are involved.Beef herds are managed extensively with a calving season during winter.Thus, highly virulent strains should somehow persist in the environment and infect susceptible calves the next year or circulate subclinically in adult cattle for long periods [11,44].The G6(IV)P [5] strains appear to have these features [11], and the lack of susceptible calves during a large part of the year could explain the low diversity of G/P genotypes and lineages observed in beef farms.In contrast, in dairy farms calves are usually born all year-round [44], and the failure of passive transfer of colostral antibodies is frequent [94].This situation continuously offers susceptible calves, which enhances the opportunity for infections by less virulent strains [11] or co-infections by different strains.Consequently, genome reassortments are more likely to occur in dairy calves, which together with point mutations, increase the genetic diversity of RVA strains [16,91].
Surveillance of G/P types in circulating RVA strains is important for vaccine development [23] and vaccination programs [21,85].Of note, G/P genotyping addresses the serotype specificity of these strains with respect to vaccines used on farms [21,85,92,95].Vaccines available in Argentina include a G6(IV)P [5] strain (UK), but few of them also contain a G10P [11] strain (B223).In our study, 86% of beef and 55.5% of dairy farms performed systematic vaccination of pregnant dams.Despite this practice, G6(IV)P [5] was the most prevalent combination among RVA strains circulating in beef farms.This finding is consistent with available data in France, where G6(IV)P [5] was prevalent either in vaccinated (76.2%) and non-vaccinated (85.7%) herds [23].Thus, vaccine-induced selective pressure is unlikely to occur in the beef herds of our region.However, the currently available vaccines may not offer protection against all circulating G/P types.Firstly, although we observed low circulation of G10P [11] strains, only cattle vaccinated with the B223 strain (G10P [11]) can develop neutralizing antibody responses against this strain if previous natural exposure has not occurred [96].In cattle and guinea pigs, antibodies induced by the vaccine strain UK (G6(IV)P [5]) are not protective against strain B223 (G10P [11]) [96].Secondly, the increase in emergent G/P combinations such as G8P [11] [this study, 95], G15P [11] [76,97], G5P [7] [42], and G24P [33] [14], or changes observed within the G6 genotype (G6P [11] strains) [this study, 21,44,92] could have an impact on vaccination programs, as the current vaccines may fail to protect [21,92].
During the study period, another interesting finding was the temporal variation of G/P genotypes in the dairy farm DR (Table 1).A single genotype combination (G6(III)P [11], G6(IV)P [5], and G10P [11]) was detected per year of sampling, suggesting yearly reintroductions.Circulating bovine RVA strains can show temporal fluctuations in their G/P genotypes [21,24].In Ireland, G6 and P [5] displayed a decrease from 2004 to 2005 and then increased in 2009, while G10 and P [11] were coupled with steady increases from 2004 to 2009 [21].According to our study, it is more likely to observe genotype fluctuations in dairy farms, and commercial bivalent vaccines (G6P [5] and G10P [11] strains) are expected to confer more protection than monovalent [21].
Up to date, a single bovine G8P [11] strain (Y136) has been sequenced for all of the 11 gene segments [39].However, only the VP7 gene was phylogenetically analyzed.The evolutionary analysis of strain 4385VT_D_BA showed that the VP3, VP4, VP6, and NSP1-NSP5 genes were very, or rather, closely related to the respective genes of RVA strains detected in calves from Argentina [76,77], Brazil [39], and Uruguay [14].However, the VP7 and VP1 genes showed a very close relationship with the respective genes of the caprine bovine-like strain 0040 (G8P [1]), which was detected in a goat kid from a dairy farm in Argentina [54].Also, the VP2 gene was rather closely related to that of strains 562 (G8P [14]) and 1115 (G8P [1]), which were obtained from Peruvian alpacas [78], as suggested by high similarity at the nt and aa level.Taken together, genomic data suggest at least two hypotheses explaining the origin of strain 4385VT_D_BA.Firstly, it could be a typical bovine RVA strain because its genotype constellation resembles those found in bovine strains [33,38,39,76], eight genes are phylogenetically derived from bovine strains, and two genes (VP7 and VP1) are closely related to those of the caprine bovine-like strain 0040, which is suspected to be a direct interspecies transmission from cattle [54].However, the VP2 gene was rather closely related to strains from alpacas, which are South American camelids that do not inhabit Argentina.Secondly, strain 4385VT_D _BA could represent a typical bovine strain with at least one gene (VP7) originating through a reassortment event, between bovine and artiodactyl bovine-like strains after interspecies transmission (e.g., cattle, goats, and guanacos) [54,76,78].This finding is not surprising because the evolutionary analysis of the E12-NSP4 genotype revealed that the guanaco strain Rio_Negro (G8P [1]) was ancestral with respect to all RVA strains from Argentina, Uruguay, and Paraguay.However, regarding VP1 and VP2 genes, more sequences from South American bovine RVA strains should be generated to obtain conclusive data.
At present, the intragenotype diversity of G8 is represented by six lineages [33,36,38,40], but a new lineage, tentatively designated G8(VII), has been proposed previously and comprises human (492SR, G8P [1]) and bovine (Y136, G8P [11]) strains from Paraguay and Brazil, respectively [39].Additionally, phylogenetic studies of human G8 strains request extensive sequence analysis of the VP7 gene in order to establish the lineages correctly [40].In this paper, we conducted a comprehensive genetic analysis of the G8 genotype, and eight lineages [G8(I) to G8(VIII)] were observed in the phylogram.In the pairwise distance frequency graph, the genetic distance values within and between lineages showed minor overlapping.However, it should be noted that previous studies and the current one have assigned certain human P [4] and P [14] strains (AU109, BP1062, and AS970) within the G8(V) lineage [33,[38][39][40].The high genetic distance values (10-12%) observed between these and other G8(V) lineage strains explain this limited overlap.Future studies should evaluate these strains as a possible new G8 lineage.Interestingly, after our analysis, the G8(VII) lineage included an alpaca strain (1115, G8P [1]) from Peru, which was reported recently [78].Moreover, the G8(VIII) lineage proposed in the current study mainly grouped human G8P [14] strains from Greece.As expected, the collection of Argentinean G8 strains from artiodactyl species [53,54,76,78], which also included the bovine study strain, grouped together in the same branch of the G8(II) lineage, supporting a common origin for VP7.
Epidemiologically, few studies have evaluated in detail the occurrence of particular lineages of G6, G10, P [5], and P [11] genotypes.In Brazil [25], France [23], Ireland [21], the Netherlands [43], and Spain [13], G6(IV) (40-100% of strains) and G6(III) (12-20% of strains) are the most prevalent intragenotype lineages among bovine strains.In this study, the ML phylogenetic analysis of the VP7-G6 sequences revealed a particular pattern.The G6(IV) lineage comprised all P [5] strains from beef calves and only two from dairy ones, whereas the G6(III) lineage included all P [11] strains from dairy calves.This segregation according to the P genotype and production system (beef or dairy) agrees with previous studies conducted in Argentina [41] and Brazil [25].It is unknown which factors cause this observation, but others also report a genetic linkage between G6(IV) and P [5], or G6(III) and P [11] [14,23,25,42].The VP4 protein of the spike has many functions including binding and fusion to cell membrane, infectivity, virulence, and neutralization (antigen) [7].Therefore, these features and host factors (i.e., breed of calves) may explain the restriction of P types according to the G6 lineage and production system.
Interestingly, the G6(IV) strains from the subregion "Cuenca del Salado" and Uruguay [14] were rather closely related, and this finding has not been observed before in South America [14,25,41,92,104].Probably, this genetic relationship could be due to the geographic proximities between this region and Uruguay.Interestingly, the G6(IV) strain 5595DR_D_BA, which was obtained from a diarrheic dairy calf, was phylogenetically closely related to the North American bovine strain B641.Up to date, none of the South American strains displayed this phylogenetic feature, at least for VP7 [25,41,44,92,104].Regarding the G10 genotype, 10 lineages (I-X) are recognized [105].The G10 strains reported herein clustered within the G10(VI) lineage together with Argentinean [41] and Uruguayan strains [14].In contrast, the Brazilian bovine G10 strains clustered apart within G10(III) and G10(IV) lineages [25].Concerning the P [5] and P [11] genotypes, most of the strains characterized in this study clustered in the P [5](I) and P [11](VI) lineages.Of note, due to the increasing number of P [11] sequences, their lineage classification might have to be updated as three strains from this study could not be included in the previously proposed six lineages (I-VI) [25,41].In this regard, the Uruguayan P [11] strains [14] together with strain 3DR_D_BA (G10P [11]) may represent a new lineage.
Although phylogenetic studies of circulating bovine RVA strains have been performed [13,14,20,22,23,25,42,43,90], none of them assessed the linkage of RVA strains according to the production system (beef or dairy) and the location of farms in production regions.The VP7 and VP8* sequences of the highly prevalent G6(IV)P [5] and G6(III)P [11] strains were analyzed; the sequences obtained in this and previous studies were included [41,44].Both genes were useful to assess the phylogenetic relationship of G6(IV)P [5] strains according to the production system and farm location.In the future, this epidemiological approach should be considered in the phylogenetic studies of RVA strains infecting domestic animals.
The spatio-temporal Bayesian inference of American bovine G6 RVA strains was performed using three VP7 gene sequence datasets.The RVA-G6 dataset, which included G6(III) and G6(IV) lineage strains together, showed that the diversification of bovine G6 RVA in the American countries would have started in the mid-1850s in the USA, almost a century before its first description as a cause of NCD in that country [45].As suggested by the phylodynamic and phylogeographic data, G6 RVA would have spread from North to South America and expanded its genetic diversity in South American countries.However, the currently available VP7 sequences from North American strains are scarce and may bias our interpretations; therefore, further sequences are needed.
The ancestral location of the G6(IV) lineage would be the USA around 1940, but two large branches in the MCCT started to evolve in Brazil and then in Argentina.This finding should be interpreted with caution because, as discussed previously, the VP7 sequences of dated strains from the USA represented only 1.6%, whereas Brazilian and Argentinean sequences comprised 47.2% and 34.6%, respectively.Introductions of RVA in Brazil from other countries in addition to the USA should not be discarded; indeed, our analyses suggest virus introductions in Brazil from Argentina and Venezuela.The age and origin of the nearest common ancestor of the G6(III) lineage were around 1980 in Argentina.This finding is expected because this lineage was first reported in Argentina [44] and thereafter in Brazil [110] and Uruguay [14].As indicated by phylogeographic data, the route of spread from Argentina to its neighboring countries seems plausible.In agreement with this finding, Argentina is suggested as the most probable origin and route of introduction for BCoV in Uruguay [111].

Conclusions
RVA infections are still a frequent cause of diarrhea in neonatal calves and can have a greater health impact on beef farms.The G/P genotypes, genetic diversity, and phylogeny of circulating bovine RVA strains can differ according to the production system and the location of farms.Thus, future studies should evaluate which factors may bias the fitness of RVA field strains in beef and dairy calves and farms.Also, our research provides further novel knowledge about the evolutionary origin of the G8 strains infecting different animal species.Finally, in American countries, the phylodynamic/phylogeographic data generated in this study support how bovine G6(IV) strains emerged first in the USA, nearly a century before its first report in calves from the same country.Conversely, bovine G6(III) strains emerged later in Argentina and were introduced in Uruguay from this neighboring country.

Figure 1 .
Figure 1.Distribution of fecal scores among sampled calves according to the clinical condition (diarrheic and non-diarrheic) and production system (beef and dairy farms).

Figure 1 .
Figure 1.Distribution of fecal scores among sampled calves according to the production system (beef and dairy farms) and clinical condition (diarrheic and non-diarrheic).

Figure 3 .
Figure 3. Phylogenetic analyses of G8 RVA strains.(A) Phylogenetic tree of the G8 genotype.The phylogram was reconstructed using the Kimura 2-parameter substitution model and the neighborjoining method.For better visualization, lineages are surrounded by black lines and strain names are abbreviated.The G8 linages are designated according to previous studies [38,40] using Roman numbers.Argentinean G8 strains from this and previous studies are highlighted with an irregular blue pentagon.(B) Comparison of genetic distances within (intra-lineage) and between (inter-lineage) the G8 lineages.Genetic distance values, expressed as percentages, were calculated using the pairwise distance algorithm and the p-distance model (MEGA 7.0).(C) Frequency graph of pairwise genetic distances.The percentage values of genetic distances are indicated in the abscissa (x-axis) and the frequencies of these values are indicated in the ordinate (y-axis).The cutoff value to best discriminate lineages is shown by a dotted line.

Figure 4 .
Figure 4. Phylogenetic trees of bovine RVA strains from Buenos Aires province, Argentina.(A) VP7 sequences of RVA strains corresponding to G6(III) and G6(IV) lineages.(B) VP8* sequences of RVA strains with P[5] and P[11] genotypes.Phylograms were reconstructed using the Kimura 2-parameter substitution model and the neighbor-joining method (MEGA 7.0 software).The phylogram of VP7-G6 lineages was rooted using strain RVA/Cow-wt/ARG/B383/1998/G15P[11].The analysis only included similar sequences and with available information.Each strain is represented with a number and symbol; white symbols correspond to beef farms, whereas black symbols correspond to dairy farms, as follows: "Cuenca del Salado"; "Depresion de Laprida"; "Mar y Sierras"; "Oeste";• "Abasto Norte"; "Abasto Sur".Bootstrap values (1000 replicates) >70% are shown as branch nodes.The name of the strain, type of production system, and province of collection are indicated.Abbreviations: D/B, Dairy or Beef; BA, Buenos Aires.(C) The map of Buenos Aires province shows RVA strains according to the production system (dairy or beef) and the location of farms in districts and production regions; each strain is represented with numbers and symbols as indicated previously.The location of the farms does not represent GPS coordinates.

Figure 5 .
Figure 5. Phylodynamic and phylogeographic reconstruction of bovine G6 RVA strains in American countries.(A) Maximum clade credibility tree (MCCT) of G6 lineage III [G6(III)] and G6 lineage IV [G6(IV)] strains together.The color of the branches represents the most probable country of origin.Posterior probability values of main clades are shown.Time of the most recent common ancestor (TMRCA) and the 95% highest posterior density (95% HPD) intervals are indicated below the tree nodes.Country abbreviations using the three-letter code (alpha-3) (ISO 3166): ARG, Argentina; BRA, Brazil; CAN, Canada; MEX, Mexico; URY, Uruguay; USA, the United States of America; VEN, Venezuela.(B) Inferred spread pattern of G6(III) and G6(IV) lineage strains in American countries.Line color represents the relative strength of connection between countries according to the Bayes Factor (BF) test: red arrows, decisive support with BF > 1000; blue arrows, strong support with 100 < BF < 1000; and green arrows, supported rates with 3 < BF <100.

Figure 5 .
Figure 5. Phylodynamic and phylogeographic reconstruction of bovine G6 RVA strains in American countries.(A) Maximum clade credibility tree (MCCT) of G6 lineage III [G6(III)] and G6 lineage IV [G6(IV)] strains together.The color of the branches represents the most probable country of origin.Posterior probability values of main clades are shown.Time of the most recent common ancestor (TMRCA) and the 95% highest posterior density (95% HPD) intervals are indicated below the tree nodes.Country abbreviations using the three-letter code (alpha-3) (ISO 3166): ARG, Argentina; BRA, Brazil; CAN, Canada; MEX, Mexico; URY, Uruguay; USA, the United States of America; VEN, Venezuela.(B) Inferred spread pattern of G6(III) and G6(IV) lineage strains in American countries.Line color represents the relative strength of connection between countries according to the Bayes Factor (BF) test: red arrows, decisive support with BF > 1000; blue arrows, strong support with 100 < BF < 1000; and green arrows, supported rates with 3 < BF <100.

Table 1 .
Genotypes and G6 lineages detected in RVA strains infecting beef and dairy calves.