Epidemiological Surveillance and Mutational Pattern Analysis of Foot-and-Mouth Disease Outbreaks in Bangladesh during 2012 – 2021

,


Introduction
Foot-and-mouth disease (FMD) contributes to huge economic losses each year in Bangladesh as well as many other endemic countries. The causative agent is an antigenically diverse RNA virus, foot-and-mouth disease virus (FMDV) which possesses a rapid transmission rate affecting clovenhoofed animals such as cattle, buffalo, sheep, goats, and pigs, etc. [1,2]. The mortality rate of FMD is relatively low and the morbidity rate can reach up to 100% [3,4]. The occurrence of FMD is affected by different risk factors such as viral, host, and environmental factors [5].
FMDV is a member of the genus Aphthovirus and the family Picornaviridae [6][7][8]. The genome is about 8.5 kb and is surrounded by an icosahedral capsid [9]. Four structural proteins form the capsid, of which VP1, VP2, and VP3 are exposed outside, and VP4 is completely internal [4]. VP1 is the most variable among the other structural proteins consisting of three major antigenic sites: sites 1, 3, and 5 [10,11]. The G-H loop, B-C loop and carboxy-terminus of VP1 conform these three major antigenic sites [12,13].
The FMDV exists as seven immunologically distinct serotypes, namely A, O, C, Asia 1, and the Southern African Territories (SAT)-1, SAT-2, and SAT-3. Each serotype further contains multiple genotypes that are usually related to the geographical region of the disease occurrence [4,7]. FMDV serotypes O, A, and Asia 1 are circulating in Bangladesh. Among them, serotype O was responsible for 82% of the outbreaks in Bangladesh [14]. Serotype C was not reported after 1990. Serotype A viruses were also reported in Bangladesh simultaneously with serotype O [14][15][16]. Circulation of FMDV Asia 1 serotype in Bangladesh is not consistent but sporadic [17].
A large FMD-susceptible livestock population of 24.7 million cattle, 26.7 million goats, 3.7 million sheep, and 1.5 million buffaloes exists in Bangladesh [18], and the occurrence of several outbreaks each year brings huge losses to the economy of Bangladesh. A study predicted that financial loss due to the FMD outbreak would be Tk. 188.57 billion or US$ 2.22 billion per year in Bangladesh [19].
The replication of the virus is erroneous and its polymerase lacks proofreading activity contributing to a wide variety of subpopulations [4,7]. Previous infection or vaccination with one FMDV serotype does not confer cross-protective immunity to another. A high level of antigenic variation leads to the frequent emergence of novel FMDV strains within a serotype in Bangladesh that may reduce the efficacy of existing vaccines. In Bangladesh, the majority of FMDV vaccination program depends on imported vaccines, and the use of imported vaccines of heterologous virus strains often leads to incomplete protection or complete vaccine failure against local strains in Bangladesh. Moreover, this vaccine provides only short-time protection (∼6 months) [20,21]. Under this circumstance, it becomes necessary to formulate FMD vaccines using circulating indigenous FMDV strains that would provide increased protection against the local virus. Again, the inadequate monitoring system, unrestricted transboundary movement of animals, and unplanned vaccination programs complicate the FMD situation in Bangladesh. All these facts necessitate constant epidemiological studies to keep track of the local strains and whether there is any emergence of newer genotypes of the virus which would be crucial to design a well-planned control strategy to prevent the disease realistically.
The current study was planned to conduct comprehensive surveillance to monitor the epidemiological situation of FMD in Bangladesh during a ten-year period from 2012 to 2021 to decipher the pattern of recent outbreaks and the effect of risk factors on the prevalence of FMD cases. The mutational trend in the VP1 sequence of circulating strains was also reported in this study. The study represents knowledge on currently circulating predominant strains as well as newly emerged FMDV strains in Bangladesh which would facilitate the selection of more appropriate FMD vaccine strains for developing more effective vaccines. This study also identified groups that are more susceptible to infection based on some risk factors such as age, gender, breed, farming system, and vaccination status of the animal, and the result would help the authorities to select which group should be prioritized more during the vaccination and control program. Therefore, the findings of the study would be valuable to design an effective control plan better suited for Bangladesh based on the prevailing FMDV situation in the country.

