Foot-and-Mouth Disease Virus Serotype O Exhibits Phenomenal Genetic Lineage Diversity in India during 2018–2022

In India, widespread foot-and-mouth disease (FMD) outbreaks occurred in 2021. The objective of this study was to identify genetic lineages and evaluate the antigenic relationships of FMD virus (FMDV) isolates gathered from outbreaks reported between 2019 and 2022. Our study shows that the lineages O/ME–SA/Ind2001e and the O/ME–SA/Cluster-2018 were both responsible for the FMD outbreaks on an epidemic scale during 2021. This observation is in contrast to earlier findings that suggested epidemic-scale FMD outbreaks in India are often connected to a single genetic lineage. Additionally, we report here the identification of the O/ME–SA/PanAsia-2/ANT10 sub-lineage in India for the first time, which was connected to two intermittent outbreaks in Jammu and Kashmir. The current study demonstrates that the O/ME–SA/ind2001e lineage has a strong presence outside of the Indian subcontinent. Furthermore, the O/ME–SA/Cluster-2018 was observed to have a wider geographic distribution than previously, and like the O/ME–SA/Ind2001d and O/ME–SA/Ind2001e lineages in the past, it may eventually spread outside of its geographic niche. For O/ME–SA/Ind2001e and O/ME–SA/Cluster-2018, the predicted substitution rate for the VP1 region was 6.737 × 10−3 and 8.257 × 10−3 nt substitutions per site per year, respectively. The time of the most recent common ancestor of the O/ME–SA/Ind2001e and O/ME–SA/Cluster-2018 strains suggests that the viruses possibly emerged during 2003–2011 and 2009–2017, respectively. Recent sightings of the O/ME–SA/PanAsia2/ANT10 virus in India and the O/ME–SA/Ind2001e virus in Pakistan point to possible cross-border transit of the viruses. The results of a two-dimensional viral neutralization test revealed that all of the field isolates were antigenically matched to the currently used Indian vaccine strain O INDR2/1975. These results suggest that the serotype O vaccine strain can protect against outbreaks brought on by all three circulating lineages.


Introduction
Foot-and-mouth Disease (FMD) is one of the main challenges to the livestock sector in India and other endemic nations. The causative agent, FMD virus (FMDV), has seven immunologically different serotypes: O, A, C, Asia 1, and Southern African Territories (SAT) 1, SAT 2, and SAT 3. These serotypes can be further divided into several topotypes, lineages, sub-lineages, and an ever-expanding ensemble of variants [1]. Three of the seven serotypes, including serotypes O, A, and Asia 1, are prevalent in India. Multiple serotypes, associated antigenic diversity, short-lived vaccine-induced immunity, rapid transmission, and other associated factors make the disease control difficult. In India, FMD outbreaks are primarily caused by serotype O (which accounts for approximately 90% of FMD outbreaks), which also dominates in other parts of the world [2]. Serotypes A and Asia1 on the other

Virus Isolation
During 2019-2022, the state FMD network laboratories collected a total of 3756 clinical samples from suspected FMD cases across various states and union territories (UTs) in India.
First, the samples were tested for serotype identification using in-house sandwich ELISA, and ELISA-negative samples were further tested using reverse transcription multiplex PCR. Subsequently, the samples were inoculated into a BHK-21 cell monolayer for virus isolation. Infected cell cultures were harvested after complete cytopathic effects were observed. The infected cell culture supernatants were used for vaccine matching and sequence-based studies. In addition, sequences were also generated directly from clinical materials for samples that could not be isolated in cell culture.

VP1 Sequencing
Following the manufacturer's instructions, the RNeasy Mini Kit (Qiagen, Hilden, Germany) was used to extract viral RNA from 138 infected cell culture supernatants and clinical samples. Reverse transcription was carried out using MMLV reverse transcriptase (Promega, Madison, WI, USA) and NK61 primer [10]. All the PCR amplification was performed using pfu DNA polymerase (Fermentas, Waltham, MA, USA). For the amplification of VP1, the primer combination of ARS4 [10] and NK61 were used. The details of PCR, the primer used for sequencing, and the thermal profile were similar to those previously mentioned [11]. The PCR products were purified using a QIAquick PCR Purification Kit (Qiagen, Germany), and the amplicons were sequenced using an ABI 3130 Genetic Analyzer (Applied Biosystems, Waltham, MA, USA).

