Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Spatio-temporal distribution and international context of bovine viral diarrhoea virus genetic diversity in France

Veterinary Research volume 55, Article number: 129 (2024) Cite this article

Abstract

Bovine viral diarrhoea (BVD) is one of the most economically damaging livestock enzootic diseases in the world. BVD aetiological agents are three pestiviruses (BVDV-1, -2 and HoBi-like pestivirus), which exhibit high genetic diversity and complex transmission cycles. This considerably hampers the management of the disease, which is why eradication plans have been implemented in several countries. In France, a national plan has been in place since 2019. Our understanding of its impact on the distribution of BVDV genotypes is limited by the availability of French genetic data. Here, we conducted a molecular epidemiology study to refine our knowledge of BVDV genetic diversity in France, characterise its international relationships, and analyse national spatio-temporal genotypic distribution. We collated 1037 BVDV-positive samples throughout France between 2011 and 2023, with a greater sampling effort in two major cattle production areas. We developed a high-throughput sequencing protocol which we used to complete the 5’UTR genotyping of this collection. We show that two main BVDV-1 genotypes, 1e and 1b, account for 88% of genotyped sequences. We also identified seven other BVDV-1 genotypes occurring at low frequencies and three BVDV-2 samples (genotype 2c). Phylogenetic analyses indicate different worldwide distribution patterns between the two main BVDV-1 genotypes. Their relative frequencies present no major changes in France since the 1990s and few variations at the national scale. We also found some degree of local spatial structuring in western France. Overall, our results demonstrate the potential of large-scale sequence-based surveillance to monitor changes in the epidemiological situation of enzootic diseases.

Introduction

Bovine viral diarrhoea (BVD) is a major infectious disease of cattle that affects virtually all livestock producing countries worldwide [1,2,3]. The disease can lead to an important reduction in zootechnical performance, resulting in significant but heterogeneous economic impact among countries, estimated at up to 687 US dollar per animal in an infected herd [4]. Its huge socio-economic burden has led seventeen European countries to implement a control plan over the past thirty years [1, 4].

BVD can be transmitted horizontally from infected animals or fomites, resulting in a transient infection with various clinical outcomes in terms of symptoms and severity. They may include digestive, respiratory, reproductive symptoms and immunosuppression, but are typically mild [5]. BVD vertical transmission can occur and result in the birth of a persistently infected (PI) calf, depending on the stage of gestation at the time of infection [6]. PI animals propagate and sustain the disease and are essential to BVD endemicity [5].

BVD etiological agents are viruses from three species of the Pestivirus genus in the Flaviviridae family: BVD virus 1 (BVDV-1), BVDV-2 and HoBi-like pestivirus, currently recognized by the nomenclature of the International Committee on Taxonomy of Viruses [7] as Pestivirus bovis, Pestivirus tauri and Pestivirus brazilense respectively. Pestiviruses are single-stranded positive-sense enveloped RNA viruses of approximately 12.3 kb. Their genomes consist of one open reading frame coding for a polyprotein, flanked by two untranslated regions (UTR) [8]. Extensive molecular epidemiology studies mainly based on 5’UTR, Npro and E2 partial regions have highlighted the subdivision of BVDV-1 into more than 24 genetic lineages (genotypes 1a to 1x) and the subdivision of BVDV-2 and HoBi-like pestivirus into 5 lineages each [2, 9,10,11].

In France, only three molecular epidemiology studies have been published, showing the predominance of BVDV-1 with two main genotypes (1e and 1b) [10, 12, 13], the sporadic presence of BVDV-2 [2, 10, 14] and the absence of HoBi-like pestivirus [10]. Our knowledge of the long-term trends of BVDV-1 genetic diversity in France and its spatial distribution is limited, although it could unveil different dynamics in the circulation of BVDV genotypes and help to understand BVD distribution. This information would be important as an eradication plan has been in place in France since August 2019 [15]. It relies on the identification and management of PI animals, either through the serological monitoring of herds by screening tank milk or blood samples, or by testing directly for viruses by qRT-PCR in ear notches from newborn calves (sample taken within twenty days of birth). In this context, monitoring BVDV genetic diversity is crucial to describe the genotypic distribution before and during the plan. It is also important to ensure that the surveillance and control methods in use, such as diagnostic tests and vaccines, remain appropriate considering the genotypes in circulation.

Here, we report BVDV 5'UTR genotyping data obtained from infected cattle in most regions of France from 2011 to 2023. We describe BVDV genetic diversity in France and its international relationships, and analyse the spatio-temporal genotypic distribution.

Materials and methods

Sample collection

French BVDV-positive samples were obtained from four different sources (Figure 1 and Additional file 1). All samples have been diagnosed by regional veterinary analysis laboratories. Positivity was assessed mainly with PCR-based tests and, in few cases, with antigen ELISA-based tests. The first dataset (called dataset 1) represents BVDV-positive samples collected by the BVD national reference laboratory (NRL) between 2011 and 2022. Samples originated from most of the cattle-producing departments (2nd French administrative level, which corresponds to Nomenclature of Territorial Units for Statistics French—NUTS 3 level). It comprises samples collected as part of routine surveillance between 2011 and 2020 and targeted surveillance between 2021 and 2022 with ten random samples from different herds per department. Routine surveillance samples are of various origins in terms of tissues and animals. Targeted surveillance samples originate from ear notches collected as part of the eradication plan in France. The second dataset (called dataset 2) comprises BVDV-positive samples collected between 2016 and 2017 in the Normandie region (French NUTS 1 level). Samples were mostly blood and were isolated from both calves and adults, PI or not, in a variety of circumstances, including cattle introductions, sanitation plans, certifications or specific diagnostics commissioned by veterinarians. The two last datasets (called datasets 3 and 4) comprise BVDV-positive ear notches isolated from newborn calves as part of the French national BVD eradication plan. They originated from Auvergne-Rhône-Alpes (AuRA, dataset 3) and Bretagne (dataset 4) regions, collected during the periods August 2019—April 2023 and September 2020—April 2023, respectively. Additionally, BVDV-positive serum samples were collected in AuRA region in 2021. Each sample of our collection is associated with a location and date of sampling (Figure 1 and Additional file 1). For dataset 1, only the department of origin and the sampling year are available. For datasets 2, 3 and 4, we were able to retrieve the farm postcode and the sampling date for most of the samples. For dataset 4, the sampling date was approximated by the birth date of the calves.

Figure 1
figure 1

Geographical and temporal distribution of Bovine viral diarrhoea virus samples in this study and cattle density. A Map of France showing department administrative subdivision. The colour gradient indicates cattle density km−2 extracted from Robinson et al. [56]. B Spatial distribution of sample collection. Black circles represent the number of samples collected in each department. The map is coloured according to French administrative regions and was plotted using the R package sf [79]. C Temporal distribution of samples. The number of samples over time is presented in two panels, according to the level of precision of the metadata available for each sample: month and year of collection (top panel) or period only (bottom panel). Bar charts are coloured according to French administrative regions.

Full size image

BVDV 5’UTR amplification and sequencing

All sequences obtained in this study are based on PCR amplification of 5’UTR using primers 326R and 324F [16]. RNA extraction step is detailed for each dataset in Additional file 2. Four distinct 5’UTR sequencing protocols were used. All 5’UTR amplicons obtained from samples in datasets 1 and 2 were sequenced using Sanger technology, as were 96 samples in dataset 3 (Additional file 2). 5’UTR sequences of the remaining samples in dataset 3 and all samples in dataset 4 were generated using a high-throughput sequencing protocol based on Illumina technology as described below.

Nucleic acids extracted from ear notches from datasets 3 and 4 (see Additional file 2 for details) were treated to remove DNAs using DNaseI (Invitrogen) at room temperature for 15 min, followed by inactivation with EDTA at 65 °C for 10 min. The resulting purified RNAs were reverse transcribed with SuperScript III reverse transcriptase (Invitrogen) using primer 326R [16] as described by the manufacturer (Invitrogen, 18,080–051) on a T100 PCR thermal cycler (BioRad). A high-throughput Illumina sequencing library preparation protocol was adapted from the 16S metagenomic sequencing protocol of the manufacturer (Illumina, 770–2018-009-B). In brief, it consists of using primers composed of a locus-specific sequence, i.e. primers 326R and 324F [16], and Illumina adapters added to primer sequence extremities (Table 1) for direct sequencing of short amplicons on an Illumina machine without DNA tagmentation. PCR amplification was performed using the modified primers and Q5 polymerase enzyme (New England BioLabs). The reaction was carried out using the temperature cycle described in Additional file 3. Samples were purified using Macherey–Nagel NucleoMag NGS Clean-up and Size select with a beads ratio of 3:2 and normalised to 0.2 ng/µL prior to the next step. Illumina dual-index barcodes were added for multiplexing amplicons with the Nextera™ XT DNA Index Kit using a limited cycle PCR (see Additional file 3). Quality and average size of the libraries were checked using an Agilent 2100 Bioanalyzer High Sensitivity DNA chips following the manufacturer’s recommendations. Each library was purified as described above, quantified in Qubit with a Qubit dsDNA HS assay kit (ThermoFisher), normalised and then mixed to obtain a pool of 96 libraries at 1 nM. BVDV libraries were diluted to 70 pM and pooled with 15–20% of PhiX (Illumina) 100 pM or a 70 pM tick mitochondrial DNA library (data not published) in order to increase nucleotide diversity. These additional libraries had a higher genetic diversity than 5’UTR amplicons and help to provide a more balanced signal for the instrument. The latter library was obtained using the same protocol as presented here. Three sequencing runs of 96 samples each were performed on an Illumina iSeq100 benchtop sequencer using 20 µL of the pooled libraries (paired-end runs and 151 cycles).

