Molecular Tracking of The Origin of Vesicular Stomatitis Outbreaks in 2004 and 2018, Ecuador

1. Introduction

Vesicular stomatitis virus (VSV) is an RNA arbovirus [1] belonging to the Rhabdoviridae family and the Vesiculovirus genus [2] and causes an acute, contagious vesicular disease. To date, two VSV serotypes have been identified: Indiana (IN) and New Jersey (NJ) [3]. Of these two serotypes, VSIV has three subtypes: VSIV 1 (Indiana 1), VSIV 3 or VSV-AV (Indiana 3 or Alagoas), and VSIV 2 (Indiana 2 or Cocal) [4]. The serotypes that affect North, Central, and South America the most are VSNJV and VSIV 1 [5], although VSIV 3 prevails in Brazil, and VSIV 2 has been reported in Brazil and Argentina [6].

The usual hosts of VSV are generally mammals, such as cattle, horses, and pigs, but also arthropods and occasionally humans. This viral infection can be considered a minor zoonosis [7]. Although vesicular stomatitis can be transmitted by direct contact and fomites, insects are believed to play a major role in transmission between animals [8]. The insects that carry VSV are known to include sandflies, black flies, and biting midges. The ability of VSV to be transmitted through vectors is a great problem, since these vectors do not function only as simple means of transport, but also as a virus reservoir. Vector control, therefore, should form an integral part of controlling vesicular stomatitis outbreaks [6,9]. In Ecuador, vesicular stomatitis is a disease whose control is mainly based on quarantines and animal movement restrictions. Occasionally, a vaccine is used to try to reduce viral spread, namely a bivalent vaccine manufactured by VECOL, Colombia, which stimulates the production of specific antibodies [10].

Vesicular stomatitis virus causes huge economic losses in cattle because of a considerable decrease in milk and meat production [11]. In addition, because its clinical manifestations are extremely similar to those of foot-and-mouth disease, outbreaks of VS result in closing borders, livestock travel restrictions, and quarantine measures [12]. Concerning the clinical signs of VS, after an initial febrile period, characteristic vesicular lesions occur on the muzzle, lips, tongue, ears, sheath, udder, ventral abdomen, and/or coronary bands of the hoof [13]. The incubation period varies from 2–6 days with an average of 3–5 days in cattle, and visible vesicular lesions can appear 24 h after infection. The first sign of illness is often excessive salivation. In humans, the incubation period varies from 24 h to 6 days [14].

Between January and December 2018, an outbreak of VS occurred in Ecuador with a total of 399 confirmed cases throughout 18 provinces. This outbreak has been described by Salinas et al. in an earlier publication [15]. The main objective of their study was to characterize the temporal and spatial dynamics of the prevalence of VS. Their study determined that the highest number of VS cases happened between weeks 35 and 44. Concerning the etiological agent, 93.23% of the cases were positive for VSV-NJ, 5.52% for VSIV, and 1.25% for both serotypes. Similarly, it was noted that the outbreaks emerged at the end of the dry season and the beginning of the rainy season. It is therefore believed that this outbreak was caused by the movement of vectors. However, the origin of the outbreak remains unknown and the rapid spread of the outbreak within a few weeks throughout the whole country over several hundreds of kilometers cannot be attributed to transmission by vectors only. This study aimed to trace the origin of the outbreak and to determine transmission patterns through the construction of phylogenetic trees of the 2018 isolates and strains from an outbreak in Ecuador in 2004 [16]. For this purpose, we used the nucleotide sequences of the hypervariable region of the phosphoprotein P of the virus that has been used extensively for VSV-NJ and VSIV phylogenetic analysis.

2. Materials and Methods

2.1. Virus Isolation in Cell Culture