Phylogenetic Analysis
In addition to the VP1 sequences of 138 strains generated in this study, we retrieved 223 Indian VP1 sequences from the GenBank and the Institute Genetic Database. The sequences were aligned using the MUSCLE tool. The mean and pairwise divergence were then computed. To further explore the differences between the genetic lineages and the Indian vaccine strain O/INDR2/1975, nucleotide and deduced amino acid sequences were compared using the BioEdit program version 7.2.5.0 [12]. To assess the evolutionary relationships among FMDV isolates, phylogenetic trees were inferred by the maximum likelihood (ML) method based on the nucleotide alignment of the VP1 sequences using the MEGA software v. 11 [13]. The ML phylogeny was produced under the Kimura 2-parameter evolution model with rate variation following a gamma distribution as determined by the model finder, and the robustness of the tree topology was assessed by bootstrap analysis with 1000 iterations.

Phylodynamic Analysis
In order to explore the evolutionary characteristics of recent FMDV isolates, 542 fulllength VP1 sequences from India and other countries were used for phylodynamic analyses. The presence of a temporal signal was examined by root-to-tip regression with Tempest v.1.5.3 software [14]. The Markov chain Monte Carlo (MCMC) method was performed in BEAST v.1.10.4 [15], and a relaxed and uncorrelated lognormal clock and exponential coalescent population prior were used to estimate the temporal phylogeny and substitution rate. Three independent runs of 200 million generations were carried out, their convergence was evaluated, and the log and tree files were then combined with the aid of Log Combiner. A maximum clade credibility (MCC) tree was summarized using Tree Annotator v1.10.4, with the burn-in option used to remove the first 10% of sampled trees, and the resulting tree was visualized by FigTree v 1.4.4. Phylogeographic analyses were performed, using an asymmetric substitution model with BSSVS options to infer asymmetric diffusion rates [16] between any pairwise location state and allowing BF calculations to verify significant diffusion rates.

Selection Pressure Analysis
Three likelihood approaches were employed to determine the positive selection pressure at certain codon sites: the single likelihood ancestor counting (SLAC) method, the fixed effects likelihood (FEL) method, and a Bayesian strategy called FUBAR. The ratio of non-synonymous (dN) to synonymous (dS) substitutions per site (ratio: dN/dS) was used to calculate the strength of selection pressure. In general, posterior probabilities > 0.9 for FUBAR and p < 0.1 for SLAC strongly imply positive selection. The Mixed Effects Model of Evolution (MEME) was used to identify the codon sites that were the subject of episodic diversifying selection. At significance levels (p 0.05), strong evidence of selection was accepted. All the analyses were carried out using the online Datamonkey web server [17].

Vaccine Matching Analysis
Monovalent bovine vaccinal serum (BVS) against Indian serotype O vaccine strain O/INDR2/1975 was obtained from the serum repository of ICAR-NIFMD, Bhubaneswar, India. Before testing, the serum was inactivated at 56 • C for 30 min in a water bath. Vaccine matching was performed using a two-dimensional virus neutralization test as described by [18]. The antibody titer was determined as the reciprocal of the last dilution of serum that neutralized 100 TCID50 in 50% of the wells. The relationship value (r1-value) was calculated as a ratio of antibody titers against field isolates to those against the vaccine strain, averaged from the two separate runs. An adequate antigenic homology between a field isolate and the vaccination strain is indicated by an r1-value of ≥0.3. On the other hand, the r1-value of less than 0.3 indicates an antigenic divergence.

