Genetic Profiles of Ten African Swine Fever Virus Strains from Outbreaks in Select Provinces of Luzon, Visayas, and Mindanao, Philippines Between 2021 and 2023

1. Introduction

Infectious diseases remain to be the primary cause of losses in food animal production. In the Philippines, recent outbreaks of African swine fever (ASF) and bird flu have been adding pressure on the animal production sector to meet the growing demand for high-value animal protein and products. Since 2005, ASFV has caused an estimated loss of at least nine million animals globally. In the Philippines, since its first detection in 2019, at least half a million animal losses were reported. Along with Indonesia, Malaysia, Thailand, and Vietnam, the Philippines is one of the five Asian countries affected by these diseases, with emerging markets in livestock production and consumption [1].

First reported in Kenya in the 1920s, ASF is a highly lethal hemorrhagic disease of domestic and wild swine, with mortality rates approaching 100% [2]. It is caused by a large, enveloped, double-stranded DNA virus of the genus Asfivirus, family Asfarviridae with a genome size of 170 to 194 kb [3,4]. The virus replicates predominantly in monocytes and macrophages of the mononuclear phagocyte system, and in the late stages of infection, other cell types may also be infected[5]. Feeding of swill or untreated food scraps to animals, and illegal movements (smuggling) of animals and pig products across boundaries are seen as the most common sources of occurrence and rapid spread of ASF[6]. ASF virus currently has a very narrow host range with no zoonotic potential, and there are no indications that this will change[7,8].

Since its re-emergence in Georgia in 2007, ASF has continued to spread across Caucasia and Eastern-Europe, reaching the Caribbean, China, and several Southeast Asian countries [9,10,11,12,13,14,15]. In the Philippines, the first outbreak was recorded in July 2019, leading to a significant reduction in the country’s swine population by around 3 million hogs. This has resulted in losses of more than PHP100 billion and an increase in the prices of pork products. As of April 2024, ASF has already spread to 17 administrative regions, with a total of 73 provinces [16]. A recent study in the Philippines pinpointed the primary factors fueling the spread of ASF. The findings highlighted swill feeding, lax farm safety protocols, and personnel movement as major contributors to the disease's transmission [17]. ASF has been considered a significant threat to worldwide pork production and currently, neither therapy nor approved vaccine exists against the disease [18]. However, several countries, including Spain, have successfully eradicated ASF. This was made possible by sufficient funding for the establishment of a network of control measures, the deployment of mobile veterinary teams, and the operation of a reference laboratory for ASF surveillance and outbreak detection [19].

Since the first report of ASF in the Philippines in 2019 and its rapid spread throughout the country, several genome sequences of locally circulating strains have already been published or uploaded to online repositories [20,21]. This study aimed to add to that collection of data and to characterize the genomes of the circulating strains of ASFV from outbreaks in select provinces of the Philippines between 2021 and 2023 by using a multi-gene-based approach. Specifically, the genomes of the field strains were characterized based on the B646L gene, CD2v serogroup, central variable region (CVR) of B602L gene (CVR_B602L), the intergenic region (IGR) between I73R and I329L genes (IGR_I73R-I329L), IGR between A179L and A137R (IGR_A179L-A137R), other regions such as O174L, K145R, Bt/Sj, J268L, ECO2, and various multigene families (MGF) such as MGF 505-5R, and MGF 505-9R and 10R regions, and were compared to other related ASFV strains. These target genes were analyzed to monitor the potential emergence of new variants and could be useful in the source tracking of strains. Genomes were also analyzed for the presence of structural variants, such as long insertion or deletion, and were also compared to representative genomes of various ASFV genotypes to infer evolutionary relatedness. To our knowledge, this is the first report on the genetic characterization of genomes of ASFV strains from outbreaks in the Philippines between 2021 and 2023 generated through a long-read sequencing platform.

2. Materials and Method

2.1. Sample Collection and Viral DNA Isolation

Whole blood samples were collected from domestic pigs exhibiting signs and symptoms of ASF by licensed veterinary consultants. The collection followed the guidelines of the Bureau of Animal Industry, Department of Agriculture (Philippines) during reported outbreaks in several locations in the Philippines (Table 1). Blood samples were stored in ice or in cool packs and sent to the BioAssets Veterinary Research and Diagnostic Laboratory (Sto. Tomas, Batangas, Philippines) for molecular diagnostics and further processing. The total DNA was extracted using the MagMax™ DNA Multi-Sample Kit (Thermo Fisher Scientific) or the IndiSpin Pathogen Kit (INDICAL Bioscience) following the manufacturer’s protocol. The quality of the DNA extracts was assessed spectrophotometrically using the DeNovix DS-11 microdrop spectrophotomer (DeNovix, USA) and visually confirmed by agarose gel electrophoresis. The extraction yield and quantity were measured using the Qubit™ 1x dsDNA HS Assay Kit (Thermo Fisher Scientific, USA) on Qubit 4.0 fluorometer (Thermo Fisher Scientific, USA). The presence of the ASF virus in the samples was confirmed either by quantitative PCR using the VetMax™ African Swine Fever Virus Detection Kit (Applied Biosystems) or the Indical virotype ASFV 2.0 PCR Kit (INDICAL Bioscience) following the manufacturer’s protocol.