To obtain enough viral material, a cell culture of the isolates, previously stored in Gibco™ DMEM (Thermo Fisher Scientific, USA) culture medium and stored in liquid nitrogen, was initiated. The viruses were grown following the protocol established by the OIE Manual of Diagnostic Tests and Vaccines for Terrestrial Animals 2021, Chapter 3.1.23: Vesicular stomatitis (Version adopted in May 2021) [17]. Briefly, the swabs with the VSV strains, stored in liquid nitrogen, were inoculated into cell culture vessels with BHK-21 cells and incubated at 37°C for 1 h [absorption]. After washing the cell cultures three times with cell culture medium, the cell culture medium was replaced with medium containing 2.5% fetal bovine serum (FBS) and incubated at 35–37°C and examined for cytopathic effect (CPE). At 24 h post-infection, a clarified medium of infected cells was frozen and used for RNA isolation. No additional strains were isolated after 48 h of incubation.

2.2. RNA Extraction and Reverse Transcriptase-PCR [RT-PCR]

For viral RNA extraction, the AccuPrep^® Viral RNA Extraction Kit (Bioneer, South Korea) was used following the protocol established by the manufacturer. To carry out the RT-PCR, two pairs of primers that amplified the VSV phosphoprotein gene (protein P) were used. The primers for VSNJV spanned a region of the gene between position 102 (P-NJ 102; 5’ GAGAGGATAAATATCTCC 3’) and 744 (P-NJ 744; 5’ GGGCATACTGAAGAATA 3’) and resulted in a fragment of 642 base pairs. The VSIV primers amplified a region between position 179 (P-Ind 179; 5’ GCAGATGATTCTGACAC 3’) and 793 (P-Ind 793; 5’ GACTCT(C/T)GCCTG(A/ G) TTGTA 3’) and resulted in a fragment of 614 bp [17].

For the RT-PCR reaction, the MegaFi™ One-Step RT-PCR kit (ABM) was used, consisting of the following concentrations for a final volume of 15 μL: 0.3 μM of each primer, 1 ng/μL of RNA, 1X of One-Step RT-PCR Buffer, and 1U/μL RT-PCR Enzyme Mix. For reverse transcription, cDNA synthesis was performed at 60°C for 15 min. The amplification procedure started with an initial denaturation at 98°C for 30 s, followed by 35 cycles at 98°C for 10 s, 52°C for 30 s, 72°C for 30 s, and a final extension of 2 min at 72°C.

RT-PCR products were analyzed by electrophoresis using 1.5% agarose gel, with SYBR Safe for visualization and 1X TBE as the buffer. Electrophoresis was run at 100 V for 30 min in a Labnet Enduro Gel XL horizontal chamber (Labnet International, Inc., USA). The gel was visualized on a Chem-iDocTM Imaging Systems (BioRad, USA) photographic documentation system and using Image LabTM software (BioRad, USA). The positive fragments were sequenced with Sanger sequencing an ABI 3500xL Genetic Analyzer (Applied Biosystems, USA) by BigDye 3.1^® capillary electrophoresis matrix.

2.3. Phylogenetic Analysis

2.3.1. Construction of Phylogenetic Trees

For the generation of a consensus sequence of the P proteins of the VSNJV and VSIV isolates, forward and reverse sequences were aligned using the UniPro Gene program. Sequence contig reconstruction was performed using Assembler by MacVector software 17.5.5. [18]. The sequence identity was confirmed by BLAST in NCBI resources. A total of 68 sequences of VSV from the United States, Mexico, Honduras, Panama, Costa Rica, Nicaragua, Guatemala, Argentina, Brazil, El Salvador, Colombia, and Ecuador were retrieved from GenBank (All accession numbers, sample origin, serotypes, and isolates are described in Supplementary Table S1) [16,19,20,21,22,23,24,25,26,27,28,29]. For the VSNJV analysis, we included as an outgroup two sequences of VSIV: an Indiana strain from Ecuador and one from Costa Rica with accession numbers EF028138.1 and EF028137.1 respectively. For VSIV analysis two of our VSNJV sequences (ON567174 and ON567120) were included as an outgroup. DNA sequences were aligned using MacVector 17.5.5 by the ClustalW algorithm with high gap creation and extension penalties by 30.0 and 10.0, respectively, searching for a strong positional homology.