Serotype Detection and Virus Isolation
During 2019-2022, a total of 532 FMD outbreaks were serotype confirmed, of which 489 were caused by serotype O, accounting for about 92% of the total FMD occurrences in the country (Table 1). Compared to 2019 and 2020, approximately a six-fold increase in the number of outbreaks was observed in 2021. The outbreaks were extensively reported from several states and UTs during 2021. Apparently, the surge in the number of outbreaks in 2021 was due to serotype O, which was responsible for 92% of the total FMD outbreaks reported. Moreover, serotype O was found to be responsible for 98%, 83%, and 93% of the FMD outbreaks, respectively, in 2019, 2020, and 2022. Serotypes A and Asia1 were associated with sporadic incidences. In total, 3756 clinical samples were processed for serotype identification using sandwich ELISA and RT-mPCR, which revealed serotype O in the majority of the FMD-positive samples (n = 1502). The clinical materials were passaged in BHK-21 cells, and FMDV could be isolated from 190 samples, of which 165 were confirmed to be serotype O.

Phylogenetic Relationships
Earlier, we reported the genetic characterization of FMDV serotype O isolates collected during 2014-2018 from India (5). In the current study, 13 isolates collected during 2018 have also been sequenced, in addition to those sequenced during 2019-2022. In total, 138 isolates sequenced in this study (Supplementary File S1) and 223 Indian sequences retrieved from the public domain and institute data bank were used for comparative analyses. Phylogenetic analysis based on the full-length VP1 region showed that all FMDV strains collected in India between 2018 and 2022 could be divided into three distinct lineages, O/ME-SA/Ind2001e, O/ME-SA/Cluster-2018, and O/ME-SA/PanAsia-2, with supporting bootstrap values of 99%, 84%, and 99%, respectively (Figures 1and S1). Out of 138 strains sequenced in this study, 78 clustered within the O/ME-SA/Ind2001e lineage, whereas 56 isolates grouped within the O/ME-SA/Cluster-2018. Interestingly, four isolates grouped together within the O/ME-SA/PanAsia-2, sharing descent with ANT10 sub-lineage, which has never been identified in India so far.
These PanAsia-2 strains were isolated from two outbreaks reported from Jammu and Kashmir in 2021, shared ancestry, and demonstrated 96% sequence homology with an isolate collected from Pakistan in 2019 ( Figure 2). Unfortunately, sequences collected after 2019 from the region are not available in the public domain for inclusion in the comparison. The inclusion of a larger number of sequences of recent origin from this group will shed more light on the possible transmission pattern. The O/ME-SA/PanAsia-2/ANT10 showed pairwise mean genetic distances of 10.1 and 12.4% from O/ME-SA/Cluster-2018 and O/ME-SA/Ind2001e, respectively, at the nucleotide level. The O/ME-SA/Cluster-2018 and the O/ME-SA/Ind2001e differed by 12.6% in the mean distance. The O/ME-SA/Cluster-2018 shared ancestry with the O/ME-SA/PanAsia-2. With the currently used Indian serotype O vaccine strain INDR2/1975, the isolates of these three genetic groups showed a mean genetic distance of 12.1 to 13.5%.  FMDV isolates (n = 138) sequenced in this study were collected from 79 outbreaks (on many occasions, more than one virus isolate from the same outbreak was sequenced). FMD strains associated with 64 outbreaks during 2021 were characterized phylogenetically, which revealed the involvement of the O/ME-SA/Ind2001e lineage in 36 outbreaks, the O/ME-SA/cluster-2018 in 26 outbreaks, and two outbreaks were caused by O/ME-SA/PanAsia-2 ( Figure 3). Though it was not possible to sequence all the FMDV outbreak strains due to the non-receipt of samples or inappropriate quality of samples leading to failure in PCR amplification, it can be presumed that the FMD epidemic observed in 2021 was due to both the O/ME-SA/Ind2001e lineage (56 percent of the outbreaks) and the O/ME-SA/cluster-2018 (41 percent of the outbreaks). On some occasions, from a single outbreak, both O/ME-SA/Ind2001e and O/ME-SA/Cluster 2018 lineage could be detected. For instance, an outbreak in the states of Karnataka (January 2021), Maharashtra (August 2021), and Jammu and Kashmir (July 2021) was caused by the simultaneous involvement of both lineages. FMDV isolates (n = 138) sequenced in this study were collected from 79 outbreaks (on many occasions, more than one virus isolate from the same outbreak was sequenced). FMD strains associated with 64 outbreaks during 2021 were characterized phylogenetically, which revealed the involvement of the O/ME-SA/Ind2001e lineage in 36 outbreaks, the O/ME-SA/cluster-2018 in 26 outbreaks, and two outbreaks were caused by O/ME-SA/PanAsia-2 ( Figure 3). Though it was not possible to sequence all the FMDV outbreak strains due to the non-receipt of samples or inappropriate quality of samples leading to failure in PCR ampli-