2.2. Library Preparation for Targeted Tiled Amplicon Sequencing

For targeted tiled amplicon sequencing, the coding-complete genomic region was amplified directly from the extracts following a tiling amplicon long-read sequencing protocol [21] with modifications. Instead of combining the primers into pools, the primer pairs were used individually to detect primer dropouts and to ensure successful amplification of all pairs prior to library preparation and sequencing. Alternate primers were designed based on the previous coding-complete genome sequence we obtained [20] for persistent primer dropouts. Sequences from up to 100 bp upstream of the forward primer and 100 bp downstream of the reverse primer (from low-performing primer pairs) were obtained and were used in subsequent targeted amplification to fill in the gaps. Alternate primers were assessed using the NCBI Primer-blast tool [22]. The PCR reaction was carried out in a 15 ul reaction volume using 1 ul of DNA (5 to 10 ng/ul), 2x PCRBio VeriFi Hot Start master mix (PCR Biosystems Ltd., UK), 200 nM primers, and the thermocycling conditions were as follows: 95°C for one min, 40 cycles of 95°C for 15 s, 15 s annealing at 60°C, and 4 min and 40 s extension at 72°C. Final extension was set at 72°C for 10 min. Amplicons were confirmed visually by agarose gel electrophoresis (0.8% w/v agarose, 3 V per cm, for up to 1 h).

After PCR, the desired amplicons per sample were pooled and purified using QIAquick PCR Purification Kit (Qiagen, USA). The purified pooled amplicons were quantified using the Qubit 1x dsDNA HS Assay Kit (Thermo Fisher Scientific, USA). A total of 200 fmol amplicon DNA per sample was used as input for the library preparation following the amplicons by ligation (SQK-LSK109) protocol for R9.4.1 flow cells (flongle or minION) with native barcoding (EXP-NBD 104 or EXP-NBD 114) expansion (Oxford Nanopore Technologies [ONT], UK). Barcoded samples were pooled equimolarly, and final DNA libraries (50 fmol for MinION or 20 fmol for flongle flow cell) were loaded. Pools of four (4) to six (6) barcoded samples were sequenced at a time on a MinION mk1b device, and the reads were obtained after 24-h run in MinKNOW Software (v. 23.11.7; ONT). The genome of Pangasinan A4 strain (2021), on the other hand, was sequenced previously using SQK-RPB-004 [20].

2.3. Genome Assembly and Annotation