Table 1 5’UTR primers adapted from Vilček et al. [16] (in bold) with Illumina adapter sequences.
Full size table

Assembly from high-throughput sequencing

Raw read quality was assessed with FastQC v0.12.1 [17] and quality reports were aggregated using multiQC v1.12 [18]. Reads shorter than 10 nucleotides were discarded. Primers, adapter sequences and flanking N bases were trimmed using cutadapt v4.2 [19]. Quality trimming was carried out using Trim-Galore v0.6.7 [20] with a quality threshold of 30. All samples were aligned to the 5’UTR of NADL reference strain (GenBank accession number M31182). The reference sequence was indexed using Burrows–Wheeler Aligner (BWA) v0.7.17 and samtools v1.16.1 [21, 22]. Reads were mapped to the reference sequence with the BWA-MEM algorithm. A Bayesian approach to variant calling was performed using FreeBayes with ploidy set at 1 and a min-coverage of 10 [23]. VCF files were filtered using VcfFilter v1.0.9 to remove any site with estimated QUAL inferior to 20. Coverage was assessed from bam files using bedtools v2.30.0 genomecov -bga option [24]. Consensus sequences were produced with bcftools v1.16 and positions with depth coverage less than 10X were masked.

Analysis was conducted using R v4.3.0 [25] and Rstudio [26]. This pipeline is implemented in Snakemake v7.24.0 [27]. The following Snakemake wrappers were included: fastqc (v1.3.2), multiqc (v1.3.2), cutadapt/pe (v1.3.2), bwa/index (v1.19.1) bwa/mem (v1.3.2) samtools/index (v1.3.2) and freebayes (v1.3.2).

Public BVDV-1 sequence collection

We collated all Pestivirus bovis (BVDV-1) genetic data publicly available from the National Center for Biotechnology Information (NCBI) on 24th August 2023. Using in-house scripts based on Python module BioPython [28], sequences were searched using the taxon identifier of Pestivirus bovis (id: 2170080) and all his descendants except BVDV-2a (id taxa: 1855269) and BVDV-2b (id taxa: 2590438) as determined by the NCBITaxa function of the ete3 module [29]. A total of 21,989 BVDV-1 sequences were retrieved, covering mainly the 5’UTR (16 933 with “5’UTR” in GenBank description). Sequences less than 100 nucleotides long were removed (n = 83), resulting in a dataset of 21 906 public sequences.

We selected a subset of 113 sequences, of which 104 represented genotypes 1a to 1x  for BVDV-1 and 9 represented genotypes 2a to 2e for BVDV-2 (Additional file 4). Some recently assigned genotypes have a conflicting letter code and were distinguished by the country of publication. These are genotypes 1l_Turkey [30] and 1l_France [12]; 1r_Italy [31] and 1r_Turkey [32]; and 1v_Turkey [33] and 1v_China [34].

Another subset of 214 publicly available sequences was selected comprising only sequences isolated in France from cattle. It included 6 sequences from Vilček et al. [13] corresponding to 5 isolates, 61 sequences from Jackova et al. [12] corresponding to 47 isolates, 145 sequences from Rivas et al. [10], as well as two others sequences: FJ493479 and JQ994205. Sequences from different genomic regions with the same isolate information were concatenated, resulting in a dataset of 199 sequences.

BVDV-1 public sequence genotype assignment

An exhaustive genotype assignment of publicly available BVDV-1 sequences was performed. Taxonomic assignment was carried out in three steps. Briefly, the 21 906 public sequences were clustered using the CD-HIT program [35] with a similarity threshold of 99%. A total of 5636 clusters were defined and the representative sequence of each cluster was retrieved. Batches of 100 representative sequences were aligned to the reference sequences of BVDV-1 genotypes (Additional file 4) using MAFFT version 7.490 [36] with add and localpair parameters. Maximum likelihood (ML) phylogenetic trees were inferred for each batch using IQ-TREE v.2.2.0 [37] with the nucleotide substitution model GTR + F + R7 and 1000 ultrafast boostrap (UFBoot) [38] (minimum correlation coefficient for UFBoot convergence: 0.95). The genotype of each representative sequence was assigned based on its inclusion in genotype clades of reference sequences using ML trees and the mrca function in the R ape package [39]. Sequences not included in genotype clades of reference sequences were assigned to the genotype of the closest reference sequence according to cophenetic distance matrix (cophenetic.phylo function in ape R package [39]) and branch statistical support (UFBoot). Genotypes were assigned to 5322 representative sequences, and by extension, to corresponding clusters, resulting in a total of 21 261 genotyped sequences. To reduce redundancy per genotype, sequences were finally clustered by country (or “unknown” in case of missing values) with a similarity threshold of 99.5% using CD-HIT program.

Phylogenetic analyses

Four multiple sequence alignments (MSA) were generated. One alignment included reference sequences, previously published French sequences and sequences of this study. The other three MSA were created specifically for genotypes 1e, 1b and 1d, and included publicly available sequences and sequences of this study assigned to each genotype. Sequences were aligned using MAFFT version 7.490 [36]. MSA were manually curated using Aliview v1.28 [40]. IQ-TREE v.2.2.0 [37, 41] was used to infer maximum likelihood phylogenies based on each of the four alignment with the best-fitting nucleotide substitution model and parameters estimated by ModelFinder embedded in IQ-TREE. Generalised Time Reversible (GTR) substitution model was selected for all datasets with an empirical state frequency (F) and a free rate heterogeneity (R) with different categories depending on the alignment (GTR + F + Rx). Ultrafast bootstrap procedure [38] and Shimodaira–Hasegawa approximate likelihood ratio test (SH-aLRT) [42] were used to assess statistical branch supports for the first dataset. Phylogenetic trees were visualised using R package ggtree [43]. A root-to-tip analysis using TempEst [44] was performed to measure the molecular clock signal contained in the 5’UTR sequence data. The linear regression between sample collection dates and root-to-tip genetic distances showed very weak correlation (R2 < 0.0005) indicating that a molecular clock approach based on 5’UTR sequences would not be robust.

Spatial and temporal statistical associations

Fisher exact tests (FET) [45] were used to test the independence between sampling periods or geographical areas and genotype relative frequencies. For large contingency tables (i.e., periods vs. genotype relative frequencies table), p-value was computed by Monte-Carlo simulation using fisher.test function of stats R package. When a significant association was identified (p-value ≤ 0.05), a post-hoc test was performed to determine the significance of the pairwise differences using fisher.multcomp function from RVAideMemoire R package [46]. P-values were adjusted using Benjamini–Hochberg correction. Spatial structure for samples with known farm postcode was assessed using Mantel correlograms, which estimate the correlation between genetic and geographic distances. Correlations were based on patristic distances between French samples with known farm postcode in the phylogenetic tree (Additional file 6) and geographic distances between centroids of corresponding postcodes. Mantel correlograms were computed using the mantel.correlog function from the vegan package, with Benjamini–Hochberg correction (options progressive = T and cutoff = F) [47]. Distance classes were determined using Sturge’s rule.

Results

Sample collection

The study included 1037 BVDV-positive samples collected in France between 2011 and April 2023 (Additional file 1). They originated from different sources (referred to as datasets 1–4, see Materials and Methods) associated with different sampling schemes. These schemes were either opportunistic (datasets 1 and 2) or random in specific geographical areas (datasets 1, 3 and 4). Most French cattle-producing regions are represented (Figures 1A, B) and we were able to collect more samples from three French regions (AuRA, Bretagne and Normandie), which are part of two main cattle production areas in France (Figures 1A, B). Due to the diversity of our sampling sources, the spatial distribution of samples is not homogeneous over time (Figure 1C). This is particularly the case for samples isolated in Normandie which were mainly collected in 2016 and 2017. The number of samples also varies greatly from one year to another, with most samples collected from 2020 onward as part of the eradication plan (Figure 1C).

5’UTR high-throughput sequencing