Finally, the phylogenetic trees of VSNJV and VSIV isolates from Ecuador were built from the 2004 and 2018 isolates by using the maximum parsimony (MP) (100 replicates, reweighted under Rescaled Consistency index, RC index) using PAUP 4.0a [30] and maximum likelihood (ML) analyses conducted in MEGA X [31]. The ML tree was corrected with the Tamura 3-parameter model [32]. The final trees for VSNJV and VSIV were edited using FigTree 1.4.4 [33]. The robustness values for all the analyses were estimated using bootstrapping with 500 pseudo-replicates and are shown as percentages. A bootstrap re-sample with 1000 pseudo-replicates did not show significant differences in the values of the resulting topologies. The strain number of the 2018 samples, the GenBank accession number, the geographic origin in Ecuador, the month of isolation, and the serotype can be found in Supplementary Table S2 together with the accession numbers of the sequence of the P gene from the outbreak in 2004 in Ecuador. We also built phylogram trees for VSNJV and VSIV with branch lengths proportional to the amount of evolutionary divergence among clades (Supplementary Figures S1 and S2).

2.3.2. Haplotype Network of the VSNJV Sequences

The number and distribution of haplotypes were determined using the DnaSP software [34]. The haplotype network was built using the PopART software package [35] with the TCS network inference method, using the statistical parsimony method based on Templeton et al. [36] and Clement et al. [37] and implemented in TCS to visualize the mutation relationships. This approach identifies the number of haplotypes in the data and the number of evolutionary steps taken from one haplotype to another. TCS is useful in disentangling relationships at the intraspecific level. Finally, Tajima’s D was computed with the same software.

3. Results

From the 399 stored bovine epithelial tissue samples, 67 virus isolates (64 VSNJV and 3 VSIV strains) were successfully propagated in cell culture, amplified with RT-PCR, and sequenced for the P protein. Our sequences were uploaded to GenBank under the accession numbers ON567111-ON567177. See also supplementary Table S2. Forward and reverse sequences were aligned using the UniPro Gene program and analyzed with the bioinformatic program MEGA-X together with the sequences of the P protein of 11 strains isolated from an outbreak in 2004, and described and deposited in GenBank by Sepúlveda et al. [16]. For this purpose, the initial alignment matrix consisted of 614 bp and was trimmed to obtain a final alignment of 410 bp VSIV matrix. The VSVNJ sequences in the alignment matrix were trimmed to 422pb from an initial alignment of 642pb. The New Jersey and Indiana strains were separately analyzed.

3.1. New-Jersey Phylogenetic Tree

The phylogenetic tree based on the nucleotide sequences of phosphoprotein gene (P) built for our 2018 and 2004 isolates from Ecuador and sequences from the Americas (Figure 1) differentiated the strains into three supported clades (A, B, and C in Figure 1). The internal monophyletic group, with Indiana strains from Ecuador and Costa Rica as outgroups, showed Clade A with North American strains (USA and Mexico), then Clade B with Mexican and Central American strains (B1 and B2 respectively), and the internal Clade C with South American strains including the Colombian strains and two internal subclades composed of Ecuadorian isolates of 2004 that were isolated in 6 different provinces (Supplementary Table S1) and three internal Ecuadorian subclades (Ecuador-Andean/Coast, Ecuador-Andean and Ecuador-Amazonian) with our strains of 2018.

The 2018 epizootic subclade is formed by three lineages: the Amazonian and Andean with strong bootstrapping support (75 and 100% respectively) and a small number of unsupported Andean and Coastal strains.

The Amazonian strains come from the provinces of Napo, Morona Santiago, Pastaza, Zamora Chinchipe, Orellana, and Sucumbios. The strains of the Andean clade were isolated in the provinces of Imbabura, Pichincha, Carchi, and Cotopaxi, and the coastal lowlands strains were isolated in the provinces of Guayas, Los Ríos, Santo Domingo de Los Tsáchilas, and Esmeraldas. See also the supplementary Figure 1 with a phylogram showing the VSNJV strains used in this publication.

