Prevalence and Molecular Characteristics of Bovine Respiratory Syncytial Virus in Beef Cattle in China

Simple Summary Bovine respiratory syncytial virus (BRSV) is an important pathogen causing cattle respiratory disease; however, the prevalence and molecular characteristics of BRSV in China remain largely unknown. The purpose of this study was to investigate the prevalence and molecular characteristics of BRSV in beef cattle with BRDC in China, and the results showed that BRSV had a wide geographical distribution, and subgroup III strains were the dominant strains in China. The Chinese strains in this study showed a unique evolutionary trend based on phylogenetic analysis of the G and F genes and genomic sequences, which contributed to a better understanding of the prevalence and genetic evolution of BRSV. Abstract Bovine respiratory syncytial virus (BRSV) is an important pathogen of the bovine respiratory disease complex (BRDC); however, its prevalence and molecular characteristics in China remain largely unknown. In this study, 788 nasal swabs from 51 beef cattle farms with BRDC outbreaks in 16 provinces and one municipality were collected from October 2020 to July 2022, and 18.65% (147/788) of samples from 23 farms across 11 provinces were detected as BRSV-positive by reverse transcription-insulated isothermal PCR (RT-iiPCR) assay. Further, 18 complete G gene sequences were classified into BRSV subgroup III, and 25 complete F gene sequences were obtained from 8 and 10 provinces. Compared to the known BRSV strains in GenBank, the G proteins and F proteins in this study shared several identical amino acid (aa) mutations. Moreover, five nearly complete genome sequences were obtained and clustered into a large branch with two America BRSV subgroup III strains (KU159366 and OM328114) rather than the sole Chinese strain (MT861050) but were located in an independent small branch. In conclusion, this study reveals that BRSV has a wide geographical distribution in China, and subgroup III strains, which have unique evolution characteristics, are the dominant strains. The results contribute to a better understanding of the prevalence and genetic evolution of BRSV.


Introduction
Bovine respiratory syncytial virus (BRSV) belongs to Orthopneumovirus, of the family Pneumoviridae, and is an important pathogen of bovine respiratory disease complex (BRDC) [1]. According to the phylogenetic analysis of G gene nucleotide sequences, BRSV can be divided into subgroups I-X. A previous study showed that different subgroup strains had antigenic differences and obvious geographical distribution characteristics [2][3][4]: subgroup I has been detected in Britain and Switzerland [5]; subgroup II has been detected in Belgium, France, Denmark, Sweden, Japan, and the Netherlands [6]; subgroup III has been detected in America, China, Turkey, and Brazil, etc. [7][8][9]; subgroup IV has been detected in Germany, Belgium, Denmark, America, and other European countries [10]; subgroups V and VI have been detected in France and Belgium [11]; subgroups VII and VIII have

RNA Extraction and cDNA Synthesis
The nasal swabs were diluted 1:5 (w/v) with phosphate-buffered saline (PBS) an centrifuged at 10,000× g for 8 min, and then filtered through 0.22 μm mesh. RNA wa extracted from 350 μL of the nasal swab suspension using RNAios Plus (TaKaRa Bio, Inc Kusatsu, Japan), according to the manufacturer's instructions. The cDNA was synthesize using the PrimeScript RT Reagent Kit, according to the manufacturer's instruction (TaKaRa Bio, Inc., Kusatsu, Japan), and stored at −20 °C.