Materials and Methods
2.1. Sample Collection. From 2012 to 2021, a total of 481 tongue or foot epithelial tissue samples were collected from the infected lesion of FMD-suspected cattle, buffalo, and pigs covering 32 different districts of Bangladesh based on the notification of the clinical history from farmers and clinical findings within the animals. All the samples were collected by registered veterinary doctors without the use of anesthesia, following quite a safe and painless process. Sample collection was carried out following the protocol approved by the Animal Experimentation Ethical Review Committee, Faculty of Biological Sciences, University of Dhaka, with the permission of the herd owner. The ethical approval number for the study protocol is Ref: 66/Biol. Sci./2018-19; Date: November 14, 2018. A questionnaire was prepared (Supplementary 1) and filled during the collection of samples which describes the history of the patients. Samples were transported from the collection site to the laboratory at 4°C within 20 hr and stored at −80°C until processing and testing. All the samples were handled and processed in biosafety level 2 laboratory facilities.
2.2. RNA Extraction and cDNA Synthesis. The homogenization of tissue samples and total RNA extraction from tissue were performed in an automated Maxwell 16 system (Promega, USA) using the Maxwell 16 total RNA purification kit (Promega, USA) according to the manufacturer's instructions. Following extraction, the RNA was reverse transcribed into complementary DNA (cDNA) using the ImProm-IITM reverse transcription system (Promega, USA).

Polymerase Chain Reaction (PCR)-Based
Amplification of VP1 and Sequencing of VP1. VP1-based PCR diagnostic assay was employed for the detection of FMDV-positive tissue samples. VP1 region of cDNA amplification was carried out using two sets of primers (VP1UF/NK61, 16F/NK61). The PCR reaction was performed using GoTaq 2× Hot Start Colorless Master Mix (Promega, USA) with either forward primer VP1UF (5′GTACTACRCSCAGTAC-3′) [22] or 16 F (5′-GAGAACTACGGWGGWGAGAC-3′) [16] and the reverse primer NK61 (5′-GACATGTCCTCCTGCATCTG-3′) [23]. The PCR products were run on 1.0% agarose gel with a 1 kb-DNA ladder (Promega, USA) for the visualization and detection of FMD-positive amplicons. Following detection, the FMDV-positive amplified PCR products were purified using the Wizard SV Gel and PCR Clean-Up System (Promega, USA) and were subjected to an automated cycle sequencing reaction using BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, USA) according to manufacturer's instructions with the same primers used in the PCR reaction. The VP1 coding sequence quality was analyzed in ABI Genetic Analyzer (Applied Biosystems, USA). Both forward and reverse sequences were assembled into a single contig using SeqMan version 7.0 (DNASTAR, Inc., Madison, WI, USA). The assembled sequences were compared with other sequences from GenBank using the basic local alignment search tool, BLAST [24], to reveal the identity of the isolated virus as well as their serotypes.
2.4. Epidemiological Analyses. The investigation of outbreaks in this study was based on the sample collected during 2012-2021 in the Microbial Genetics and Bioinformatics Laboratory, University of Dhaka [14,25,26]. We have followed a research-based data collection, and samples were collected upon communicating with outbreak information from the farmers as Bangladesh lacks systematic surveillance of FMD outbreaks. It is possible that not all outbreak information could be collected during the specified period of time. There could be a lack of samples from each and every incidence of outbreak cases. It should also be noted that the collection of the sample and reporting of the outbreak were interrupted by the pandemic situation of COVID-19 in Bangladesh during 2020-2021. Morbidity, mortality, and case fatality rates were calculated based on the population size (n = 3,580) included in this study. The influence of various risk factors such as season, age, gender, breed, and farming system on the occurrence of FMD cases was investigated and statistically verified by the chi-square (χ 2 ) test at a 5% significance level using Statistical Package for Social Science, SPSS 26.0 for Windows (SPSS Inc., Chicago, IL, USA). The results of χ 2 test were presented in Supplementary 2.
Representative VP1 sequences of Bangladeshi isolates from 2012 to 2021 from our laboratory and also from reported sequences by other researchers in Bangladesh were included in the phylogenetic study to determine the genotype based on the clade formation in MEGA11 [30]. VP1 sequences of Bangladeshi isolates included in this study are listed in Supplementary 2.
The consensus VP1 coding sequences (complete 1D region) of local FMDV isolates were aligned using the ClustalW program [31] with the related gene sequences from GenBank. Phylogenetic Neighbor-Joining [32] trees were constructed (bootstrap replicates 1,000) based on the Kimura-2 parameter model [33]. A discrete Gamma distribution of a value of 1 was used to model evolutionary rate differences among sites. Fewer than 5% alignment gaps, missing data, and ambiguous bases were allowed at any position, as all the sequences were not completely aligned on the full range.
2.6. Analysis of VP1 Amino Acid Variations. The VP1 coding sequences of FMDV isolates reported in Bangladesh were used to analyze the overall variations among all FMDV serotype O, A, and Asia 1 viruses. The amino acid was translated based on the standard genetic code after codon-based alignment with MEGA11 software [30]. In the mutational pattern analysis, sequences under Ind-2001BD2 sublineage and also from MYMBD21 sublineage were excluded because sequences reported under both of these two sublineages had no amino acid variations. 3.1.1. Risk Factor Analysis. The occurrence of FMD outbreaks varied in different seasons. It was found that 51% (36/71) of the outbreaks occurred during the rainy season (July-October), followed by 28% (20/71) outbreaks during winter (November-February) and 21% (15/71) during summer (March-June) ( Figure 3). In the χ 2 test at 5% significance level, significant seasonal influence on the occurrence of the FMD outbreaks was observed where p<0:001 was detected (Supplementary 2).