A high-throughput 5’UTR sequencing protocol was applied on 282 samples isolated from ear notches, included in datasets 3 and 4 (Additional file 1). These samples have diverse viral loads, with qPCR cycle threshold values (Ct) ranging from 18.28 to 34.89. The number of reads mapped per sample ranged from 3 to 438,597 (median of 22 548.5). The mean depth coverage varied from 1.3 to 213 795 (median of 9368.7). Most of the samples (270 out of 282, 96%) covered > 80% of the reference sequence. The remaining 12 samples had low coverage (< 65%), including 7 samples with no coverage, and were therefore not included in further analyses. These samples were sequenced in the same sequencing run (run 3) which showed a slight decrease in the number of mapped reads per sample compared to the other sequencing runs (Additional file 5A) and a much lower percentage of mapped reads (Additional file 5B). The lower efficiency of this third run was likely due to a high concentration of primer dimer. This is consistent with the mean read lengths of the three runs, which are 141 nucleotides (nt) (run 1), 142 nt (run 2) and 79 nt (run 3). Overall, a correlation is observed between the percentage of mapped reads and initial viral loads (Additional file 5B).

BVDV genetic diversity in France

Of the 1,025 BVD-positive samples successfully sequenced in this study, only three samples were identified as BVDV-2 (genotype 2c), with all the others belonging to BVDV-1 (Additional file 1). This is consistent with previous studies highlighting the rarity and absence of BVDV-2 and Hobi-like pestivirus in France, respectively [10, 12,13,14]. Consequently, we focused exclusively on BVDV-1 in the subsequent analyses. The 5’UTR sequences of the 1,022 successfully sequenced BVDV-1 samples were aligned with 199 previously published French sequences and 104 reference sequences defining BVDV-1 genotypes (including 9 from France representing genotypes 1x and 1l_France). A maximum likelihood phylogenetic tree was inferred from this alignment of 1316 sequences (Additional file 6). The tree showed that the French sequences of this study were mostly divided into five distinct genotypes, including 580 sequences in 1e (57%), 327 in 1b (33%), 71 in 1d (7%), 13 in 1a (1.3%), 12 in 1l_France (1.2%). The 19 remaining sequences fell into 1k, 1i, 1r, 1s clades and into undetermined intermediate genotypes (see Additional files 1 and 6). Previously published French sequences were distributed throughout the tree, suggesting that we did not uncover any major unreported BVDV-1 genetic diversity in France. Only genotypes 1k and 1i have not previously been reported in France, but they account only for four sequences. Three samples (BVDV1Fr0209, BVDV1Fr0262 and BVDV1Fr0916, see Additional file 1, highlighted by an asterisk on Additional file 6) were identical or almost identical to vaccine strains KE-9 (genotype 1b; Bovela, Boehringer Ingelheim) or Oregon C24 (genotype 1a; MUCOSIFFA, CEVA Santé Animale). An epidemiological survey confirmed that two out of the three samples were related to vaccinated animals (NRL for BVD, personal communication), consistent with the detection of these vaccine strains elsewhere [48]. No information was available for the last sample. The phylogeny shows an overall lack of support and a low phylogenetic resolution with numerous polytomies, especially for samples with close sampling dates and location in 1e clade.

French genetic diversity of BVDV-1 in a global context

The international context of the French BVDV-1 genetic diversity was investigated for the main genotypes identified in France (i.e. 1e, 1b, 1d, 1a and 1l; Figure 2, Additional files 7 and 8). We genotyped all publicly available BVDV-1 sequences on the NCBI database (n successfully genotyped = 21,261) and explored the most common location of samplings. Genotypes 1a and 1b were present in almost all cattle-producing countries, while genotype 1d has been mainly found in Europe (Additional files 7A and B). Genotype 1e has historically been sampled mainly in Western Europe, and genotype 1l_France has been identified almost exclusively in France (Additional file 7A and B).

Figure 2
figure 2

Phylogenetic relationships between French and international sequences of BVDV-1 genotypes 1e and 1b. Maximum likelihood phylogenetic trees inferred from French and international sequences of genotypes (A) 1e and (B) 1b. Tips are coloured according to the sample continent or country of origin and sizes represent the number of closely related publicly available sequences (99.5% identity) from the same location. European countries with less than ten sequences were grouped under the location “Europe” in 1e phylogeny.

Full size image

The phylogenetic tree showing relationships between French and international sequences (Figure 2A) revealed that genotype 1e can be divided into four monophyletic groups called clades 1e.1 to 1e.4 for clarity. Most of the French genetic diversity of genotype 1e fell into clade 1e.1 with 488 sequences (75.3% of total 1e French sequences), while clades 1e.2, 1e.3 and 1e.4 comprised 47 (7.3%), 35 (5.4%) and 67 (10.3%) French sequences respectively. Clade 1e.1 was characterized by an overall low phylogenetic resolution compared to the three other 1e clades, with in particular an important polytomy at internal nodes. French sequences predominate (54.2% of total sequences in clade 1e.1) and are distributed between internal and external nodes. Clade 1e.2 included mostly sequences sampled in Switzerland between 2008 and 2012 as part of the Swiss national eradication plan [49]. French sequences from several regions were scattered throughout the clade, mainly at external nodes. Clade 1e.3, which corresponded to genotype 1x  in Rivas et al. [10], and clade 1e.4 contained mostly French sequences distributed at external nodes. Overall, genotype 1e phylogeny showed a marked difference between clade 1e.1 and the rest of the tree, suggesting a primary contribution of this clade to BVD disease in France.

The phylogeny of genotype 1b (Figure 2B) including French and international sequences can be divided into two monophyletic groups called clades 1b.1 and 1b.2 for clarity, with most of the French sequences falling into clade 1b.2 (241 sequences, 58.6% of total 1b French sequences). Clade 1b.1 contained sequences from most of the continents, while clade 1b.2 sequences were comparatively more sampled in Europe. French sequences were distributed throughout the tree into dozens of monophyletic groups at external nodes. International sequences fell into both internal and external nodes. The tree showed a clear mix of French and international sequences, suggesting a high rate of international migration, including several introductions into France that have contributed greatly to BVD disease.

The phylogenetic tree corresponding to genotype 1d (Additional file 8) showed three monophyletic groups, called 1d.1, 1d.2 and 1d.3 for clarity. Clade 1d.2 includes most of the French sequences (70 sequences, 90.9% of total 1d French sequences), which clustered in one subclade that was not particularly closely related to international sequences. The rest of the French sequences are scattered throughout the tree. Given the small number of 1a and 1l_France sequences identified, no phylogenetic trees were reconstructed with international sequences for these genotypes.

Spatial and temporal distribution of BVDV-1 genotypes in France

We investigated whether the distribution of BVDV-1 genotypes was structured in space and time. We collated all genotyping data available for France, including the present work and previous studies [10, 12, 13] ranging from 1994 to 2023 (Figure 3A), and examined the potential change in genotype relative frequency over this time interval. Periods were defined on the basis of the least precise available temporal information. A Fisher exact test (FET) comparing the proportion of genotypes (that is, 1e, 1b and all other genotypes grouped as “others”) including all periods was significant (p-value = 0.036), but a post-hoc analysis revealed no significant pairwise comparisons between periods. This analysis suggests that there has been no significant shift in the proportion of BVDV-1 genotypes in France since the 1990s, with a higher degree of confidence for the period 2018–2023 due to a more extensive national sampling (Figure 1C and Rivas et al. [10]).

Figure 3
figure 3

Temporal and spatial distribution of BVDV-1 genotypes in France. A Proportion of genotypes over time in France. Genotyping data include current and previous studies (1994–1998: Vilček et al., 2001 [13], 2004–2005: Jackova et al., 2008 [12]; 2018–2020: Rivas et al., 2022 [10]). Rare genotypes (1r, 1s, 1k, 1i) were grouped. Numbers at the top of bar charts indicate the total number of samples for each period. B Proportion of each genotype per French department over the 2011–2023 period. Genotyping data includes current study and those from Rivas et al. [10]. Pie charts denote genotype proportions per department and sizes are proportional to the number of samples. Black boxes show the location of the three focused regions: Bretagne, Normandie and Auvergne-Rhône-Alpes.

Full size image