Amplification of Complete G and F Gene Sequences from Clinical Samples
A pair of primers (G-F: 5′-GAACCATCAACCAATCAAGT-3′; G-R: 5

RNA Extraction and cDNA Synthesis
The nasal swabs were diluted 1:5 (w/v) with phosphate-buffered saline (PBS) and centrifuged at 10,000× g for 8 min, and then filtered through 0.22 μm mesh. RNA was extracted from 350 μL of the nasal swab suspension using RNAios Plus (TaKaRa Bio, Inc., Kusatsu, Japan), according to the manufacturer's instructions. The cDNA was synthesized using the PrimeScript RT Reagent Kit, according to the manufacturer's instructions (TaKaRa Bio, Inc., Kusatsu, Japan), and stored at −20 °C.

Amplification of Complete G and F Gene Sequences from Clinical Samples
A pair of primers (G-F: 5′-GAACCATCAACCAATCAAGT-3′; G-R: 5′-CGCCATCCTTATTTGCC-3′) was designed for the amplification of the 1086bp sequence located at nt 4476-5562 in the reference genome USII/S1 (KU159366), which contains the complete G gene sequences. Another pair of primers (F-F: 5′-AGAGAGCACCAA-Indicates the province or municipality of the sample collected in this study,

RNA Extraction and cDNA Synthesis
The nasal swabs were diluted 1:5 (w/v) with phosphate-buffered saline (PBS) and centrifuged at 10,000× g for 8 min, and then filtered through 0.22 μm mesh. RNA was extracted from 350 μL of the nasal swab suspension using RNAios Plus (TaKaRa Bio, Inc., Kusatsu, Japan), according to the manufacturer's instructions. The cDNA was synthesized using the PrimeScript RT Reagent Kit, according to the manufacturer's instructions (TaKaRa Bio, Inc., Kusatsu, Japan), and stored at −20 °C.

Amplification of Complete G and F Gene Sequences from Clinical Samples
A pair of primers (G-F: 5′-GAACCATCAACCAATCAAGT-3′; G-R: 5′-CGCCATCCTTATTTGCC-3′) was designed for the amplification of the 1086bp sequence located at nt 4476-5562 in the reference genome USII/S1 (KU159366), which contains the complete G gene sequences. Another pair of primers (F-F: 5′-AGAGAGCACCAA-Indicates provinces with BRSV positive-sample distribution.

RNA Extraction and cDNA Synthesis
The nasal swabs were diluted 1:5 (w/v) with phosphate-buffered saline (PBS) and centrifuged at 10,000× g for 8 min, and then filtered through 0.22 µm mesh. RNA was extracted from 350 µL of the nasal swab suspension using RNAios Plus (TaKaRa Bio, Inc., Kusatsu, Japan), according to the manufacturer's instructions. The cDNA was synthesized using the PrimeScript RT Reagent Kit, according to the manufacturer's instructions (TaKaRa Bio, Inc., Kusatsu, Japan), and stored at −20 • C.

Amplification of Complete G and F Gene Sequences from Clinical Samples
A pair of primers (G-F: 5 -GAACCATCAACCAATCAAGT-3 ; G-R: 5 -CGCCATCCTTA TTTGCC-3 ) was designed for the amplification of the 1086bp sequence located at nt 4476-5562 in the reference genome USII/S1 (KU159366), which contains the complete G gene sequences. Another pair of primers (F-F: 5 -AGAGAGCACCAAGCAGAGC-3 ; F-R: 5 -CATATTTGCAGGGATTTCTTC-3 ) was designed for the amplification of the 2273bp sequence located at nt 5263-7526 in the reference genome USII/S1 (KU159366), which contains the complete F gene sequences. Primers were synthesized by Tsingke Biotechnology Co. (Chengdu, China). All amplification products were purified and cloned into the pMD19-T simple vector (TaKaRa Bio, Inc., Kusatsu, Japan) and sequenced Tsingke Biotechnology Co. (Chengdu, China) in both directions.

Genome Amplification
Ten pairs of primers ( Table 2) were designed for the amplification of the genome of BRSV according to the complete genome sequence of BRSV in GenBank. Primers were synthesized by Tsingke Biotechnology Co. (Chengdu, China). All amplification products were purified and cloned into the pMD19-T simple vector (TaKaRa Bio, Inc., Kusatsu, Japan) and sequenced Tsingke Biotechnology Co. (Chengdu, China) in both directions. The nucleotide sequence was assembled using SeqMan software (Version 7.0, DNA Star, Madison, WI, USA). Putative ORFs and their corresponding amino acids were predicted using the ORF Finder tool (http://www.ncbi.nlm.nih.gov/gorf/gorf.html), accessed on 12 September 2022.

Detection of BRSV
Among 788 clinical samples, 147 samples (18.65%; 0-37%) were detected as BRSVpositive by RT-iiPCR; the positive samples were distributed across 23 cattle farms in 11 provinces. Detailed information on the detection results in different provinces is shown in Table 1 and Figure 1.

Molecular Characterization of the G Gene Sequences
In total, 18 complete G gene sequences (GenBank accession number: OP137017-OP137034) were amplified from clinical samples from 18 farms in eight provinces (Sichuan, Inner Mongolia, Ningxia, Gansu, Shanxi, Henan, Hebei, and Qinghai). The 18 G gene sequences were 792 bp in length, encoding 264 amino acids, and shared 96.6-100% nt identity (94.7-100% aa identity) with each other; moreover, they shared 80.7-95.9% nt identity (69.1-97.1% aa identity) with the known subgroup III strains in GenBank. Interestingly, the 18 G proteins shared seven identical aa mutations, compared with all 55 known subgroup III G proteins in the GenBank ( Table 3). In addition, compared with the sole Chinese subgroup III strain (DQ strain, MT861050) in GenBank, the 18 G proteins, which were identical in length to some subgroup III strains in GenBank, had 18 continuous nucleotide inserts (nt 5465-5483 in the reference genome USII/S1, KU159366), resulting in 6 amino acid inserts and another 21 identical aa mutations. The details of the amino acid mutations in the strains considered in this study are shown in Table 4.
A neighbor-joining phylogenetic tree ( Figure 2) based on all available complete G nucleotide sequences of BRSV in GenBank indicated that the 18 strains in this study clustered into a large branch with two American subgroup III strains (USII/S1, KU159366, and BRSV\KS\467\2021, OM328114), rather than the sole Chinese strain (DQ strain, subgroup III, MT861050), but were located on an independent small branch.
A neighbor-joining phylogenetic tree ( Figure 3) based on all the available complete F nucleotide sequences of BRSV in GenBank indicated that the 25 complete F gene sequences clustered into a large branch with two known American subgroup III strains (USII/S1, KU159366, and BRSV\KS\467\2021, OM328114), rather than the sole Chinese strain (DQ strain, subgroup III, MT861050), but were located in an independent small branch.

Genomic Characterization of BRSV Strains
The five nearly complete genome sequences (GenBank accession number: OP137030-OP137034) were successfully amplified from clinical samples obtained in five provinces  (Table 5). With the exception of the identical aa mutations in the G and F proteins described above, the five strains also shared four identical aa mutations, the P protein (A74S), SH protein (N81A), and L protein (F1247C and C954Y), compared to thirteen known complete genome sequences in GenBank (Figure 2).
A neighbor-joining phylogenetic tree ( Figure 4) based on all the available nucleotide sequences of the BRSV complete genome in GenBank indicated that five strains clustered into a large branch with two known American subgroups III strains (USII/S1, KU159366, and BRSV\KS\467\2021, OM328114), rather than the sole Chinese strain (DQ strain, subgroup III, MT861050), but were located in an independent small branch. No recombination event was found in the five nearly complete genome sequences.

Genomic Characterization of BRSV Strains
The five nearly complete genome sequences (GenBank accession number: OP137030-OP137034) were successfully amplified from clinical samples obtained in five provinces  (Table 5). With the exception of the identical aa mutations in the G and F proteins described above, the five strains also shared four identical aa mutations, the P protein (A74S), SH protein (N81A), and L protein (F1247C and C954Y), compared to thirteen known complete genome sequences in Gen-Bank (Figure 2).    Note: The nucleotide homology of the 5 strains in this study with all 13 known BRSV strains in GenBank were indicated by bold and underlined, and amino acid homology was not shown.

Prevalence of BRSV in China
In recent years, beef cattle breeding has experienced rapid growth in China, and BRDC has become a major risk in beef cattle. There are multiple causative agents leading to BRDC, among which BRSV is an important pathogen of BRDC; however, the prevalence and molecular characteristics of BRSV in China remain largely unknown. In this study, 18.65% clinical samples of cattle with BRDC were detected as BRSV-positive; the positive samples were distributed in 23 cattle farms across eleven provinces, and the geographical distance between the two farthest provinces was more than 2000 km, which suggests that BRSV has a wide geographical distribution in China. Further, subgroup III strains, which are distributed across America, Turkey, Brazil, and Italy, etc. [3,7,9,23,24], were found to be the dominant BRSV strains in China in this study, and they have unique evolutionary characteristics. The subgroup Ⅲ strains have caused several outbreaks of respiratory dis-

Prevalence of BRSV in China
In recent years, beef cattle breeding has experienced rapid growth in China, and BRDC has become a major risk in beef cattle. There are multiple causative agents leading to BRDC, among which BRSV is an important pathogen of BRDC; however, the prevalence and molecular characteristics of BRSV in China remain largely unknown. In this study, 18.65% clinical samples of cattle with BRDC were detected as BRSV-positive; the positive samples were distributed in 23 cattle farms across eleven provinces, and the geographical distance between the two farthest provinces was more than 2000 km, which suggests that BRSV has a wide geographical distribution in China. Further, subgroup III strains, which are distributed across America, Turkey, Brazil, and Italy, etc. [3,7,9,23,24], were found to be the dominant BRSV strains in China in this study, and they have unique evolutionary characteristics. The subgroup III strains have caused several outbreaks of respiratory disease in America, which have brought huge economic losses to the cattle industry [7]; therefore, more attention should be paid to BRSV in China's beef cattle industry. A previous study showed that vaccination was an effective method to prevent BRSV, and there were antigenic differences among different subgroups [2,3]. However, there is still no commercial vaccine for BRSV in China; the results of this study could contribute to BRSV vaccine development in China.

Molecular Characterization of the G Gene Sequences
The G protein was mainly involved in receptor binding and the adsorption process [17] and in the production of antibodies [25]. The G protein consists of three domains: the cytoplasmic domain (1-37aa), the transmembrane domain (38-65aa), and the extracellular domain (66-257aa) [26]. The extracellular domain has a central hydrophobic region (CHR; 158-189aa), which is highly conservative, and 174-185aa in the CHR are immunedominant [27]. In this study, 18 G proteins shared seven identical aa mutations (Table 3) compared with all subgroup III G proteins in GenBank, and seven aa mutations were located in the extracellular domain. Interestingly, the 18 G proteins had an identical aa mutation (L180S) in the immune-dominant region (174-185aa) of CHR. According to the linear epitopes of the four aa (180aa,183aa, 184aa, and 205aa) of the G protein, BRSV can be divided into four antigen subgroups: A, AB, B, and unclassified [28,29]. Among them, Leu180 and Thr205 were in subgroup A, Leu180 and Ala205 were in subgroup AB, and Pro180 (Ser183 and Pro184) and Ala205 were in subgroup B. Moreover, 180aa mutation may affect the determination of the antigen subgroup of the 18 strains. The effect of the amino acid mutation of the 18 strains on G protein antigenicity needs further study.

Molecular Characterization of the F Gene Sequences
The F protein is involved in the adaptive immune response to stimulate the production of neutralizing antibodies and can promote the entry of viral particles into the host cells and mediate the fusion of infected cells to form syncytium [14][15][16]. The F protein has three hydrophobic peptides, including an amino-terminal signal peptide region (1-26aa), a site for proteolytic cleavage (131-136aa), and a hydrophobic transmembrane anchor sequence (522-549aa). Through the analysis of 25 F proteins in this study, it was found that 25 strains had five amino acid deletions due to initiation codon mutation (ATG-ACG), and an identical aa mutation (C25G), compared with other known strains in GenBank; these mutations were located in the amino-terminal signal peptide region. At present, the effect of the amino-terminal signal peptide region on the detailed function of the F protein is unclear, but the mutation of the initial codon of the F protein of Newcastle Disease Virus (ND), a member of the paramyxovirus family, will affect the virulence of the virus [30]. The effect of initial codon mutation on BRSV's virulence needs further study. In addition, the hydrophobic transmembrane anchor region (522-549aa) of the F protein involves the anchoring of the F protein to the cell membrane [31]. Interestingly, four amino acid mutations (V533A, V535A, V537M, and V544A) were found in three strains from Gansu Province in this region. Whether the amino acid mutation in this region will affect the anchoring of the F protein and host cells needs further study.

Conclusions
In conclusion, this study confirmed that BRSV has a wide geographical distribution and that subgroup III strains are the dominant strains in China. The Chinese strains considered in this study showed a unique evolutionary trend based on the phylogenetic analysis of the G, F, and genome sequences, contributing to a better understanding of the prevalence and genetic evolution of BRSV.  Institutional Review Board Statement: All tested nasal swabs samples were delivered to our laboratory for pathogenic diagnosis by the local farmers, and no other animal experiments were involved in this study.