3.2. Indiana Phylogenetic Tree

Figure 2 shows the phylogeny of VSIV strains from Ecuador in comparison to additional sequences available from the Americas. Again, strains are region-specific in topology. The basal Clade A is comprised of northern and southern Brazil, the inner Clade B is formed by sequences from Argentina and southern Brazil, and Clade C by sequences from North, Central, and northern South America. Clade A corresponds to VSIV 2, Clade B to VSIV 3, and Clade C to VSIV 1. The strains from Ecuador are paraphyletic and have two different close ancestors to Clade C. The sequences from the Amazonian 2004 and 2018 (Loja and Orellana Provinces) strains have a Colombian ancestor, while the Andean-Coastal sequences (Pichincha and Santo Domingo de Los Tsáchilas) are located basally in Clade C and are closely related to the Brazilian and Argentinian Clade B. See also the Supplementary Figure S2 with a phylogram showing a much lower divergence between the VSIV strains in comparison with the VSNJV strains shown in Supplementary Figure S1.

3.3. Haplotype Network Topology of VSV-NJ

The parsimony haplotype network (Figure 3) shows a structured haplotype geographic distribution. Low down in the network, the same phylogenetic tree sequence transformation series is observed from haplotypes of North America, and Central America, finishing with Colombian haplotypes to the left side. One of them (H4) is shared with the Ecuadorian sequence of the 2004 outbreak, with two additional unique sequences from Ecuador: H6 and H8 with 3 and 6 mutations respectively. Then, high above in the network, the Colombian H3, which connects with the main Ecuadorian subnetwork as follows: 5, 7, and 10 mutations with Ecuador’s 2004 Amazonian outbreak (center in purple, from Orellana and Napo), Ecuador’s 2018 Amazonian outbreak (left in magenta) and 2018 Andean outbreak (right in yellow) respectively. Two unique Coastal haplotypes (H32 and H29 in orange) emerge from one mutation each from the most frequent haplotype in the Andes (H30). The 2018 outbreak haplotypes H27-Amazonian (left, magenta) and H36-Sierra (right, yellow) are closely and directly connected to the intermediate branch of the 2004 outbreak (H9, purple) in Ecuador using 3 and 4 mutations respectively.

4. Discussion

Vesicular stomatitis is known to be an endemic disease in North, Central, and South America, and outbreaks of the disease in Ecuador occur frequently. Between 2005 and 2014, Ecuador reported 182 outbreaks of VS to OIE-WAHIS (World Organization for Animal Health – Animal Health Information System). At that time, VS was a notifiable disease. However, due to the mild, self-limiting nature of the disease and improbable international spread through the trade of animals, VS was de-listed by the OIE as a reportable animal disease. Nonetheless, small VSV outbreaks have been reported every year in Ecuador, with a notable increase in 2018 with 399 registered outbreaks. Three factors likely favored the rapid spread of the disease in that year: the presence of vectors, climatic conditions, and, although self-quarantine measures were taken, the transport of cattle to auctions or between farms in the Andes highlands. In the next paragraphs, we discuss what are the likely causes (mechanisms) underlying transmission.

4.1. VSNJV Epidemiology

Concerning the VSNJV isolates, we show that two clades were present during the 2018 outbreak and these clades are geographically distinguishable from each other: one with three lineages in the Amazon and the other that infects animals on the Pacific coast and in the Andean highlands. This last clade consists of strains that are all closely related. The presence of two clades is possibly related to the fact that these regions are geographically separated by the Andes mountains that run from north to south, dividing the Ecuadorian mainland into three regions, from the west to the east: the Pacific coast, the Andean highlands, and the Amazon basin (see Figure 4). Due to this geographic barrier, different virus populations have diverged, and thus different transmission patterns [35].