Phylodynamic Analyses
In this analysis, 542 serotype O VP1 sequences, comprising 361 Indian and 181 foreign sequences obtained from GenBank (accessed on 24 January 2023), were included. India is predicted to be the ancestral root state for O/ME-SA/Ind2001e and O/ME-SA/Cluster-2018 in the current analysis with good statistical support (root state posterior probabilities = 0.99). For the O/ME-SA/PanAsia-2/ANT10 sub-lineage, Pakistan had the highest root-state posterior probabilities of 0.99 ( Figure 4). The O/ME-SA/ind2001e lineage has been reported from Bangladesh, Bhutan, Nepal, Myanmar, Thailand, the United Arab Emirates, Saudi Arabia, Russia, Mauritius, China, South Korea, and Jordan during 2015-2018 [6,19]. The O/ME-SA/Ind2001e has firmly established itself outside of the Indian subcontinent, as this study further demonstrates. Recently, outbreaks of O/ME-SA/Ind2001e lineage were reported in Pakistan in 2019 [20]. The outbreaks in Pakistan were due to three phylogenetically distinct groups closely related to strains circulating in Nepal, India, and

Lineage O/ME-SA/Ind2001e
The dataset contains 432 sequences, 288 of which were collected in India and the rest from other countries. The isolates of the O/ME-SA/Ind2001e lineage showed a maximum and mean genetic divergence of 4.7 and 3.0% at the nucleotide level and 3.8 and 2.0% at the amino acid level, respectively. Overall, in the nucleotide alignment, 339 sites were found to be invariable and 300 sites showed polymorphisms. Out of 300 polymorphic nucleotide sites, 220 were parsimony informative and 80 were singleton variable sites. The substitution rate for the VP1 region was estimated to be 6.737

Lineage O/ME-SA/Cluster-2018
The dataset comprised VP1 sequences from 73 isolates, of which 69 isolates were sampled in India and only four isolates were collected in Bangladesh in 2021. The pairwise nucleotide and amino acid divergence among the cluster 2018 was determined to be 2.3-6.8% and 1.4-3.7%, respectively, with a mean divergence of 3 and 1%. Overall, in the nucleotide alignment, 510 sites were found to be invariable and 129 sites showed polymorphisms. Out of 129 nucleotide sites, 83 were parsimony informative and 46 were singleton variable sites. In total, 99 codon positions out of 213 were found conserved. The dN/dS ratio for O/ME-SA/Cluster-2018 was estimated to be 0.127, indicating strong purifying selection and evidence for purifying selection was found at nine codon positions.

Vaccine Matching
In total, 59 serotype O field isolates were subjected to a vaccine-matching exercise. Field isolates were selected based on their geographical location, collection date, genetic group, and adaptation to the BHK-21 cell culture. The isolates selected represent the O/ME-SA/Ind2001e lineage (n = 30), O/ME-SA/Clustre-2018 (n = 27), and O/ME-SA/Pa-nAsia-2/ANT10 (n = 2). All the isolates, irrespective of the genetic groups to which they belong, showed antigenic match (r-value > 0.3) with the currently used vaccine strain INDR2/1975 (Table 2). Both the isolates of O/ME-SA/PanAsia-2/ANT10 also showed antigenic homology with the vaccine strain.