Results
Various risk factors (age, gender, breed, farming system, and vaccination status) were considered in this study to observe the effect of these factors on the FMD morbidity rate.

Distribution of FMD Serotypes over Years (2012-2021).
Among the 71 outbreaks, 85% (60 out of 71) of the outbreaks were caused by serotype O, which was the prevalent serotype, whereas serotype A and Asia 1 were responsible for 11% (8 out of 71) and 4% (4 out of 71) of the outbreaks, respectively.
From FMDV VP1 coding sequences reported in Bangladesh, it was revealed that serotype O, A, and Asia 1 was circulating during 2021-2021. Serotype O was detected each year from 2012 to 2021 [14,26]. Serotype A was also present in the past 10 years but was not in circulation in 2015, 2018, and 2021 [16,29]. Asia 1 was reported from 2012 to 2013 and then in 2018 [17,22,34]. No serotype C    case was detected in Bangladesh ( Figure 6) during the timespan.

Phylogenetic Analysis of FMDV Serotypes.
Only representative FMDV VP1 sequences were considered in the phylogenetic study for the ease of data presentation from sequences reported from our laboratory [14,16,25,[27][28][29]22] as well as from other researchers of Bangladesh [34,35].  (Figure 9). No Asia 1 serotype was found since 2018. Sporadic circulation pattern of the Asia 1 serotype in Bangladesh was also observed in previous studies [17,22]. During 2012-2021, a few novel strains were introduced in FMDV circulation of Bangladesh, which is illustrated in Figure 10.  143, 148, 190, 194, 196, and 209. In the case of Asia 1, N47S, T50V, M146L, and M211L were unique. Mutations at the 43rd and 48th residues were common to all three serotypes O, A, and Asia 1. Table 2 represents the mutational pattern for serotypes O, A, and Asia 1 observed from 2012 to 2021. Here colored boxes (yellow for serotype O, blue for serotype A and black for serotype Asia 1) indicated the presence of a particular mutation in a given year, and white box represented the absence of the particular mutation.
Among serotype O isolates, more than two amino acid substitutions were found in positions: 43,138,140,142,197,198,201