Raw reads were basecalled using the super accurate model (Dorado v. 7.3.9; ONT) and were demultiplexed in MinKNOW Software (v. 23.11.7; ONT) with default parameters. Basecalled reads per sample were combined in one fastq file, which was used as input to a snakemake pipeline called LILO (https://github.com/amandawarr/Lilo) [21]. After obtaining the output from LILO, basecalled reads were mapped against the resulting scaffolds to assess the assembly quality using minimap2 (https://github.com/lh3/minimap2) [23]. Sniffles2 [24] was used to detect structural variants (SVs) in the assembled genomes. For genomes that could not be assembled using LILO, ViralWasm-Consensus [25] was used as an alternate assembler. Genome coverage was determined using bedtools [26] (v. 2.31.1), minimap2 (v. 2.26-r1175) and samtools [27] (v. 1.19.2).

To annotate the assembled genomes, GenBank format of the coding sequences from ASFV strain Georgia 2007/1 (Accession: NC_044959.2) was obtained from NCBI Nucleotide database and compared against the assembled genomes in Genome Annotation Transfer Utility (GATU) tool [28] using default blastn and needle parameters.

2.4. Genetic Characterization and Phylogenetic Analysis

Genomes of the field strains were characterized based on sequences of 13 genes and regions, namely, B646L gene, EP402R for CD2v serogroup, central variable region (CVR) of B602L gene (CVR_B602L) [29], the intergenic region (IGR) between I73R and I329L genes (IGR_I73R-I329L) [30], IGR between A179L and A137R (IGR_A179L-A137R) [31], other regions such as CP204L [32], O174L, K145R [33], Bt/Sj and J268L [34], ECO2 (IGR between I329L and I215L and the partial I215L gene) [35], MGF 505-5R, and MGF 505 9R and 10R regions [13]. Furthermore, these were compared to sequences of several related ASFV strains. Additionally, to confirm the B646L genotype of the strains, the latest African swine fever virus p72 genotyping tool and biotyping [36,37] was employed.

Representative genomes from various genotypes obtained from NCBI Datasets (https://www.ncbi.nlm.nih.gov/datasets/) were downloaded (accessed on 01 July 2024) and were compared against the assembled genomes. The entire genome sequences were aligned using the MAFFT [38] (v. 7.520) in Unipro UGENE software [39] (v. 48.1) and maximum likelihood trees were constructed in IQTREE [40] (v. 2.3.1) using ModelFinder [41] and ultrafast bootstrap [42] with 1,000 replicates. Phylogenetic trees were visualized using TreeViewer [43] (v.2.2.0).

3. Results and Discussion

A total of ten (10) representative ASF-positive whole blood samples between 2021 and 2023 were obtained from six (6) provinces (Figure 1). Of these samples, the most recent was from Mindoro Oriental, which was a sample from one of the earliest cases in November 2023 in the province. Samples from Negros Occidental 2023 cases were obtained within a few months prior to the outbreak in Mindoro Oriental. The coding-complete genomes of the 10 strains were sequenced and assembled with lengths ranging from 187,609 bp for Bataan 2022 strain to about 189,500 bp for the rest of the strains (Table 1). Owing to the limitations of the tiled amplicon sequencing method we employed, we were unable to obtain sequences for the terminal inverted repeat regions at both ends of the genome [21]. The assembled genome of ASFV/Philippines/Pangasinan/A4/2021 (ON963982.2) is 192,265 bp in length, which is longer than the assemblies obtained from tiled amplicon sequencing. Shotgun whole-genome sequencing of ASFV-positive blood sample was performed using a different library preparation protocol, resulting to significantly lower mean coverage (21x) compared to tiled amplicon sequencing (1,022x up to 5,576x). The increased coverage of the tiled amplicons are said to produce more accurate assembly than shotgun sequencing directly from extracted DNA [21]. All the 10 coding-complete genomes were classified as Biotype 2 together with Georgia 2007/1 and other 120 genotype II isolates [37].

A maximum-likelihood phylogenomic tree was constructed to further analyze the genome of ASFVs in the Philippines from the period covered. The tree revealed that the ASFV genomes from the Philippines clustered together with p72 genotype II ASFVs, which include strains detected in Asia and in Europe (Figure 2).

The complete sequences of genes and genomic regions selected for a multi-gene-based characterization of the Philippine strains were obtained from annotation and through blastn homology search with Georgia 2007/1 as the reference strain (GenBank Accession No.: FR682468.2). The complete B646L and EP402R sequences of the 10 Philippine strains obtained from the assembled genomes were 100% identical to the corresponding regions of Georgia 2007/1. All the strains belonged to p72 genotype II and CD2v serogroup 8 (Figure 3 and Figure 4, and Table 2). Among the Eurasian countries, p72 genotype II is the most common circulating genotype [14]. Based on several local ASFV surveillance studies using p72 gene as target, all of the strains so far detected in the Philippines were p72 genotype II [17,44,45]. In South Korea and in many Asian countries, the most frequently detected p72 genotype II isolates were classified as CD2v serogroup 8 [14,15]. Furthermore, the strains were characterized as belonging to CVR1 Georgia variant type based on the CVR regions in the B602L gene (100% identical to Georgia 2007/1 strain).

Amino acid sequence alignment of the tetrameric tandem repeat sequences (TRS) of the B602L CVR of the Philippine strains showed a 10 tandem amino acid repeat sequence pattern (BNDBNDBNAA), which is 100% identical to the CVR of Georgia 2007/1 (Figure 5). There is no report yet of other CVR types and variants in the country and in other Asian countries [14].

The ten local strains were IGR_I73R-I329L II variants containing an additional TRS (TATATAGGAA) pattern (Figure 6). Based on the repetition number of TRS, IGR_I73R-I329L variants can be classified as IGR I (two copies), IGR II (three copies), IGR III (four copies), and IGR IV (five copies) [30]. In Europe and in Asia, IGR II is the most common genotype. Among the neighboring Asian countries between 2018 and 2023, several IGR I variants were reported in China in 2018, in Vietnam in 2019, and in South Korea in 2019 and in 2023. IGR III variants were likewise detected between 2019 and 2021 in China, Vietnam, and in South Korea. IGR IV variants were only detected in Vietnam in 2021 [14]. IGR_I73R-I329L was considered as a genetic marker for p72 genotype II intragenotypic strain discrimination and was applied in source tracking and tracing of ASFV strains in Eastern Europe [30,33].

For IGR_A179L-A137R, all 10 strains had the same number of TRS (two repetitions of ‘GATACAATTGT’) as in Georgia 2007/1 strain. There was no deletion detected but a C-to-T substitution were present at the 143rd and 144th positions in all four (4) Negros Occidental strains and in the Mindoro Oriental strain (Figure 7). Interestingly, these strains were all from 2023 outbreaks and were reported just several weeks apart. It is possible that these strains came from a common source.

In terms of IGR_{MGF 505 9R/10R} analysis (Table 2), all 10 strains were identical to Georgia 2007/1 and were classified as MGF-1 based on the insertion of a 17-nt TRS (GATAGTAGTTCAGTTAA) [46]. Recent expansion of MGF variants were based on the number and type of TRS found near the 9R and 10R genes as variation of the 17-nt TRS ‘AGTAGTTCAGTTAAGAT’ and ‘AGTTCATTTAAGTCAAT’, respectively. Among the IGR_{MGF 505 9R/10R} variants, MGF-1 (with ABBCD__EFGHHH pattern) is the largest group, comprising almost 90% of strains from all sampled countries in one report. Other IGR_{MGF 505 9R/10R} variants (MGF-2 to MGF-8) have been detected in European countries such as in Russia, Romania, Lithuania, Latvia, and in Poland [35].

The O174L gene can also be used in tandem with other genes or regions in strain tracing and source tracking [35]. Variants of this gene are divided into three types: variant I, which is 100% identical to Georgia 2007/1; variant I with SNP; and variant II, with a 14-nt TRS (CAGTAGTGATTTT) insertion. All the 10 Philippine strains belong to O174L variant I. Furthermore, K145R and MGF 505-5R genotyping showed that all the Philippine strains are 100% identical to Georgia 2007/1 and are therefore considered as variant I (Table 2). Variant II would have C65167A in K145R and G38332A in MGF 505-5R [33].

The region spanning IGR_I329L-I215L and partial I215L gene, which is named as ECO2, has been used in ASFV strain tracking in Eastern Europe [35]. ECO2 variants can be grouped into four (4) variants: ECO2-I variant is 100% identical to Georgia 2007/1, while ECO2-II, ECO2-III, and ECOII-IV variants have the SNP at the 62nd position in I215L gene region, as A498G, and as G446A, respectively. All 10 Philippine strains were ECO2-I variants (Table 2). In Europe and in Asia, most of the ASFV strains sampled were ECO2-I variants, while several ECO2-III and ECO2-IV variants detected recently in China [35].

The intergenic Bt/Sj region of all 10 Philippine strains was identical to Georgia 2007/1, while the J268L and CP204L genes of all but the Batangas 2021 strain (Philippines/BTG2021KSU1-1/2021) were 100% homologous to the reference strain (Table 2). The Batangas 2021 strain has G-to-A substitution at the 144th position of J268L gene and an A-to-G substitution at the 210th position of CP204L gene relative to the Georgia 2007/1 sequence.

The genomes of the 10 Philippine strains were also analyzed for presence of structural variant (short or long indels). A 1.9-kb deletion was detected in the 5’-end of the genome of Bataan 2022 strain (Philippines/BAN20221-4/2022), which corresponded to the region between 17,000th and 18,939th positions of the Georgia 2007/1 genome (Figure 8). The deletion affected the MGF 360-6L, MGF 360-4L and ASFV G ACD 00300 genes.

Members of the MGF 360 and MGF 530 were observed to be important in virus cell tropism and might be required for efficient virus replication in macrophages [47,48]. In the field isolate OURT88/3 and in the tissue-culture adapted BA71V isolates, deletion of 6 members of MGF 360 and up to 2 members of MGF 530 resulted to attenuation [49]. On the other hand, the function or role of ASFV G ACD 00300 gene is yet to be elucidated. The animal source of the samples from Bataan province showed acute clinical symptoms suggesting that the strain that infected the animals was not attenuated despite having a deletion that affected MGF 360-6L and MGF 360-4L.

The analysis of ASFV whole-genome sequence of strains remains to be the most accurate method for tracing the source of the virus and in understanding its spatiotemporal evolution. However, this method can be time-consuming and can be limited by the availability of comparable strain data. Our in-house tiling amplicon long-read sequencing approach allowed us to obtain coding-complete ASFV genomes that can be subjected to a multi-gene based characterization. Analyzing select gene markers allows for source tracking, strain tracing, and variants detection, and enables differentiation between ASFV strains to better understand their introduction and spread. As shown in our data, the circulating genotype II ASFV strains in the Philippines have already shown evidence of evolution into several variants depending on the genetic markers analyzed. Further study should focus on understanding various ASFV subgroups that could already be present in the Philippines through analysis of additional genetic markers and generation of more whole-genome sequences of strains from other outbreaks. These information on strains from select provinces in the country are timely, especially now that the Philippines has authorized the use of modified live ASFV vaccines.