The proportion of BVDV-1 genotypes by department in France is shown on Figure 3B. Of the previous French genotyping data, only that of Rivas et al. [10] was included along with the present data, as it contains information on the department of sampling and is included in the time-frame of this study (2011–2023). Genotypes 1e and 1b were the main genotypes in most of the well-sampled regions (n ≥ 30, mean relative frequencies of 1e and 1b were 57% ± 6% and 32% ± 8%, respectively). We compared differences in the relative frequencies of genotypes between regions using a FET. Test was significant (p-value = 0.0005), and post-hoc tests indicated also 27 significant pairwise comparisons (Additional file 9A), more than half of which (17 out of 27) included Grand-Est or Bourgogne-Franche-Comté (BCF) regions, both located in the North-East part of France. The majority of significant comparisons involved “others” genotypes and 1e or 1b, suggesting that the ratio of the main genotype proportions was stable across well-sampled regions of France. The same test was applied within regions for 24 departments with at least ten samples localized in seven regions: Normandie, Bretagne, AuRA, Occitanie, Grand Est, Pays de la Loire and BCF. FET was significant only for BCF region (p-value = 0.0005). Post-hoc tests showed significant differences between three departments: Yonne (code: 89), Doubs (25) and Jura (39) (Additional file 9B). This suggests few variations in the relative frequency of genotypes on a national scale, with some regional specificities such as in the North-East of France. Nevertheless, we cannot exclude the existence of significant divergent trends that are not captured by our dataset.

Spatial structure of BVDV-1 genetic diversity

The geographical structuring of BVDV-1 genetic diversity was investigated at a regional scale for three well-sampled regions of the study (Bretagne, Normandie and AuRA, i.e. datasets 2, 3 and 4, Additional file 1). Only samples with known farm postcodes were included in the subsequent analyses. Maps of genotype proportions by farm postcode for these regions (Figures 4A, 4B and 4C) did not show strong geographical structuring of genotypes at the scale of a region or a department, although local clusters of similar genotypes can be observed, notably for Normandie. We therefore used Mantel correlograms to explore population structure relationship to space (Figure 4D). The correlations were very weak overall, and the correlograms shapes vary according to the region under consideration. Pairwise distances were significantly correlated for Normandie region (Figures 4B, D), positively in the first distance classes (up to 50 km) and negatively for longer distances (between 50 and 180 km). This describes an increasing genetic dissimilarity between pairs of samples as the geographical distance increases. Significant correlations were also observed for Bretagne region (Figures 4A, D) for smaller geographical scales: positively first (0–17 km), then negatively (34–51 km). AuRA region (Figures 4C, D) correlogram showed significant correlations for two distances classes, 12–25 km and 151–164 km, but cannot be interpreted. The lower magnitude of autocorrelation observed in AuRA region could be due to the longer time frame, that is four years of observations compared to two (Normandie) or three (Bretagne) years. To rule out this hypothesis, Bretagne and AuRA datasets were truncated to two years (2020–2022 for Bretagne and 2019–2021 for AuRA), but results showed similar trends compared to complete dataset results (Additional file 10). These results suggest that spatial structuring of BVDV-1 genetic diversity exist at local scales but vary depending on the region.

Figure 4
figure 4

Spatial distribution and local structure of BVDV-1 genotypes in three focused French regions. Proportions of BVDV-1 genotypes at local scale in (A) Bretagne, (B) Normandie (B) and (C) Auvergne-Rhône-Alpes. Pie charts denote genotype proportions per farm postal code and sizes are proportional to the number of samples. A small amount of random variation was added to each pie chart location to ensure anonymity. Grey lines outline department borders. D Mantel correlograms measuring correlations between phylogenetic and geographical distances for the three focused regions. Each dot is positioned in the middle of a distance class. Statistically significant correlations are indicated by a square.

Full size image

Discussion

Efficient high-throughput sequencing protocol for 5'UTR sequencing

In this study, we report a large number of BVDV 5'UTR genetic data (> 1000 samples) obtained from infected cattle in France between 2011 and 2023. This dataset includes sequences generated with a high-throughput sequencing protocol targeting Pestivirus 5’UTR that we have developed. This protocol is fast, reliable and can work on an Illumina lab-benchtop sequencing machine such as the iSeq100. With this high-throughput sequencing protocol, we successfully sequenced 270 out of 282 samples covering the complete targeted loci length. The small size of the 5'UTR and the efficiency of the Vilček et al. primers [16], make it possible to use this protocol to obtain robust sequencing when using low amounts of poor-quality samples, such as ear notches stored at −20 °C for several months. Given the mean coverage depth obtained here, the number of samples that can be multiplexed in one sequencing run is potentially high (at least 96 samples). These are two advantages compared to Sanger sequencing [50], which is the commonly used method to sequence 5’UTR BVDV amplicons [10, 48, 51]. Moreover, it requires only basic laboratory equipment to prepare the libraries, and is relatively inexpensive as no DNA tagmentation step during the sequencing library preparation is required. Finally, the conserved domains in the 5'UTR and the universality of the primers used enable monitoring of genetic diversity without a priori on virus species or genotype as it is possible to detect several pestiviruses [16]. Although the cost per sample is almost twice that of Sanger, the high sequencing success rate and the high level of multiplexing make our sequencing approach cost-competitive. This could facilitate future large-scale diversity studies of Pestivirus-related diseases.

Circulation of two main genotypes of BVDV-1 in France

Several Pestivirus species capable of infecting cattle are of interest for health management: P. bovis (BVDV-1), P. tauri (BVDV-2), P. brazilense (HoBi-like pestivirus), and P. ovis (Border Disease Virus) [52]. However, publication on Pestivirus genetic diversity in France is scarce [10, 12, 13]. In this study, we collected sequences of BVDV-positive samples with a fairly representative sampling of cattle distribution in France. In particular, we acquired samples from major cattle production areas: Normandie, Auvergne-Rhône-Alpes (AuRA) and Bretagne. To date, this type of extensive sampling for molecular typing has only been carried out and published in Switzerland [51] and the United-Kingdom [53]. Despite its low phylogenetic resolution [9, 54], 5’UTR has been used in the vast majority of previous taxonomic studies, which facilitate comparisons. Taxonomic assignment of the samples of this study based on partial 5’UTR sequences shows the presence of two major BVDV-1 genotypes (1e and 1b), three minor genotypes (1d, 1a and 1l_France, subsequently termed 1l) and several marginal genotypes (1s, 1r, 1i and 1k) in France. This genetic diversity is consistent with previous observations [10, 12, 13], with the exception of 1k (Bourgogne-Franche-Comté, 2011/2012) and 1i (Normandie, 2016/2017) which have not been reported to date in France. Genotypes 1r (Normandie, 2016/2017) and 1s (AuRA, 2019) were previously observed in the East of France (in 2018 and 2019 respectively [10]). Of note, clade 1x defined in Rivas et al. [10] based on 5’UTR sequences clustered here among 1e sequences. The weak branch supports observed in the phylogenetic tree, combined with the absence of clear criteria for clade delimitation, suggest that genome-based tree reconstructions would be needed to clarify the status of this clade. Only BVDV-1 and -2 (genotype 2c) were identified here. We did not find isolates belonging to the Pestivirus ovis species, although most diagnostic tests can detect both BVDV and BDV [55]. This could partially be explained by our sampling scheme, which intended to characterize the broad genetic structure of circulating BVDV. Indeed, the detection of P. ovis would have required a more specific sampling design with a focus on mixed cattle-sheep farming systems; and higher sheep densities are in general found in regions that were less sampled in our study, such as the South of France [52, 56].

Distinct international context between main French BVDV-1 genotypes

An exhaustive taxonomic assignment of publicly available 5’UTR sequences, combined with phylogenetic tree reconstructions, was performed to provide a better understanding of international virus flows. Genotype assignment analysis was complex due to the short length of studied sequences (dataset is composed of 75% 5’UTR sequences), errors in sequences and low phylogenetic signal of 5’UTR marker. These limits prevented us from constructing a well-supported phylogenetic tree directly from all publicly available sequences and using standard genetic clustering algorithm to confirm genotype assignment. Nevertheless, our results were overall consistent with the literature review of BVDV-1 genetic diversity by Yeşilbağ and collaborators [2]. Some inconsistencies could be due to differences with genotype assignment method in previous papers as reviewed in Yeşilbağ et al. [2], genetic diversity described but not publicly available such as 1l sample from Italy [57] or missing country annotation in GenBank (e.g., only nine publicly available sequences are annotated with Belgium as country of collection, whereas twenty are listed in paper by Yeşilbağ et al. [2]).

Each country seems to have its own genotype combination. France is characterized by a low number of genotypes and a predominance of the genotype 1e. This genotype has mainly been sampled in France and neighbouring countries, such as Switzerland [49] or Italy [57, 58]. Most of the French sequences group into one main clade (1e.1), that may have contributed primary to BVD infections in France. The potential contribution of under-sampled geographical areas to this clade cannot be ruled out due to the sampling bias in the publicly available dataset. Genotype 1b has been sampled in almost all cattle-producing countries. French sequences were distributed into different lineages within genotype 1b phylogeny, and grouped with sequences from all continents. This suggests several introduction events into France which gave rise to numerous transmission chains that contributed greatly to BVD disease. These results highlight different international context between genotypes 1e and 1b. The presence of rare genotypes (1s, 1k, 1r and 1i) in France is probably due to migrations from under-sampled French region or other countries. However, it is not clear whether their presence is due to the cryptic maintenance of transmission chains or to one-off events.