Regarding the transmission of VSV in 2018, the analysis showed several small and independent outbreaks in the Amazon region and the spread of another outbreak caused by a homogeneous genetic variant in the eastern coastlands and the Andes. It is well known a large number of stomatitis infections are caused by vectors [3,39,40] and the presence of several small, genetically variable, and frequent outbreaks each year in the Amazon region [41] indicates transmission from different hosts, most probably insects or perhaps mammals [42,43,44] and not due to contact between cows. Cow density is the lowest in this part of Ecuador, and most dairy farms are small-scale with only a few cows and little contact between cattle farms. We hypothesize that the outbreak in the eastern coastal lowlands and the Andes is most likely transmitted by the movement of livestock. The outbreak in this area took place over a relatively short time of 3 months (September-November see supplementary Table S2) and spread over several hundreds of kilometers. Cow density is high in this region and cattle trade is much more common. The highest density of cattle in Ecuador is found in the coastal and highland regions with 42.4% and 48.4% of livestock respectively, while the remaining 9.2 % is found in the Amazon [45]. If insects were involved in this outbreak, the genetic variation of the virus would have been more diverse due to the genetic adaptation of VSV to different insect reservoirs [40]. Moreover, insects are rare, especially in the Andes region with a much colder climate than the Amazon region, thus separating the Amazon region from the Andes. A low density of cattle in the Amazon means fewer contacts between animals, thus transmission is mainly via insect bites. In the coastal region, cattle trade and movement are more frequent, increasing the risk of the virus infection being transmitted between animals and facilitating a rapid spread of the virus over hundreds of kilometers.

Concerning the origin of VSV in Ecuador, a comparison of the VSNJV sequences of Ecuador and other American countries shows that strains from Colombia and Ecuador have a high genetic homology, with more similarity in the provinces close to the border between the two countries. The basal position of the Colombian sequences in the phylogenetic tree provides some evidence that the virus, at some time before 2004, entered Ecuador from Colombia. The network topology suggests a Colombian origin for the 2004 outbreak (H4 shared, and H6, H8 closely connected), and an Ecuadorian 2004 variant, as a close ancestor of the 2018 outbreak. This is shown in the intermediate branch of H9 and H7 Amazonian strains that connect with the 2018 Ecuadorian subnetworks See Figure 3, Of course, the haplotype network alone is not indicative of enzootic transmission in the Amazon. Many other factors such as geographical isolation and vector species play a role. The complete list of haplotypes and the countries where they were isolated is available in Supplementary Table S3.

4.2. VSIV Epidemiology

Concerning the epidemiology and transmission of VSIV strains, only limited information can be deducted from our analysis. We sequenced the P gene of only 3 strains from the 2018 outbreak. Only one other strain from Ecuador, from the outbreak in 2004, is available in GenBank. The phylogenetic tree topology of our Amazonian sequences is closely related to a Colombia sequence. The Andean sequences, basal in the clade, are derivates from Brazilian-Argentinian sequences. This suggests two different origins for the virus from Ecuador that emerged simultaneously in 2018, but also an unknown enzootic transmission cycle as a precursor of the epizootics.

The origin of the strains closely related to the Colombia sequence can be explained by the close economical relationship between the two countries (cattle trade) and sharing of insect reservoirs. However, the simultaneous occurrence of two lineages in the 2018 outbreak in the Andean mountains in the South-Central-North America clade requires further studies.

4.3. Transmission by Vectors

Concerning the transmission of VSV by vectors, multiple studies have shown that different arthropods – such as Lutzomyia sand flies, Simulium black flies, and Culicoides biting midges – are competent vectors in the transmission of VSV-NJ. Suspect vectors with a very low transmission probability include Aedes mosquitoes. [42,43]. The virus is also maintained in sandflies by transovarial transmission and in Culicoides by venereal transmission [44,45,46]. Moreover, experimental studies have shown that cows exposed to the bites of these arthropods are capable of generating antibodies against VSV [46,47].