Vaccine Matching
In total, 59 serotype O field isolates were subjected to a vaccine-matching exercise. Field isolates were selected based on their geographical location, collection date, genetic group, and adaptation to the BHK-21 cell culture. The isolates selected represent the O/ME-SA/Ind2001e lineage (n = 30), O/ME-SA/Clustre-2018 (n = 27), and O/ME-SA/PanAsia-2/ANT10 (n = 2). All the isolates, irrespective of the genetic groups to which they belong, showed antigenic match (r-value > 0.3) with the currently used vaccine strain INDR2/1975 (Table 2). Both the isolates of O/ME-SA/PanAsia-2/ANT10 also showed antigenic homology with the vaccine strain.

Discussion
The prevalence of FMD in India is a major hurdle to the growth of the livestock industry due to its adverse impact on productivity, and trade in livestock and livestock products. In India, a uniform vaccine strain policy and standard vaccination strategy have been implemented countrywide. Under this program, all cattle and buffaloes are vaccinated bi-annually with an inactivated trivalent FMD vaccine for protection against FMD. Complex epidemiology of the disease poses a serious challenge to its control in FMD-endemic countries. In India, FMD is primarily caused by FMDV serotype O. FMD outbreaks on an epizootic scale were recorded in India in 2021. In this study, genetic and antigenic characterization of FMDV serotype O viruses isolated from India during 2018-2022 is reported. The VP1 coding sequence of 138 serotype O isolates and antigenic relationships of 59 isolates with the Indian vaccine strain were determined. By combining reference sequences from India and other parts of the world, the molecular epidemiology of FMDV serotype O was studied using Bayesian phylogeographic inference and the maximum likelihood method.
Our analyses revealed that all serotype O viruses collected in India belonged to the ME-SA topotype. The majority of the isolates were grouped within the O/ME-SA/Ind2001e lineage and the O/ME-SA/cluster-2018, representing 56% and 41% outbreaks, respectively. The lineage O/ME-SA/Ind2001e has circulated extensively in South East Asia, the Middle East, and East Asia since its first detection in Nepal in 2012 and has caused a series of epidemics in several countries outside of Pool 2 [6]. In India, the O/ME-SA/Ind2001e lineage was first detected in 2015, in spite of intensive surveillance [5].  [6].
The O/ME-SA/Cluster-2018 was first detected in 2018 in India in a limited number of outbreaks, and subsequently increase in circulation has been noticed in 2021 and 2022. After its first detection in the state of Maharashtra, O/ME-SA/Cluster-2018 lineage has spread to six more Indian states. Almost 40% of the FMD outbreaks observed in 2021 were caused by this lineage. Outside India, the lineage was detected in Bangladesh in 2021, which indicates that the lineage is progressively undergoing geographic expansion. The Bayesian phylogeographic analysis provided strong support for India as the country of origin for O/ME-SA/Ind2001e and O/ME-SA/cluster-2018. Nonetheless, it should be highlighted that sample biases that affect phylogeographic reconstruction quality, especially in areas with sparse data, may have an impact on the precise origins of epidemics [6].
The recent documentation of O/ME-SA/PanAsia-2/ANT10 in India and O/ME-SA/Ind2001e in Pakistan clearly indicates virus exchange between the two countries. Though transmission and incursion of FMDV were reported in countries such as Bangladesh and Nepal, the exchange of FMDV between India and Pakistan is buzzing as the two countries do not indulge in direct livestock trade [20], and also have strict border control. Since the lineage was detected in Jammu and Kashmir, a UT bordering Pakistan, the possibility of airborne transmission cannot be excluded. The inclusion of recent sequences of PanAsia-2 will give more insights into the evolutionary and transmission events. FMDV-O topotype ME-SA and lineage PanAsia-2 are the extensively distributed lineages in Iraq, Iran, Afghanistan, and Pakistan, and the dominant sub-lineage of serotype O that is commonly detected in Pakistan is O/ME-SA/PanAsia-2/ANT10 [21]. Animal movements and international trade are considered important risk factors for the emergence of various exotic FMDVs, including serotypes and lineages that are not included in vaccination plans [22]. On many occasions, virus sequences from different regions and states of India clustered together closely, indicating frequent inter-state transmission of FMDV. This chain of transmission events is mainly caused by the unregulated movement of infected animals, as reported in previous studies [4,5]. This altogether highlights the huge risk potential of the transboundary nature of FMDV and its ability to easily spread across the country.
In India, higher FMD incidences on the epidemic scale occur cyclically after every few (3)(4)(5) years, and generally, a new genetic lineage of serotype O is associated with such an upsurge. Surprisingly, the FMD epidemic observed in 2021 was caused by two lineages, including O/ME-SA/Ind2001e and O/ME-SA/Cluster-2018. Another interesting observation is that the previous three epidemics occurred in a gap of 5-6 years, but the 2021 epidemic happened within three years after the last epidemic in 2018. This could be due to the fact that during 2019 and 2020, vaccinations were not practiced due to COVID-19-linked lockdown and movement restrictions. The FMD epidemics in Sri Lanka generally appear every four to six years, although the evolutionary dynamics of the patterns of appearance remain unclear [23].
FMDV evolves primarily by point mutation, and the size of this mutant swarm has a significant impact on the virus's ability to adapt, spread, and cause disease [24]. The appearance of new virus strains can pose a challenge to disease control strategies in endemic settings, especially when the in-use vaccine strain fails to offer sufficient crossprotection [25]. In the vaccine matching analysis, all the serotype O isolates showed an r-value of >0.3, which indicates a perfect antigenic match with the field isolates. The antigenic relationship between the field strain and the vaccine strain established through vaccine matching studies is a critical determinant of the efficacy of the vaccine strain used. However, other important factors such as regularity of vaccination schedule, vaccine coverage, antigenic mass in the vaccine, cold chain logistics from the factory to field, and biosecurity breaches at the outbreak site, including unrestricted animal movement, might contribute to the number and extensiveness of outbreaks in the country. Generally, serotype O vaccines are broadly cross-reactive, often exhibiting in vitro protection against a number of genotypes. For instance, in South America, the O/Campos vaccine has been used for routine vaccination for more than five decades and still matches the viruses circulating in the region [26]. The Indian type O vaccine strain (O/IND/R2/1975) has also been reported to have antigenic similarities with viruses from East Africa, mainly Eritrea, Ethiopia, Kenya, and Sudan [27]. The strain O/INDR2/1975 has been in use for more than four decades in India and is still able to provide cross-protection to different genetic variants within the ME-SA topotype circulating in the country.