Migration rates between countries and their temporal dynamics are important to reconstruct transmission routes and understand the impact of cattle movements from abroad to BVD circulation dynamics in France. Due to the low phylogenetic resolution of 5’UTR marker, it was not possible to further explore the temporal dynamics of this clade in the present study. BVDV-1 could have been introduced into France on several occasions through different routes: mainly live cattle imports, closer contact between herds at borders, or import of contaminated semen. Changes in observed genetic diversity through the migration of one or more genotype could be expected as a result of large cattle movement from another geographical area with a different genetic diversity profile [59, 60]. However, the small number of live cattle imported into France [61] and the low probability of importing BVD-positive cattle put the current impact of this route into perspective. Indeed, only a handful of animals imported during the Switzerland eradication plan (2008–2012) was BVD-positive [62], even if this observation is to be considered with caution as “Trojan cows” (pregnant dams carrying a PI foetus) number can be underestimated [63, 64]. Close contacts at borders are hard to quantify without further data, but could be low as only Belgium and Luxembourg share a common border with France without a natural barrier such as the Rhine for Germany or mountains (Alps or Pyrenees) for Italy, Switzerland and Spain. Common grazing areas can be found in both mountain ranges but the extent of infection made possible by transhumance remains to be determined. German sequences of genotype 1d [64] did not particularly cluster with French sequences (Additional file 8). On the contrary, French and Swiss sequences are found relatively interspersed in clades 1e.1 and 1e.2 (Figure 2A), which supports the hypothesis of small-scale but frequent contacts between herds. A similar pattern could be expected with Italy due to the high proportion of 1e in North Italy [57, 58], but less sequences were retrieved or successfully genotyped. Finally, infections through semen are limited but possible [65], and imports could have impacted virus entry rates more likely decades ago as there were historical biosecurity recommendations but no mandatory controls on frozen semen [66]. Effective migration rate (i.e. migration events resulting in at least one local infection) is assumed to be much lower than real migration rate as PI presence in a herd tends to result in its seroconversion, providing herd immunity. This should limit the risk of introduction of another BVDV strain and facilitate the local maintenance of genotypes for years [67]. In addition, cross-reactivity have been demonstrated among BVDV-1 strains [68], and to a lesser extent between BVDV species [69]. Thus, effective migration rate into France should be low due to herd immunity for BVDV-1 genotypes and presumably higher for BVDV-2 and Hobi-like pestivirus. It is rather the lower international prevalence that may explain the rarity of these pestiviruses in France. However, phylogeny structure may reflect biases in available data and miss nodes related to key migration events. Publicly available sequences are biased geographically (some countries or areas are under-sampled), temporally (most data come from studies conducted over limited periods of time, as opposed to regular monitoring), in terms of local context (e.g. eradication plan or mitigation measures, such as vaccination) and sampling scheme. As a result, the clustering of French and international samples cannot be interpreted as direct evidence of migration as geographical, temporal and potentially genetic structures are confounded in the data, and the temporal aspect is not taken into account.

No major changes in temporal and spatial distribution of BVDV-1 genotypes in France

Statistical analysis suggests that there have been no clear changes in relative frequencies of genotypes 1e and 1b since the 1990s in France. This may indicate that the rate of introduction of BVDV-1 into France might not be high enough to replace established genotypes, or that such replacement through migration only occurred before the first molecular characterisation of circulating genotypes in France [13]. Furthermore, we showed that the relative frequencies of both main genotypes did not present major disparities, especially for the main cattle-producing regions. It could be the result of the great potential for virus dispersal through the dense cattle trade network. In addition, the presumably low effective migration rate could limit the introduction of other genotypes and slow the potential changes of genotype frequencies in the metapopulation of cattle herds.

France is currently implementing an eradication plan, which will presumably have an impact on the distribution of BVDV-1 genetic diversity. The plan should increase the removal of PI animals and therefore the rate of viral population extinction in subpopulations. Genetic diversity loss would then be expected, through an increase of genetic drift. The probability of losing a genotype would be directly linked to its frequency: the rarest genotypes would have the highest probability of being eliminated, while highly prevalent genotypes would take longer to be eradicated. This could explain why the main genotype tend to become even more predominant in the context of an eradication plan as seen in Switzerland [49] and, more strikingly, in Germany [48]. However, this effect is not yet visible in France, probably because of the gradual introduction of control measures and the latency in the implementation of the eradication plan, as well as the slow temporal dynamics of BVDV-1.

Region-specific local spatial clustering in France

The spatial structuring of BVDV-1 genetic diversity is overall weak. Low positive Mantel correlations between genetic and geographic distances are significant only at short-distance classes for Normandie and Bretagne region (Figure 4D). This is consistent with what was observed previously: no correlation was reported in Switzerland [70], only at farm level in Austria [71] and there was evidence for local spread in nearby farms in England and Wales [72]. These studies have suggested the importance of cattle movements in the absence of geographical structuring. More generally, trade movements seem to be a driver of BVDV lineages distributions [58]. However, trade movements data could not be obtained for this study. Differences among the three French regions could further be explained by several factors: regional husbandry practices, location, the sampling scheme, and applied control measures. Western France (Normandie and Bretagne regions) has historically had very different productions systems and cattle density (mostly dairy and high cattle density) to those in AuRA (relatively balanced between dairy and beef and moderate cattle density) [73, 74]. These two factors have been shown to be risk factors [75,76,77] and could explain differences in spatial structuring between western France and AuRA. However, our sampling scheme was not suitable for specifically investigating the effect of production type (milk or beef) on genetic diversity, as dairy areas were over-represented in the data. Furthermore, the central geographical location of AuRA region in France could also explain the low spatial structuring in this region due to the potential border effect, involving cattle exchange with surrounding departments. As coastal regions, this effect is less expected for Bretagne and Normandie regions. The higher positive correlations observed in Normandie could be partly explained by the sampling scheme, as samples from the same herd were more likely to be included in this dataset than in Bretagne and AuRA datasets. Indeed, it has been shown that viruses sampled in the same herd for the same outbreak are almost genetically identical [63], which could increase correlation signal for short-distance classes. Nevertheless, our results suggest that, in Normandie region during the studied period (2016—2017), BVDV-1 genotypes could persist locally through transmission among neighbouring herds with few long-distance dispersal events. The positive short-distance correlation observed in Bretagne could be linked to compliance with control measures since 1998 in Bretagne [78], which could limit the risk of virus introduction when cattle are purchased from another herd and reinforce the importance of contacts with neighbours. This is consistent with the work by Qi et al. [76] which highlights the importance of neighbourhood contacts in virus diffusion in a modelling study calibrated on Bretagne data.

Conclusion

This study presents an in-depth spatio-temporal analysis of BVDV-1 genetic diversity in France. It reveals distinct international patterns between the two main BVDV-1 circulating genotypes. Yet, we did not identify major changes in the proportions of these genotypes in France, either in space or over time. This study also shows that barcoding data generation, such as 5’UTR sequences, can be efficiently implemented to genetically monitor variants of an enzootic pathogen. However, the level of phylogenetic resolution associated with 5'UTR sequence data precludes advanced spatio-temporal analyses required to unravel BVDV-1 complex evolutionary and epidemiological dynamics. Genome-based studies would enable to explore the dynamics of the virus in greater details at different spatial and temporal scales, and make a better contribution to BVD disease management.

Availability of data and materials

All 5’UTR sequences from this study are available on NCBI under the accession numbers: PP513311 to PP514335. Farm postal codes for Bretagne and AuRA regions are available on request from the authors and after consultation with the GDS.