In Ecuador, VSNJV was isolated in 1976 in the mosquito Mansonia indubitans [47,48], which breeds in lagoons with aquatic floating plants very common in the South American lowlands [53]. This supports the idea of several insect reservoirs present in the Ecuadorian Amazon causing regular, although not always detected, enzootic transmission, with eventual emerging epizootics as in the 2004 and 2018 Amazonian outbreaks. Most probably several genetic variants of the virus exist in the Amazon. The outbreaks took place in counties several hundred kilometers separated indicating that different gene pools of the virus could exist. For Influenza A it has been shown that host and geographic barriers shape the competition, extinction, and coexistence patterns of H1N1 lineages and clades. [56] Also, like the Venezuelan equine encephalitis virus, a different and enzootic vector can be present in the Amazon or the interface between forest and agriculture [49,50]. It has been demonstrated that the genetic diversity of the West Nile virus depended on the mosquito species [57]. If different hosts for the virus exist in the Amazon these hosts most probably will transmit different lineages.

Therefore, systematic surveillance of potential sylvan vertebrates and insect vectors should be performed in those areas close to the epizootic areas to elucidate the mechanisms of emergence and re-emergence of VSV. The virus could also have spread thanks to the interaction with vertebrate animals, such as bats, deer, monkeys, and herbivorous rodents, which tend to serve as amplifiers since in these organisms sustainable levels of viremia can occur [42,43,48,49,50].

5. Conclusions

The two clades in the VSNJV phylogenetic tree of the 2018 outbreaks suggest a different transmission pattern in the Andean and Amazon regions. In the Amazon, the genetic diversity of the sequenced virus suggests transmission by vectors as a possible source of infection. In this region vector control measurements, to limit vector exposure and dispersion, are required to control viral spread. Vector control is very important for the control of virus spread as has been shown for the management of the spread of VSV in equines in the USA [57]

In the Andean region, a group of genetically identical viruses was found, suggesting transmission of the virus through cattle movement. Because of the relatively cold climate in this region and the low density of insects, transmission from other reservoirs plays a most probably minor role. A future clonal spread of the virus can be avoided when strict cattle movement restrictions are established at the first detection of infected animals.

6. Recommendations

South American countries lack surveillance of VSV vectors and knowledge of reservoirs and vectors. We recommend carrying out vector and sylvatic vertebrate reservoir research and surveillance in tropical and subtropical areas where wild reservoirs can be a source of new outbreaks.

In Ecuador, VS is a disease whose control is mainly based on self-quarantine. Also, a commercial bivalent vaccine (VECOL, Bogota, Colombia) is occasionally used in Ecuador to reduce viral spread. However, no efficacy studies are available for this vaccine. Therefore, studies should be planned to show the efficacy of this vaccine in overcoming the spread of the virus. Although there is limited data on VECOL vaccine efficiency, other studies have shown protective immunity in cattle vaccinated on a commercial scale with inactivated, bivalent vesicular stomatitis vaccines [54]. It is also important to address how the vaccination status of neighboring countries might contribute to VSV mutations and spread in Ecuador. However, there is no public information about vaccination frequency in neighboring countries.

7. Limitations of the Study

Our analysis is based on 67 strains, only a small fraction of the 399 registered outbreaks in 2018. It is also possible that many more outbreaks took place in that year that was not reported. And thus, our sample is not a fair representation of the research population, possibly introducing a bias in our results.

To obtain enough viral material for this study, a cell culture of the VSV isolates, previously stored in a culture medium and stored in liquid nitrogen, was initiated. Cell culture propagation of the virus can introduce mutations that might not reflect the ancestral state of the original isolates. We examined the literature and it has been shown that, that only four single-nucleotide substitutions in 77,500 bases were found in the VSV genes after cell propagation (a mutation frequency  of  5.16×10⁻⁵) [55], and our study is based on the sequencing of only about 27,470 bp. More indications that no mutations were introduced during cell propagation are the sequences of the strains circulating in the Andes. All 35 VSV strains, isolated on different time points and again grown in cell culture, amplified, and sequenced, had an identical sequence of the partially amplified P protein. We, therefore, believe that cell propagation did not introduce a bias in our study.