Conclusions
The increased number of outbreaks of FMD in 2021 was associated with two lineages, viz., O/ME-SA/Ind2001e and O/ME-SA/Cluster-2018. Proportionately based on the sequences analyzed, it can be presumed in simplistic terms that 60% of the FMD outbreaks might have been caused by the O/ME-SA/Ind2001e lineage, with 40% caused by the O/ME-SA/Cluster-2018. Two isolated FMD outbreaks in Jammu and Kashmir was caused by O/ME-SA/PanAsia-2/ANT10, and seems to be self-limiting. The PanAsia-2 strain could not be detected further in any of the FMD outbreaks recorded subsequently. The O/ME-SA/Cluster-2018 appears to be the next dominant lineage, as previously anticipated [5], and it may eventually spread outside of its geographic niche of origin, much like the O/ME-SA/Ind2001d and e lineages in the past. Irrespective of genetic diversity, vaccine matching analyses established clear antigenic homology between the field isolates and the currently used Indian vaccine strain.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/v15071529/s1, File S1. The history of FMDV serotype O field isolates determined in this study. The GenBank accession numbers are provided for all the individual isolates. Figure   Institutional Review Board Statement: Not applicable as no animals were handled in this study.

Informed Consent Statement: Not applicable.
Data Availability Statement: All required data are available as texts and figures in the main text of the article or in the Supplementary Materials. The sequence data sets generated during this research are publicly available at NCBI GenBank.