References

  1. Richter V, Kattwinkel E, Firth CL, Marschik T, Dangelmaier M, Trauffler M, Obritzhauser W, Baumgartner W, Käsbohrer A, Pinior B (2019) Mapping the global prevalence of bovine viral diarrhoea virus infection and its associated mitigation programmes. Vet Rec 184:711

    Article  PubMed  PubMed Central  Google Scholar 

  2. Yeşilbağ K, Alpay G, Becher P (2017) Variability and global distribution of subgenotypes of bovine viral diarrhea virus. Viruses 9:128

    Article  PubMed  PubMed Central  Google Scholar 

  3. Scharnböck B, Roch F-F, Richter V, Funke C, Firth CL, Obritzhauser W, Baumgartner W, Käsbohrer A, Pinior B (2018) A meta-analysis of bovine viral diarrhoea virus (BVDV) prevalences in the global cattle population. Sci Rep 8:14420

    Article  PubMed  PubMed Central  Google Scholar 

  4. Richter V, Lebl K, Baumgartner W, Obritzhauser W, Käsbohrer A, Pinior B (2017) A systematic worldwide review of the direct monetary losses in cattle due to bovine viral diarrhoea virus infection. Vet J 220:80–87

    Article  PubMed  Google Scholar 

  5. Evans CA, Pinior B, Larska M, Graham D, Schweizer M, Guidarini C, Decaro N, Ridpath J, Gates MC (2019) Global knowledge gaps in the prevention and control of bovine viral diarrhoea (BVD) virus. Transbound Emerg Dis 66:640–652

    Article  PubMed  Google Scholar 

  6. McClurkin AW, Littledike ET, Cutlip RC, Frank GH, Coria MF, Bolin SR (1984) Production of cattle immunotolerant to bovine viral diarrhea virus. Can J Comp Med 48:156–161

    CAS  PubMed  PubMed Central  Google Scholar 

  7. Walker PJ, Siddell SG, Lefkowitz EJ, Mushegian AR, Adriaenssens EM, Alfenas-Zerbini P, Dempsey DM, Dutilh BE, García ML, Curtis Hendrickson R, Junglen S, Krupovic M, Kuhn JH, Lambert AJ, Łobocka M, Oksanen HM, Orton RJ, Robertson DL, Rubino L, Sabanadzovic S, Simmonds P, Smith DB, Suzuki N, Van Doorslaer K, Vandamme A-M, Varsani A, Zerbini FM (2022) Recent changes to virus taxonomy ratified by the international committee on taxonomy of viruses (2022). Arch Virol 167:2429–2440

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Tautz N, Tews BA, Meyers G (2015) The molecular biology of pestiviruses. In: Kielian M, Maramorosch K, Mettenleiter TC (eds) Advances in Virus Research. Elsevier, Amsterdam, pp 47–160

    Google Scholar 

  9. de Oliveira PSB, Silva Júnior JVJ, Weiblen R, Flores EF (2021) Subtyping bovine viral diarrhea virus (BVDV): which viral gene to choose? Infect Genet Evol 92:104891

    Article  PubMed  Google Scholar 

  10. Rivas J, Hasanaj A, Deblon C, Gisbert P, Garigliany M-M (2022) Genetic diversity of Bovine Viral Diarrhea Virus in cattle in France between 2018 and 2020. Front Vet Sci 9:1028866

    Article  PubMed  PubMed Central  Google Scholar 

  11. Kalaiyarasu S, Mishra N, Jayalakshmi K, Selvaraj P, Sudhakar SB, Jhade SK, Sood R, Premalatha N, Singh VP (2022) Molecular characterization of recent HoBi-like pestivirus isolates from cattle showing mucosal disease-like signs in India reveals emergence of a novel genetic lineage. Transbound Emerg Dis 69:308–326

    Article  CAS  PubMed  Google Scholar 

  12. Jackova A, Novackova M, Pelletier C, Audeval C, Gueneau E, Haffar A, Petit E, Rehby L, Vilcek S (2008) The extended genetic diversity of BVDV-1: typing of BVDV isolates from France. Vet Res Commun 32:7–11

    Article  CAS  PubMed  Google Scholar 

  13. Vilček Š, Paton DJ, Durkovic B, Strojny L, Ibata G, Moussa A, Loitsch A, Rossmanith W, Vega S, Scicluna MT, Palfi V (2001) Bovine viral diarrhoea virus genotype 1 can be separated into at least eleven genetic groups. Arch Virol 146:99–115

    Article  PubMed  Google Scholar 

  14. Viard G (2022) Description atypique de BVD type 2 dans un troupeau de vaches laitières. JNGTV GTV2022 517–526. https://www2.sngtv.org/article-bulletin/description-atypique-de-bvd-type-2-dans-un-troupeau-de-vaches-laitieres/. Accessed March 2024.

  15. Ministère de l’Agriculture et de l’Alimentation (2019) Arrêté du 31 juillet 2019, fixant des mesures de surveillance et de lutte contre la maladie des muqueuses/diarrhée virale bovine (BVD). Journal officiel de la République Française 0177:45

    Google Scholar 

  16. Vilček Š, Herring AJ, Herring JA, Nettleton PF, Lowings JP, Paton DJ (1994) Pestiviruses isolated from pigs, cattle and sheep can be allocated into at least three genogroups using polymerase chain reaction and restriction endonuclease analysis. Arch Virol 136:309–323

    Article  PubMed  Google Scholar 

  17. Andrews (2010) FastQC: a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc. Accessed March 2022.

  18. Ewels P, Magnusson M, Lundin S, Käller M (2016) MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32:3047–3048

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Martin M (2011) Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J 17:10

    Article  Google Scholar 

  20. Krueger F, James F, Ewels P, Afyounian E, Schuster-Boeckler B FelixKrueger/TrimGalore: v0.6.7 - DOI via Zenodo. 10.5281/ZENODO.5127899. 2021.

  21. Li H (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv. https://doiorg.publicaciones.saludcastillayleon.es/10.48550/ARXIV.1303.3997

    Article  Google Scholar 

  22. Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, Whitwham A, Keane T, McCarthy SA, Davies RM, Li H (2021) Twelve years of SAMtools and BCFtools. Giga Sci 10:008

    Article  Google Scholar 

  23. Garrison E, Marth G (2012) Haplotype-based variant detection from short-read sequencing. arXiv preprint

  24. Quinlan AR, Hall IM (2010) BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26:841–842

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. R Core Team (2023) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. https://www.R-project.org/. Accessed Dec 2023.

  26. RStudio Team (2023) RStudio: Integrated Development Environment for R. http://www.rstudio.com/. Accessed Dec 2023.

  27. Mölder F, Jablonski KP, Letcher B, Hall MB, Tomkins-Tinch CH, Sochat V, Forster J, Lee S, Twardziok SO, Kanitz A, Wilm A, Holtgrewe M, Rahmann S, Nahnsen S, Köster J (2021) Sustainable data analysis with Snakemake. F1000Res 10:33

    Article  PubMed  PubMed Central  Google Scholar 

  28. Cock PJA, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, Friedberg I, Hamelryck T, Kauff F, Wilczynski B, de Hoon MJL (2009) Biopython: freely available python tools for computational molecular biology and bioinformatics. Bioinformatics 25:1422–1423

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Huerta-Cepas J, Serra F, Bork P (2016) ETE 3: reconstruction, analysis, and visualization of phylogenomic data. Mol Biol Evol 33:1635–1638

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Yeşilbağ K, Christine F, Barbara B-W, Zeki Y, Feray A, Aykut O, Ibrahim B, Sibilina Cedillo R, Heinz-Jürgen T, Matthias K (2008) Genetic heterogeneity of bovine viral diarrhoea virus (BVDV) isolates from Turkey: Identification of a new subgroup in BVDV-1. Vet Microbiol 130:258–267

    Article  Google Scholar 

  31. Giammarioli M, Ceglie L, Rossi E, Bazzucchi M, Casciari C, Petrini S, De Mia GM (2015) Increased genetic diversity of BVDV-1: recent findings and implications thereof. Virus Genes 50:147–151

    Article  CAS  PubMed  Google Scholar 

  32. Yeşilbağ K, Förster C, Ozyiğit MO, Alpay G, Tuncer P, Thiel H-J, König M (2014) Characterisation of bovine viral diarrhoea virus (BVDV) isolates from an outbreak with haemorrhagic enteritis and severe pneumonia. Vet Microbiol 169:42–49

    Article  PubMed  Google Scholar 

  33. Oğuzoğlu TÇ, Koç BT, Coşkun N, Doğan F, Duran-Yelken S (2019) Endless variety for bovine virus diarrhea viruses: new members of a novel subgroup into pestivirus a from Turkey. Trop Anim Health Prod 51:1083–1087

    Article  PubMed  Google Scholar 

  34. Deng M, Chen N, Guidarini C, Xu Z, Zhang J, Cai L, Yuan S, Sun Y, Metcalfe L (2020) Prevalence and genetic diversity of bovine viral diarrhea virus in dairy herds of China. Vet Microbiol 242:108565

    Article  CAS  PubMed  Google Scholar 

  35. Fu L, Niu B, Zhu Z, Wu S, Li W (2012) CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28:3150–3152

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Katoh K, Standley DM (2013) MAFFT multiple sequence alignment software version 7: omprovements in performance and usability. Mol Biol Evol 30:772–780

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Nguyen L-T, Schmidt HA, von Haeseler A, Minh BQ (2015) IQ-TREE: a aast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol 32:268–274

    Article  CAS  PubMed  Google Scholar 

  38. Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Vinh LS (2018) UFBoot2: improving the ultrafast bootstrap approximation. Mol Biol Evol 35:518–522

    Article  CAS  PubMed  Google Scholar 

  39. Paradis E, Schliep K (2019) ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35:526–528

    Article  CAS  PubMed  Google Scholar 

  40. Larsson A (2014) AliView: a fast and lightweight alignment viewer and editor for large datasets. Bioinformatics 30:3276–3278

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS (2017) ModelFinder: fast model selection for accurate phylogenetic estimates. Nat Methods 14:587–589

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O (2010) New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol 59:307–321

    Article  CAS  PubMed  Google Scholar 

  43. Yu G, Smith DK, Zhu H, Guan Y, Lam TT (2017) ggtree : an r package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol 8:28–36

    Article  Google Scholar 

  44. Rambaut A, Lam TT, Max Carvalho L, Pybus OG (2016) Exploring the temporal structure of heterochronous sequences using TempEst (formerly Path-O-Gen). Virus Evol 2:007

    Article  Google Scholar 

  45. Fisher RA (1935) The design of experiments. Oliver & Boyd, Oxford

    Google Scholar 

  46. Herve M (2023) RVAideMemoire: Testing and Plotting Procedures for Biostatistics. R package version 0.9–83–3, https://CRAN.R-project.org/package=RVAideMemoire. Accessed Sept 2023.

  47. Oksanen J, Simpson GL, Blanchet FG, Kindt R, Legendre P, Minchin PR, O’Hara RB, Solymos P, Stevens MHH, Szoecs E, Wagner H, Barbour M, Bedward M, Bolker B, Borcard D, Carvalho G, Chirico M, Caceres MD, Durand S, Evangelista HBA, FitzJohn R, Friendly M, Furneaux B, Hannigan G, Hill MO, Lahti L, McGlinn D, Ouellette M-H, Cunha ER, Smith T, Stier A, Braak CJFT, Weedon J (2022) vegan Community Ecology Package. R package version 2.6–4.

  48. Wernike K, Pfaff F, Beer M (2024) “Fading out” - genomic epidemiology of the last persistently infected BVDV cattle in Germany. Front Vet Sci 10:1339248

    Article  PubMed  PubMed Central  Google Scholar 

  49. Stalder H, Bachofen C, Schweizer M, Zanoni R, Sauerländer D, Peterhans E (2018) Traces of history conserved over 600 years in the geographic distribution of genetic variants of an RNA virus: Bovine viral diarrhea virus in Switzerland. PLoS One 13:e0207604

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Hawkins GA (2017) Analysis of human genetic avriations using DNA sequencing. In: Jalali M, Saldanha FY, Jalali M (eds) Basic Science Methods for Clinical Researchers. Elsevier, Amsterdam, pp 77–98

    Chapter  Google Scholar 

  51. Stalder H, Hug C, Zanoni R, Vogt H-R, Peterhans E, Schweizer M, Bachofen C (2016) A nationwide database linking information on the hosts with sequence data of their virus strains: a useful tool for the eradication of bovine viral diarrhea (BVD) in Switzerland. Virus Res 218:49–56

    Article  CAS  PubMed  Google Scholar 

  52. Braun U, Hilbe M, Peterhans E, Schweizer M (2019) Border disease in cattle. Vet J 246:12–20

    Article  PubMed  Google Scholar 

  53. Russell GC, Grant DM, Lycett S, Bachofen C, Caldow GL, Burr PD, Davie K, Ambrose N, Gunn GJ, Zadoks RN (2017) Analysis of bovine viral diarrhoea virus: Biobank and sequence database to support eradication in Scotland. Vet Rec 180:447–447

    Article  CAS  PubMed  Google Scholar 

  54. Becher P, Orlich M, Shannon AD, Horner G, König M, Thiel HJ (1997) Phylogenetic analysis of pestiviruses from domestic and wild ruminants. J Gen Virol 78:1357–1366

    Article  CAS  PubMed  Google Scholar 

  55. Bouzalas IG, Gelasakis AI, Chassalevris T, Apostolidi ED, Pappas F, Ekateriniadou L, Boukouvala E, Zdragas A (2023) Circulation of pestiviruses in small ruminants from Greece: first molecular identification of border disease virus. Vaccines 11:918

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Robinson TP, Wint GRW, Conchedda G, Van Boeckel TP, Ercoli V, Palamara E, Cinardi G, D’Aietti L, Hay SI, Gilbert M (2014) Mapping the global distribution of livestock. PLoS One 9:e96084

    Article  PubMed  PubMed Central  Google Scholar 

  57. Luzzago C, Lauzi S, Ebranati E, Giammarioli M, Moreno A, Cannella V, Masoero L, Canelli E, Guercio A, Caruso C, Ciccozzi M, De Mia GM, Acutis PL, Zehender G, Peletto S (2014) Extended genetic diversity of bovine viral diarrhea virus and frequency of genotypes and subtypes in cattle in Italy between 1995 and 2013. Biomed Res Int 2014:147145

    Article  PubMed  PubMed Central  Google Scholar 

  58. Ebranati E, Lauzi S, Cerutti F, Caruso C, Masoero L, Moreno A, De Mia GM, Peletto S, Zehender G, Luzzago C (2018) Highlighting priority areas for bovine viral diarrhea control in Italy: a phylogeographic approach. Infect Genet Evol 58:258–268

    Article  PubMed  Google Scholar 

  59. Strong R, Errington J, Cook R, Ross-Smith N, Wakeley P, Steinbach F (2013) Increased phylogenetic diversity of bovine viral diarrhoea virus type 1 isolates in England and Wales since 2001. Vet Microbiol 162:315–320

    Article  CAS  PubMed  Google Scholar 

  60. Luzzago C, Ebranati E, Sassera D, Lo Presti A, Lauzi S, Gabanelli E, Ciccozzi M, Zehender G (2012) Spatial and temporal reconstruction of bovine viral diarrhea virus genotype 1 dispersion in Italy. Infect Genet Evol 12:324–331

    Article  PubMed  Google Scholar 

  61. FranceAgriMer (2021) Compétitivité de la filière française bovin. https://www.franceagrimer.fr/content/download/66954/document/20210625-COMPETITIVITE_CAHIER_BOVIN.pdf. Accessed March 2023.

  62. Schweizer M, Stalder H, Haslebacher A, Grisiger M, Schwermer H, Di Labio E (2021) Eradication of bovine viral diarrhoea (BVD) in cattle in Switzerland: lessons Taught by the complex biology of the virus. Front Vet Sci 8:702730

    Article  PubMed  PubMed Central  Google Scholar 

  63. Lindberg ALE (2003) Bovine viral diarrhoea virus infections and its control. A Rev Vet Q 25:1–16

    Article  CAS  Google Scholar 

  64. Wernike K, Gethmann J, Schirrmeier H, Schröder R, Conraths F, Beer M (2017) Six years (2011–2016) of mandatory nationwide bovine viral diarrhea control in Germany—a success story. Pathogens 6:50

    Article  PubMed  PubMed Central  Google Scholar 

  65. Gard JA, Givens MD, Stringfellow DA (2007) Bovine viral diarrhea virus (BVDV): Epidemiologic concerns relative to semen and embryos. Theriogenology 68:434–442

    Article  CAS  PubMed  Google Scholar 

  66. De Ruigh L, Bosch J, Brus M, Landman B, Merton J (2006) Ways to improve the biosecurity of bovine semen. Reprod Domestic Animals 41:268–274

    Article  Google Scholar 

  67. Courcoul A, Ezanno P (2010) Modelling the spread of bovine viral diarrhoea virus (BVDV) in a managed metapopulation of cattle herds. Vet Microbiol 142:119–128

    Article  PubMed  Google Scholar 

  68. Alpay G, Yeşilbağ K (2015) Serological relationships among subgroups in bovine viral diarrhea virus genotype 1 (BVDV-1). Vet Microbiol 175:1–6

    Article  PubMed  Google Scholar 

  69. Bauermann FV, Flores EF, Ridpath JF (2012) Antigenic relationships between Bovine viral diarrhea virus 1 and 2 and HoBi virus: possible impacts on diagnosis and control. J Vet Diagn Invest 24:253–261

    Article  PubMed  Google Scholar 

  70. Stalder HP, Meier Ph, Pfaffen G, Wageck-Canal C, Rüfenacht J, Schaller P, Bachofen C, Marti S, Vogt HR, Peterhans E (2005) Genetic heterogeneity of pestiviruses of ruminants in Switzerland. Prev Vet Med 72:37–41

    Article  CAS  PubMed  Google Scholar 

  71. Hornberg A, Fernández SR, Vogl C, Vilcek S, Matt M, Fink M, Köfer J, Schöpf K (2009) Genetic diversity of pestivirus isolates in cattle from Western Austria. Vet Microbiol 135:205–213

    Article  CAS  PubMed  Google Scholar 

  72. Vilcek S (1999) Genetic typing of bovine pestiviruses from England and Wales. Vet Microbiol 69:227–237

    Article  CAS  PubMed  Google Scholar 

  73. Agreste (2022) La Bretagne agricole et agroalimentaire. https://draaf.bretagne.agriculture.gouv.fr/IMG/pdf/14_essentiel_synthese_bretagne.pdf. Accessed on 3 April 2024.

  74. Rapey H, Gendron P-J, Healy S, Hiriart-Durruty M, Veny N, Miquel M, Bonestebe M, Dumont B (2018) La diversité de l’élevage de ruminants au sein des territoires. L’exemple de la région Auvergne-Rhône-Alpes. Économie rurale 365:89–102

    Article  Google Scholar 

  75. Iotti B, Valdano E, Savini L, Candeloro L, Giovannini A, Rosati S, Colizza V, Giacobini M (2019) Farm productive contexts and the dynamics of bovine viral diarrhea (BVD) transmission. Prev Vet Med 165:23–33

    Article  PubMed  Google Scholar 

  76. Qi L, Beaunée G, Arnoux S, Dutta BL, Joly A, Vergu E, Ezanno P (2019) Neighbourhood contacts and trade movements drive the regional spread of bovine viral diarrhoea virus (BVDV). Vet Res 50:30

    Article  PubMed  PubMed Central  Google Scholar 

  77. Ersbøll AK, Ersbøll BK, Houe H, Alban L, Kjeldsen AM (2010) Spatial modelling of the between-herd infection dynamics of bovine virus diarrhoea virus (BVDV) in dairy herds in Denmark. Prev Vet Med 97:83–89

    Article  PubMed  Google Scholar 

  78. Joly A, Fourichon C, Beaudeau F (2005) Description and first results of a BVDV control scheme in Brittany (western France). Prev Vet Med 72:209–213

    Article  PubMed  Google Scholar 

  79. Pebesma E (2018) Simple features for r: standardized support for spatial vector data. The R Journal 10:439

    Article  Google Scholar 

Download references

Acknowledgements

The authors thank all the actors involved in the project, without whom such an extensive sample set would not have been possible: veterinarians and farmers involved in samples collection, veterinary analysis laboratories enrolled in BVD diagnostics and GDS network (see corresponding departments listed in Additional file 1). We would particularly like to thank: GDS Ardèche, GDS Haute-Loire, GDS Loire GDS, Puys-de-Dôme, GDS Rhône and diagnostic laboratory Terana, especially Sylvain Naulot, Stéphane Aulagnon and the PCR team in Terana Loire for Auvergne-Rhone-Alpes region; GDS Bretagne (Alain Joly and Loïc Maurin) and Segilab for Bretagne region; GDS Calvados, GDS Eure, GDS Manche, and GDS Orne for Normandie region.

Funding

This work was funded by GDS France, ANSES, Animal Health department of INRAE through young researcher grant (JT), MERSI fund (CL), LABEO and Région Nouvelle-Aquitaine.

Author information

Author notes

  1. Julien Thézé and Guy Kouokam have equal contribution.

Authors and Affiliations

  1. Université Clermont Auvergne, INRAE, VetAgro Sup, UMR EPIA, Saint-Genès-Champanelle, France

    Claire Lescoat, Séverine Barry, Xavier Bailly & Julien Thézé

  2. Laboratoire LABÉO, Caen, France

    Delphine Perrotte, Élise Oden, Fabienne Benoit & Pierre-Hugues Pitel

  3. Ruminant Disease and Welfare Unit, National Reference Laboratory for BVD, ANSES, Ploufragan-Plouzané-Niort Laboratory, Niort, France

    Valentin Herbet, Marc Tabouret & Guy Kouokam

  4. Oniris, INRAE, UMR BIOEPAR, Nantes, France

    Gaël Beaunée

  5. ANSES, Sophia Antipolis Laboratory, Sophia Antipolis, France

    Véronique Duquesne

  6. French Federation of Animal Health Protection Groups (GDS France), Paris, France

    Guy Kouokam

Authors
  1. Claire Lescoat
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Delphine Perrotte
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Séverine Barry
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Élise Oden
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Valentin Herbet
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Gaël Beaunée
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Marc Tabouret
    View author publications

    You can also search for this author inPubMed Google Scholar

  8. Fabienne Benoit
    View author publications

    You can also search for this author inPubMed Google Scholar

  9. Pierre-Hugues Pitel
    View author publications

    You can also search for this author inPubMed Google Scholar

  10. Véronique Duquesne
    View author publications

    You can also search for this author inPubMed Google Scholar

  11. Xavier Bailly
    View author publications

    You can also search for this author inPubMed Google Scholar

  12. Julien Thézé
    View author publications

    You can also search for this author inPubMed Google Scholar

  13. Guy Kouokam
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

Conceptualization and validation, XB, GK, CL, and J.T. Resources, SB, GK, CL, DP. Methodology, SB. Investigation and data curation, GB, SB, VH, GK, CL, EO and DP. Formal analysis, software, and visualization, CL and DP. Writing – original draft, all authors. Writing – review and editing, GB, XB, GK, CL, and JT. Funding acquisition, GB, VD, JT, MT. Supervision and project administration, FB, XB, GK, P-HP and JT. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Julien Thézé.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Handling editor: Vincent Béringue.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1. Sample information.

Additional file 2. RNA extraction and BVDV 5’UTR sequencing details using Sanger protocol for datasets 1, 2 and 3.

Additional file 3. PCR conditions for (A) QIAGEN OneStep RT-PCR kit, (B) Q5 enzyme (New England Biolabs) and (C) Kapa Hifi enzyme in Nextera™ XT DNA Index Kit.

Additional file 4. Reference sequences used to assign a genotype.

Additional file 5. Validation of the 5’UTR high-throughput sequencing protocol.

(A) Relationship between the number of reads mapped and viral load (cycle threshold – Ct) per sample. Dots are coloured according to the sequencing run. (B) Relationship between percentage of reads mapped and viral loads (Cts) per sample. A polynomial regression was fitted for each run and the resulting R2 are displayed.

Additional file 6. Genetic diversity of BVDV-1 in France.

Maximum likelihood phylogenetic tree was inferred from 1316 French isolates and international reference strains. Tips coloured in yellow and blue denote previously published French and reference sequences, respectively. Asterisks indicate samples that are identical or almost identical to vaccine strains. The tree was midpoint rooted. Solid bars at right correspond to BVDV-1 genotype delimitations based on reference strains. Statistical support for nodes were assessed using an ultrafast bootstrap (UFBoot, 1000 replicates) and a Shimodaira–Hasegawa approximate likelihood-ratio test (SH-aLRT, 1000 replicates). Diamonds at internal nodes denote well-supported clades with UFboot ≥ 95% and SH-aLRT >= 80% as suggested by IQ-TREE documentation.

Additional file 7. Genotype and spatio-temporal distribution of publicly available BVDV-1 sequences.

(A) Location of publicly available sequences assigned to the main BVDV-1 genotypes circulating in France (1e, 1b, 1d, 1a and 1l). Each panel represent a genotype. Bar charts are coloured according to the continent of origin. The number of sequences per location is log-transformed. (B) Year of collection of publicly available sequences assigned to the main BVDV-1 genotypes circulating in France. Each panel represent a genotype. Bar charts are coloured according to the continent of origin. The number of sequences per year are log-transformed.

Additional file 8. Phylogenetic relationships between French and international sequences of BVDV-1 genotype 1d.

Maximum likelihood phylogenetic tree inferred from French and international sequences of genotype 1d. Tips are coloured according to the sample continent or country of origin and their sizes represent the number of closely related publicly available sequences (99.5% identity) from the same location. European countries with less than ten sequences were grouped under the location “Europe” in the phylogeny.

Additional file 9. Regional disparities in BVDV-1 genotype proportions.

Result of Fisher exact test (FET) post-hoc pairwise comparisons of genotype proportions (1e, 1b and other genotypes) between locations. (A) Tests were performed on all regions with a least 30 sampled sequences. AuRa and BCF codes correspond to Auvergne-Rhône-Alpes and Bourgogne-Franche-Comté regions, respectively. (B) Tests were performed on all Bourgogne-France-Comté departments with at least ten sampled sequences. Number in brackets correspond to department code.

Additional file 10. Spatial distribution and local structure of BVDV-1 genotypes in two focused French regions with reduced temporal amplitude.

Proportions of BVDV-1 genotypes at local scale in (A) Bretagne, (B) and Auvergne-Rhône-Alpes. Dataset are truncated to include only two years of sampling for each region. Pie charts denote genotype proportions per farm postal code and sizes are proportional to the number of samples. A small amount of random variation was added to each pie chart location to ensure anonymity. Grey lines outline department borders. (C) Mantel correlograms measuring correlations between phylogenetic and geographical distances for the three focused regions. Each dot is positioned in the middle of a distance class. Statistically significant correlations are indicated by a square.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lescoat, C., Perrotte, D., Barry, S. et al. Spatio-temporal distribution and international context of bovine viral diarrhoea virus genetic diversity in France. Vet Res 55, 129 (2024). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13567-024-01377-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13567-024-01377-9

Keywords

Veterinary Research

ISSN: 1297-9716