Discussion
In Bangladesh, FMD is endemic and considered a significant threat to the cattle industry. The time period of 2012-2021  was epidemiologically critical, showing changes in the circulation pattern of FMDV in Bangladesh with the emergence of multiple FMDV variants. The FMD cases included in the study were 19% fatal, with 10.4% mortality (Figure 2).
Occurrence of the highest percentage (51%) of FMD outbreaks during the rainy season might be due to relatively higher humidity and lower temperature that facilitated the survival of the virus, while higher temperature during summer might be responsible for the lower number of outbreaks. The effect of climate was also found in a previous study by Rahman et al. [41].
The morbidity rate was higher in the adult population than in calves and young cattle because older animals are more frequently get exposed to the multiple serotypes of FMDV over time at common grazing fields or markets via more contact with animals of other herds [42]. On the other hand, calves and young animals are mostly kept in home settings resulting in less contact with other animals. The relatively lower prevalence in younger animal groups (<2 years) might be due to the low frequency of contact with other herds as well as the presence of passive maternal immunity [42][43][44][45]. But mortality and fatality rates were higher in calves which might be due to the lack of immunity against multiple serotypes among calves, whereas adults become immune to multiple serotypes via frequent exposure ( Figure 4). The reason behind more susceptibility in male animals (59.4%) than females (50.8%) might be due to more exposure of male cattle in the field settings for agricultural purposes compared to females. Working in the field increases the possibility of skin injury, and thus these cattle might easily become susceptible to infection whenever they come into close contact with other infected cattle [46,47]. The higher    Transboundary and Emerging Diseases      susceptibility of cases was observed among local or indigenous breeds (55.7%) than that of crossbred cattle (53%) as in Bangladesh cross breed cattle are reared following more safety measures and regular vaccination in commercial farms for meat and milk production, whereas indigenous ones are kept in local households and extensively used for cultivation with no regular vaccination provided by the local farmers [46][47][48]. Under the intensive farming system, animals are kept in closed settings following proper animal health guidelines, whereas in semi-intensive farming, animals are allowed to graze on the field sometimes, which might have resulted in the higher prevalence of FMD in semi-intensive farms [46][47][48][49]. From the χ 2 test at a 5% significance level, it was confirmed that season, age, gender, and farming system were significant factors for the occurrence of FMD (p<0:05). But the association of breed with FMD cases was not found significant in the statistical (χ 2 ) test (p>0:05). Another important factor confirmed by the χ 2 test was the vaccination   [14,36,50]. In 2013, another novel Ind-2001BD2 sublineage was detected, but its circulation did not continue in the subsequent years [14]. Under serotype Asia 1, G-VIII lineage was circulating during 2012-2013 after 1996, and then in 2018, Asia 1 serotype emerged as a novel lineage, G-IX (BD-18) in Bangladesh (Figure 9) [17]. Since 2018, no other cases of Asia 1 serotype were reported confirming the sporadic circulation pattern of this serotype. In 2021, a new lineage SA-2018 intruded in Bangladesh and was first reported from our Microbial Genetics and Bioinformatics Laboratory at the Department of Microbiology of the University of Dhaka [26]. This lineage was first reported from India in 2018 [51]. In VP1-based phylogeny, the isolates under SA-2018 lineage (OP320455.1-OP320458.1) formed a separate clade emerging from the Indian isolates indicating an evolutionary relationship with those Indian isolates (Figure 7), which might have resulted from unrestricted cattle movement or international trade at borders with India [21,50]. The association of unregulated animal movement with the virus transmission across the country or among neighboring countries and the emergence of exotic FMDV strains were evident in previous studies [50][51][52]. In our previous study, these isolates were confirmed as a novel sublineage, MYMBD21, under SA-2018, lineage showing a 5%-6% nucleotide distance from the Indian isolates [26]. In the case of serotype A, only one lineage G-VII under ASIA topotype was detected in Bangladesh ( Figure 8). Some isolates under serotype O were uncharacterized (BD_BAU_ML1_2013, BD_BAU_ML2_2013, BD_SI_5_2013, O/BAN/BLRI/450.2/ 2018), which require detailed VP1-based analyses (Figure 7) [35].
From the amino acid substitution analyses, it was evident that most of the changes occurred in the antigenic sites, the G-H loop, B-C loop, and C-terminal loop. The RGD motif containing Arg-Gly-Asp was found to be conserved in all serotypes as it is an important recognition element in integrindependent cell adhesion processes [4,39,53].
Many single amino acid substitutions were observed in the VP1 sequences of FMDV. Amino acid residues at 43-46, 48, 142, 198, 202, and 210 positions were variable in more than one serotype (Table 1). Amino acid substitutions K45Q, T142A of serotype O, E194K of serotype A, and E202K of serotype Asia 1 sequences were also apparent in previously conducted studies [54,[55][56][57]. Frequently reported exchanges were reported in amino acid residues-142, 194, and 210 in a previous study [4], and substitutions were also detected in those positions for Bangladeshi isolates included in this study (Table 1). Amino acid residue 142 is located within the immunodominant epitope, G-H loop adjacent to RGD motif, and substitutions at this site can alter the structural conformation of the G-H loop of VP1 depending on the type of amino acid exchange [58]. Again, amino acid variations at positions-134, 143, and 158 within G-H loop; positions-196, 198, 210, and 213 within C-terminal were detected in other studies [54, 55, 57, 59, 60-62, 63, 64], and these positions were also variable among Bangladeshi isolates in this study (Table 1). Substitutions at the 210th amino acid residue detected in serotype O and Asia 1 can cause failure in the VP1-2A product [4,63,65].
Amino acid exchanges in which a negatively charged amino acid was replaced by a positively charged amino acid had occurred in amino acid residues 52, 138, 194, 202. In residue 52, negatively charged aspartic acid (D) was substituted by positively charged lysine (K). Residues 138, 194, and 202 had common amino acid substitutionsnegatively charged glutamic acid (E) was substituted by positively charged lysine (K). Another common case of amino acid exchange occurred in amino acid residue-44, 46, 140, where uncharged, polar asparagine (N) was substituted by polar, negatively charged aspartic acid (D) ( Table 1). Position 194 was found as a variable site which is very close to 195-197 residues that form one of the walls of the heparan sulfate (HS) binding site. Introduction of positive charge by E194K mutation in serotype A can affect the HS receptor binding to the virus [4,66].
Within the serotype O, Ind-2001d contains more genetically diverse viruses than the other type O sublineages found in Bangladesh [14]. Significant mutations S139G, N140A, P142T, and P158T were found between Ind-2001d and Ind-2001e sublineages which were significant for the G-H loop displacement between these two sublineages [14]. In Table 2 Serotype A isolates from 2019 to 2020 showed more variations in the B-C loop and C-terminus. Asparagine (N) is converted into negatively charged, aspartate (D) at the 44th residue located in the antigenic site 3 within the B-C loop and can affect the antigenicity [4,66]. E194K (negative to positive charge conversion) and Q198R (uncharged to positive charge conversion) mutations at the C-terminal site can affect antigenic site 2 and thereby modulate the heparan binding site as these two positions are adjacent to one of the walls of the heparan binding site [4,66,67].
Within the G-VIII lineage of Asia 1 isolates, few mutations at position-44, 202, 210 were found, among which A44T demonstrated polarity change and E202K (glutamate, E to lysine, K) referred to negative charge conversion into positive charge. Among the unique mutations of G-IX lineage, charge and polarity shifts also occurred. Uncharged alanine (A) was converted into negatively charged glutamic acid (E) at the 44th residue, and changes in polarity were observed at position 50 where polar threonine (T) was substituted by nonpolar valine (V) in the B-C loop affecting the antigenic site. These mutations might be responsible for the lineage turnover from G-VIII to G-IX.
Apart from the mutational analysis presented in this study, further serological assays are required to confirm any changes in the antigenicity of the recently circulating FMDV strains.
Several significant mutations accumulated within the circulating FMDV virus in Bangladesh from 2012 to 2021. Some mutations at significant antigenic sites contributed to the frequent emergence of novel variants within 10-year period that challenge the effectiveness of the existing vaccines and available treatment options. For the effective control of the disease, regular monitoring of FMD outbreaks, risk factor analysis, and planning to imply effective measures should be given emphasis.  [20,25]. This vaccine is not marketed yet, and the use of traditional imported vaccines is still going on causing incomplete protection, which, in turn, giving rise to the emergence of newer FMDV variants. For instance, a new sublineage, MYMBD21 of SA-2018 lineage emerged in Bangladesh in 2021. There was a considerable VP1 nucleotide (12%-13%) and amino acid sequence (5%) divergence of MYMBD21 isolate from both the current field vaccine strain, O/India/R2/75 (accession number: AF204276.1) and proposed local vaccine strain, BAN/TA/Dh-301/2016 (accession number: KY077628.1). Again, 3D analysis revealed the displacement of the critical antigenic site, G-H loop of VP1 of MYMBD21 strain from existing vaccine strains that raised the possibility of vaccine escape and confirmation of which requires further serological analyses [26]. However, the isolates (OP320455.1-OP320458.1) belonging to this novel variant were stored in the seed bank of Microbial Genetics and Bioinformatics Laboratory, Department of Microbiology, University of Dhaka for subsequent serological analysis to assess its suitability as vaccine strain and for developing an immediate vaccine, if necessary, in future. Thus, regular monitoring regarding FMD virus evolution must be carried out to mitigate the problem of vaccine failure.
This study reported overall FMD epidemiological pattern and risk factor-based analysis during 2012-2021, providing the track of circulating and emerging FMD strains, identifying high-risk groups, analyzing mutations and variable regions within the VP1 capsid protein. The findings of the study would be valuable for better understanding the current FMD situation, which would facilitate adopting effective preventive actions which will lead to achieve the implementation of the OIE/FAO prescribed Progressive Control Pathway for FMD road map in Bangladesh.

Conclusions
One country should have its own surveillance system to plan control strategies of any infectious disease based on the disease pattern prevailing in the particular region. But Bangladesh lacks constant epidemiological surveillance of FMD, which complicates the management of the disease. To investigate the FMD situation in Bangladesh, this study analyzed FMD outbreaks in 10 years timespan. Among different groups of animals, the older male cattle reared under semiintensive farming system, indigenous and nonvaccinated cattle were found to be more prone to infection. These groups should be prioritized more while adopting any vaccination or other control strategies of FMD in Bangladesh. Detection of Ind-2001e sublineage of serotype O as the predominant genotype in recent years suggests its use as an indigenous vaccine strain. There was a recent evolution of a novel variant, MYMBD21 sublineage under SA-2018 lineage in Bangladesh that presented the possibility of vaccine escape and thus was stored in the seed bank for further serological analyses. Overall, this study provides prerequisite knowledge of FMD epidemiology in Bangladesh that would facilitate designing and implementing more appropriate control and preventive approaches, particularly for Bangladesh to fight against the highly infectious disease.

Data Availability
Data supporting the conclusions of the study can be found in the supplementary materials.

Ethical Approval
Sample collection and processing followed the protocol approved by the Animal Experimentation Ethical Review Committee (AEERC), Faculty of Biological Sciences, University of Dhaka (Ref: 66