# Use of an Alignment-Free Method for the Geographical Discrimination of GTPVs Based on the GPCR Sequences

Tesfaye Rufael Chibssa, Yang Liu, Melaku Sombo, Jacqueline Kasiiti Lichoti, Janchivdorj Erdenebaatar, Bazartseren Boldbaatar, Reingard Grabherr, Tirumala Bharani K Settypalli, Francisco J Berguido, Angelika Loitsch, et al.

► To cite this version:

Tesfaye Rufael Chibssa, Yang Liu, Melaku Sombo, Jacqueline Kasiiti Lichoti, Janchivdorj Erdenebaatar, et al. Use of an Alignment-Free Method for the Geographical Discrimination of GTPVs Based on the GPCR Sequences. Microorganisms, 2021, 9 (4), pp.855. 10.3390/microorganisms9040855. [https://doi.org/10.33069/microorganisms.9040855](https://doi.org/10.33069/microorganisms.9040855)

HAL Id: hal-03306901

https://hal.inrae.fr/hal-03306901

Submitted on 29 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Use of an Alignment-Free Method for the Geographical Discrimination of GTPVs Based on the GPCR Sequences

Tesfaye Rufael Chibssa<sup>1,2,3</sup>, Yang Liu<sup>1,●</sup>, Melaku Sombo<sup>2</sup>, Jacqueline Kasiiti Lichoti<sup>4</sup>, Janchivdorj Erdenebaatar<sup>5</sup>, Bazartseren Boldbaatar<sup>5</sup>, Reingard Grabherr<sup>3,●</sup>, Tirumala Bharani K. Settypalli<sup>1</sup>, Francisco J. Berguido<sup>1</sup>, Angelika Loitsch<sup>6</sup>, Delesa Damena<sup>2</sup>, Giovanni Cattoli<sup>1</sup>, Adama Diallo<sup>7,8</sup> and Charles Euloge Lamien<sup>1,*<sup>●</sup></sup>

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Chibssa, T.R.; Liu, Y.; Sombo, M.; Lichoti, J.K.; Erdenebaatar, J.; Boldbaatar, B.; Grabherr, R.; Settypalli, T.B.K.; Berguido, F.J.; Loitsch, A.; et al. Use of an Alignment-Free Method for the Geographical Discrimination of GTPVs Based on the GPCR Sequences. *Microorganisms* **2021**, *9*, 855. https://doi.org/10.3390/microorganisms9040855

Academic Editor: Shawn Babiuk

Received: 29 March 2021

Accepted: 14 April 2021

Published: 16 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** Goatpox virus (GTPV) belongs to the genus Capripoxvirus, together with sheeppox virus (SPPV) and lumpy skin disease virus (LSDV). GTPV primarily affects sheep, goats and some wild ruminants. Although GTPV is only present in Africa and Asia, the recent spread of LSDV in Europe and Asia shows capripoxviruses could escape their traditional geographical regions to cause severe outbreaks in new areas. Therefore, it is crucial to develop effective source tracing of capripoxvirus infections. Earlier, conventional phylogenetic methods, based on limited samples, identified three different nucleotide sequence profiles in the G-protein-coupled chemokine receptor (GPCR) gene of GTPVs. However, this method did not differentiate GTPV strains by their geographical origins. We have sequenced the GPCR gene of additional GTPVs and analyzed them with publicly available sequences, using conventional alignment-based methods and an alignment-free approach exploiting k-mer frequencies. Using the alignment-free method, we can now classify GTPVs based on their geographical origin: African GTPVs and Asian GTPVs, which further split into Western and Central Asian (WCA) GTPVs and Eastern and Southern Asian (ESA) GTPVs. This approach will help determine the source of introduction in GTPV emergence in disease-free regions and detect the importation of additional strains in disease-endemic areas.

**Keywords:** GTPV; alignment-free algorithm; k-mer; G-protein-coupled chemokine receptor

## 1. Introduction

Goatpox virus (GTPV) is a DNA virus belonging to the genus Capripoxvirus within the family Poxviridae, together with sheeppox virus (SPPV) and lumpy skin disease virus (LSDV) [1]. GTPV affects goats, sheep and small wild ruminants [2,3].

<footnote n="1">1 Animal Production and Health Laboratory, Joint FAO/IAEA Division of Nuclear Techniques in Food and Agriculture, Department of Nuclear Sciences and Applications, International Atomic Energy Agency, Friedenstrasse 1, A-2444 Seibersdorf, Austria; chibssasafo@gmail.com (T.R.C.); lydialiu329@hotmail.com (Y.L.); T.B.K.Settypalli@iaea.org (T.B.K.S.); F.Berguido@iaea.org (F.J.B.); G.Cattoli@iaea.org (G.C.)</footnote>

<footnote n="2">2 National Animal Health Diagnostic and Investigation Centre (NAHDIC), Sebeta P.O. Box 04, Ethiopia; sombomelaku@yahoo.com (M.S.); delesa_damenaa@yahoo.com (D.D.)</footnote>

<footnote n="3">3 Institute of Biotechnology, University of Natural Resources and Life Sciences (BOKU), Muthgasse 18, 1190 Vienna, Austria; reingard.grabherr@boku.ac.at</footnote>

<footnote n="4">4 Directorate of Veterinary Services, Ministry of Agriculture, Livestock and Fisheries, Private Bag, Nairobi 00625, Kenya; kasiiti.orengo@gmail.com</footnote>

<footnote n="5">5 Institute of Veterinary Medicine, Mongolian University of Life Sciences, Zaisan, Ulaanbaatar 17024, Mongolia; erdene64j@gmail.com (J.E.); boldoomglvet@yahoo.com (B.B.)</footnote>

<footnote n="6">6 Institute for Veterinary Disease Control, Austrian Agency for Health and Food Safety (AGES), Robert Koch-Gasse 17, 2340 Modling, Austria; angelika.loitsch@ages.at</footnote>

<footnote n="7">7 UMR CIRAD INRA, Animal, Santé, Territoires, Risques et Ecosystèmes (ASTRE), CEDEX 05, 34398 Montpellier, France; adama.diallo@cirad.fr</footnote>

<footnote n="8">8 Laboratoire National d'Elevage et de Recherches Vétérinaires, Institut Sénégalais de Recherches Agricoles (ISRA), Hann, Dakar BP 2057, Senegal</footnote>

<footnote>* Correspondence: C.Lamien@iaea.org; Tel.: +43-1-2600-28314</footnote>
Although goatpox is a disease of goats, the causative agent can be GTPV or SPPV. Similarly, both GTPV and SPPV can cause sheeppox and, unfortunately, it is not possible to differentiate infections of small ruminants by GTPV or SPPV either clinically or serologically [4] due to cross-reactivity between capripoxviruses (CaPVs). Hence, molecular methods are increasingly being used to genotype and classify CaPVs and undertake molecular epidemiological studies [2,5,6].

Together, GTPV and SPPV cause significant economic losses in small ruminants in endemic regions of Africa and Asia and together they make up one of the major impediments to the improvement of goat and sheep productions in the affected areas [7].

Though GTPV is currently confined to Africa and Asia, recent outbreaks of LSDV and SPPV [7] showed that CaPVs could escape their traditional geographical confinements into new areas and cause severe outbreaks. During an incursion of GTPV into a new geographical area, it would be of great importance to identify the source of introduction to implement proper control measures. The identification of a suitable molecular target for discriminating GTPV strains based on their geographical origins can complement field investigations to determine a potential route of introduction during GTPV outbreaks in disease-free areas. This can also help detect the introduction of new strains in disease-endemic areas.

Previous studies showed that the G-protein-coupled chemokine receptor (GPCR gene), a host range gene, was suitable to classify CaPVs into one of the three species by real-time PCR [5] and by phylogenetic reconstruction [2,8–11]. Using the GPCR gene, Le Goff et al. (2009) [2] could distinguish three GTPV sub-clusters comprising isolates from different geographical locations, without being able to define the geographical boundaries clearly. The limited number of samples and the absence of isolates from some geographical locations made it impossible to draw definite conclusions.

The GPCR later became a widespread target for molecular characterization of CaPVs [8–11]. Thus, a large number of sequences were produced by several research groups and made available in public databases.

Here, we have sequenced the GPCR genes of 12 GTPVs and analyzed them together with 37 GTPV sequences available in public databases.

Previous analyses using the GPCR relied on conventional sequence alignments and statistical methods for phylogenetic reconstructions. However, those methods usually exclude gap positions or perform poorly when gaps are present in the sequence alignments [12–14]. In many cases, these gaps could represent essential information to enable the classification of viruses. For instance, the GPCR sequence alignments showed the presence of a 21-nucleotide deletion in all SPPVs and several GTPVs, as well as a 12-nucleotide deletion in a few isolates of LSDV [2]. The removal of gaps during analysis can result in the loss of information that is essential for comparison.

To avoid this loss of information, alignment-free methods were proposed, including those using the k-mer frequencies to compute the pairwise distances between sequences [13–15].

Here, we have used the k-mer frequencies-based distance calculation and heatmaps to compare the GPCR sequences of GTPVs from various geographical locations and classify them according to their geographical origins.

## 2. Materials and Methods

### 2.1. Sequencing of New Isolates

Twelve isolates (Table 1) of GTPVs from Ethiopia, Kenya, Ghana and Mongolia were acquired and sequenced for this study. Total DNA was extracted from clinical samples and cell culture supernatants and the GPCR gene was amplified as previously described [8]. All pathological samples and viral isolates were handled within the biosafety level 3 facility of the Institute for Veterinary Disease Control, Austrian Agency for Health and Food Safety, Austria.

The positive PCR products were purified and sequenced at LGC Genomics (Berlin, Germany). The sequence data were assembled using Vector NTI 11.5 software (Invitrogen,
Waltham, MA, USA). The generated GPCR sequences have been deposited in GenBank under accession numbers MN161836 to MN161847.

<table><caption>**Table 1.** Goatpox viruses used in the comparative analysis.</caption><thead><tr><th>Isolate</th><th>Year of Collection</th><th>Species</th><th>Country</th><th>Region</th><th>Accession Number</th></tr></thead><tbody><tr><td>* Awi/O13/2011</td><td>2011</td><td>Sheep</td><td>Ethiopia</td><td>Africa</td><td>MN161836</td></tr><tr><td>* Bale/O14/2007</td><td>2007</td><td>Sheep</td><td>Ethiopia</td><td>Africa</td><td>MN161837</td></tr><tr><td>* Giner/O15/2007</td><td>2007</td><td>Goat</td><td>Ethiopia</td><td>Africa</td><td>MN161838</td></tr><tr><td>* Towele/O17/2013</td><td>2013</td><td>Goat</td><td>Ethiopia</td><td>Africa</td><td>MN161839</td></tr><tr><td>* Halasya/G18/2013</td><td>2013</td><td>Goat</td><td>Ethiopia</td><td>Africa</td><td>MN161840</td></tr><tr><td>* Kitengela/O58/2011</td><td>2011</td><td>Sheep</td><td>Kenya</td><td>Africa</td><td>MN161841</td></tr><tr><td>* Kitengela/O59/2011</td><td>2011</td><td>Goat</td><td>Kenya</td><td>Africa</td><td>MN161842</td></tr><tr><td>* Kiambu/G143/2009</td><td>2009</td><td>Goat</td><td>Kenya</td><td>Africa</td><td>MN161843</td></tr><tr><td>* Ghana</td><td>Unknown</td><td>Goat</td><td>Ghana</td><td>Africa</td><td>MN161844</td></tr><tr><td>Assosa/G01/2010</td><td>2010</td><td>Goat</td><td>Ethiopia</td><td>Africa</td><td>KP663696</td></tr><tr><td>Chagni/O06/2012</td><td>2012</td><td>Sheep</td><td>Ethiopia</td><td>Africa</td><td>KP663700</td></tr><tr><td>Chagni/G03/2012</td><td>2012</td><td>Goat</td><td>Ethiopia</td><td>Africa</td><td>KP663699</td></tr><tr><td>Chagni/G02/2012</td><td>2012</td><td>Goat</td><td>Ethiopia</td><td>Africa</td><td>KP663698</td></tr><tr><td>Chagni/G01/2012</td><td>2012</td><td>Goat</td><td>Ethiopia</td><td>Africa</td><td>KP663697</td></tr><tr><td>NVI/G01/2009)</td><td>2009</td><td>Goat</td><td>Ethiopia</td><td>Africa</td><td>KP663709</td></tr><tr><td>Yemen/83</td><td>1983</td><td>Goat</td><td>Yemen</td><td>West Asia</td><td>FJ869362</td></tr><tr><td>Metekel/O01/2010</td><td>2010</td><td>Sheep</td><td>Ethiopia</td><td>Africa</td><td>KP663705</td></tr><tr><td>Akaki/O01/2008</td><td>2008</td><td>Sheep</td><td>Ethiopia</td><td>Africa</td><td>KP663692</td></tr><tr><td>Isiolo_SP</td><td>1959</td><td>Sheep</td><td>Kenya</td><td>Africa</td><td>KJ818279</td></tr><tr><td>Chad-VC6</td><td>Unknown</td><td>Goat</td><td>Chad</td><td>Africa</td><td>FJ869392</td></tr><tr><td>Burkina/Benogo/3A</td><td>Unknown</td><td>Goat</td><td>Burkina Faso</td><td>Africa</td><td>FJ869353</td></tr><tr><td>DesseI</td><td>Unknown</td><td>Goat</td><td>Unknown</td><td>Africa</td><td>FJ869354</td></tr><tr><td>Chad/VC8</td><td>Unknown</td><td>Goat</td><td>Chad</td><td>Africa</td><td>FJ869363</td></tr><tr><td>Kedong_SP</td><td>1955</td><td>Sheep</td><td>Kenya</td><td>Africa</td><td>KJ818280</td></tr><tr><td>G20-LKV</td><td>2000</td><td>Goat</td><td>Kazakhstan</td><td>Central Asia</td><td>AY077836</td></tr><tr><td>SA2/2017</td><td>2017</td><td>Goat</td><td>Saudi Arabia</td><td>West Asia</td><td>MG232389</td></tr><tr><td>Iraq/61/Gorgan</td><td>1961</td><td>Goat</td><td>Iraq</td><td>West Asia</td><td>FJ869357</td></tr><tr><td>Pellor</td><td>2000</td><td>Goat</td><td>Kazakhstan</td><td>Central Asia</td><td>AY077835</td></tr><tr><td>Turkey/98-Denizli</td><td>1998</td><td>Goat</td><td>Turkey</td><td>West Asia</td><td>FJ869356</td></tr><tr><td>Maharashtra/12</td><td>2008</td><td>Goat</td><td>India</td><td>South Asia</td><td>KF495246</td></tr><tr><td>Maharashtra/Sheep/19</td><td>2010</td><td>Sheep</td><td>India</td><td>South Asia</td><td>KF495252</td></tr><tr><td>Mukteswar/1946_P-4</td><td>1946</td><td>Goat</td><td>India</td><td>South Asia</td><td>KF495245</td></tr><tr><td>* MOG/GP/L/5/08</td><td>2008</td><td>Goat</td><td>Mongolia</td><td>East Asia</td><td>MN161845</td></tr><tr><td>* MOG/GP/T/4/08</td><td>2008</td><td>Goat</td><td>Mongolia</td><td>East Asia</td><td>MN161846</td></tr><tr><td>* MOG/GP/T/6/08</td><td>2008</td><td>Goat</td><td>Mongolia</td><td>East Asia</td><td>MN161847</td></tr><tr><td>Akola/2008_P-2</td><td>2008</td><td>Goat</td><td>India</td><td>South Asia</td><td>KF495244</td></tr><tr><td>Sambalpur/1982_P-9</td><td>1982</td><td>Goat</td><td>India</td><td>South Asia</td><td>KF495243</td></tr><tr><td>Uttarkashi/1978</td><td>1978</td><td>Goat</td><td>India</td><td>South Asia</td><td>KF495242</td></tr><tr><td>Jammu_Kashmir/Goat/27</td><td>2013</td><td>Goat</td><td>India</td><td>South Asia</td><td>KF495249</td></tr><tr><td>Jammu_Kashmir/Sheep/27</td><td>2013</td><td>Sheep</td><td>India</td><td>South Asia</td><td>KF495248</td></tr><tr><td>GS-V1</td><td>2011</td><td>Unknown</td><td>China</td><td>East Asia</td><td>JQ310667</td></tr><tr><td>India/83</td><td>1983</td><td>Goat</td><td>India</td><td>South Asia</td><td>FJ869358</td></tr><tr><td>SPPV_Oman/84</td><td>1984</td><td>Sheep</td><td>Oman</td><td>West Asia</td><td>FJ869390</td></tr><tr><td>Oman/84</td><td>1984</td><td>Goat</td><td>Oman</td><td>West Asia</td><td>FJ869359</td></tr><tr><td>Bangladesh/86</td><td>1986</td><td>Goat</td><td>Bangladesh</td><td>South Asia</td><td>FJ869355</td></tr><tr><td>FZ</td><td>2012</td><td>Goat</td><td>China</td><td>East Asia</td><td>KC951854</td></tr><tr><td>AV41</td><td>2018</td><td>Goat</td><td>China</td><td>East Asia</td><td>MH381810</td></tr><tr><td>2015/Gansu/China</td><td>2015</td><td>Unknown</td><td>China</td><td>East Asia</td><td>MG817383</td></tr><tr><td>HuB/2009/China</td><td>2009</td><td>Goat</td><td>China</td><td>East Asia</td><td>JQ310672</td></tr></tbody></table>

## 2.2. Collection of GTPVs for Alignment-Free Method

The sequences generated for twelve new isolates of this study were analyzed with publicly available GTPVs sequences. Nucleotide sequences and deduced amino acid sequences were aligned using the ClustalW algorithm implemented in MEGA7 [16]. The

<footnote>* GTPV sequenced for this study.</footnote>
complete GPCR gene sequences of 37 additional GTPVs from various regions were retrieved from GenBank. For phylogenetic reconstructions, sequences representative of LSDV and SPPV were also retrieved from GenBank and included in the analysis.

Firstly, the GPCR nucleotides sequences were aligned using the ClustalW algorithm (codon option) implemented in MEGA 7. The aligned sequences in FASTA format were converted to nexus format using Seaview. Secondly, the Bayesian phylogenetic inference was performed with BEAST v1.8.4 [17]. The BEAUti module was used to generate BEAST files from the sequence alignments in nexus format. The HKY substitution +G nucleotide substitution was used with a UPGMA starting tree. Then, the Markov Chain Monte Carlo method was run with BEAST, for 10,000,000 generations with a sample taken each 10,000 generations. The TRACER programme was used to inspect the log files and determine the optimum number of burn-in based on the Effective Sample Sizes (ESS > 200). TreeAnnotator was used to generate the Maximum Clade Credibility (MCC) after discarding the 5% burn-in. Finally, the tree was visualized with associated meta-data using the ggtree package in R version 3.5.2 [18].

## 2.3. *Distances Calculation and Visualization*

Only the sequences that clustered within GTPVs group in the phylogenetic analysis were used in the alignment-free analysis.

The multi-fasta file with the GPCR gene sequences of 49 GTPVs was loaded into the pairwise analysis interface of Heatmapper (http://www.heatmapper.ca/pairwise/, accessed on 14 May 2019). Heatmapper [15] was used to compute k-mer frequencies for each sequence and produce the distance matrix table using an alignment-free algorithm. Briefly, k-mer size of 3 was used to compute the k-mer frequencies for each sequence; then the Euclidian distance measurement method was used to calculate the pairwise distances between each pair of sequences. The distance matrix was visualized as a heatmap with "gplots" package in R version 3.5.2 and the Complete-linkage clustering method was implemented to re-order the sequences.

Additionally, three alignment-based algorithms were used to generate distance matrices (p-distance, Maximum Composite Likelihood (MCL) and Kimura 2-parameter, all implemented in MEGA software) and produce heatmaps. To examine the influence of the gaps, the positions corresponding to the 21-nucleotide deletion of East and South Asian GTPVs were manually removed and the pattern of distance matrix with alignment-free algorithm was analyzed.

# 3. Results

## 3.1. *Sequencing and Phylogenetic Reconstructions*

We have successfully sequenced the full GPCR gene for all 12 samples collected from both sheep and goats from Ethiopia, Ghana, Kenya and Mongolia and compared that to GTPVs retrieved from GenBank.

The maximum credibility tree of the GPCR sequences showed that all 12 sequences of this study belonged to GTPVs (Figure 1). In the phylogenetic tree, 37 additional sequences retrieved from GenBank clustered with GTPVs. GTPVs formed two distinct groups: one with African isolates and a second one with Asian isolates (Figure 1). The Asian GTPVs produced two sub-clusters: sub-cluster 1 contained GTPVs from Eastern and Southern Asia (ESA) mainly and sub-cluster 2 consisted of isolates of Western and Central Asia (WCA). Among the newly sequenced GTPVs, one isolate from Ghana, five isolates from Ethiopia and three from Kenya, clustered with African GTPVs, while the three GTPVs from Mongolia, clustered with ESA GTPVs (Figure 1).

Some isolates retrieved from GenBank clustered outside the groups matching their geographical origin. SPPV Oman/84 (FJ869390), GTPV Oman/84 (FJ869359) and GTPV SA2/2017 (MG232389), from Western Asia, clustered with ESA GTPVs, while GTPV Yemen/83(FJ869362) also from Western Asia, clustered with African GTPVs (Figure 1).
To further determine whether the GPCR gene sequences could help classify GTPVs based on their geographical origins, we selected all 49 GTPVs (37 from GenBank and 12 from this study) for alignment-free sequence comparison.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Maximum clade credibility (MCC) tree based on of the complete GPCR complete gene sequences of capripoxviruses. The posterior probabilities are plotted as respective nodes labels. The GTPV sequences are highlighted based their geographical origins (purple for Africa, green for West and Central Asia, blue for East and South Asia). Note GTPV Yemen clustering with African GTPV and GTPV Oman and SPPV Oman, clustering with East and South Asia GTPVs.</figcaption></figure>

## 3.2. *Distances Calculation and Visualization Using Alignment-Free Methods*

In the heatmap, based on the Euclidian distances calculated from the k-mer frequencies of each pair of sequences, we observed two main clusters: cluster 1 comprised isolates from Asia and cluster 2 comprised isolates from Africa (Figure 2). Cluster 1 with Asian GTPVs further split into two sub-clusters: sub-cluster 1.1 comprising ESA GTPVs and sub-cluster 1.2 containing WCA GTPVs (Figure 2). In general, all GTPVs clustered according to their geographical origins, except GTPV Yemen/83 (FJ869390) from Western Asia, located within cluster 2 of African GTPVs and GTPV Oman/84(FJ869359) and SPPV Oman/84 (FJ869390),
also from Western Asia, located in sub-cluster 1.1 of ESA GTPVs (Figure 2). Contrasting with the result of the phylogenetic reconstruction, GTPV SA2/2017 (MG232389) clustered within sub-cluster 1.1 of WCA GTPVs (Figure 2).

<figure><img src="image_6.png" /></figure>

**Figure 2.** The heatmap for the nucleotides k-mers frequencies variations of 49 GTPVs. The Complete-linkage clustering method was used to re-order the sequences. The vertical side colors indicate the origin of the samples (red for Africa, orange for West and Central Asia, yellow for East and South Asia). The GTPVs recovered from goats are shown in blue, those from sheep in red and those with unknown host in black. The clusters and sub-clusters (Scl) are indicated with boxes.

Heatmaps based on distances calculated from multiple sequence alignments using various substitution models (p-distance, MCL and k2P) produced similar patterns, however with lower resolution. While the two blocks of Asian and African GTPV were apparent, the separation between the two sub-groups of Asian isolates was less prominent and the GTPV SA2/2017 (MG232389) isolate clustered within the ESA GTPVs sub-block 2.1 (Supplement Figure S1). Likewise, after removing the positions containing gaps before calculating the Euclidian distances from the k-mer frequencies, the resulting heatmap, appeared to have a lower resolution, with GTPV SA2/2017 (MG232389) clustering within the ESA GTPVs sub-block 1.2 (Supplement Figure S2).

### 3.3. *Amino Acids Profiles of GTPVs' GPCR*

The comparative analysis of the amino acid sequence alignments showed 10 conserved amino acid (aa) differences between African and WCA GTPVs and 17 aa differences (7 aa deletion and 10 aa substitution) between African and ESA GTPVs out of which, nine were common to all Asian GTPVs. The WCA GTPVs differed from ESA GTPVs by seven amino acids, missing in the ESA GTPVs (Supplement Figure S3). In summary, based on the amino acid profile, we can distinguish three different clusters of GTPVs: the African GTPVs, the WCA GTPVs and the ESA GTPVs (Supplement Figure S3). Most of the 49 GTPVs of this study had a GTPV profile that matched the geographical origin except SPPV Oman/84 (FJ869390) and GTPV Oman/84 (FJ869359), which presented the ESA GTPV profile and GTPV Yemen/83 (FJ869390), with the African GTPV profile. Based on the amino acid sequence, GTPV SA2/2017 (MG232389) clustered within the WCA GTPVs sub-group, similar to the result obtained using the heatmap based on the Euclidian distances calculated from the k-mer frequencies.
GTPV SA2/2017 was the only Western Asian isolate sharing an additional amino acid similarity with ESA GTPVs: H was present at position 375 instead of R.

## 4. Discussion

We analyzed the GPCR sequences of GTPVs from various locations, including 12 newly sequenced isolates from Africa and Asia, using multiple sequence alignments and alignment-free methods.

The analysis of the maximum credibility tree of the GPCR gene sequences confirmed that the newly sequenced isolates belong to GTPVs and enabled the identification of 37 additional GTPVs included for alignment-free analysis. Based on the phylogenetic analysis, the alignment-free sequence comparisons and the information attached to the sequences, we identified two main clusters: the African GTPVs and the Asian GTPVs, which further subdivided into ESA GTPVs and WCA GTPVs.

In the phylogenetic reconstruction, four isolates, SPPV Oman/84, known to be a GTPV [2], GTPV Oman/84, GTPV Yemen/83 and GTPV SA2/2017, all from Western Asia, clustered outside the groups corresponding to their geographical origin. Contrasting with the phylogenetic analysis, the alignment-free method based on the Euclidian distances calculated from the k-mer frequencies [15], correctly classified GTPV SA2/2017 within the WCA GTPVs, consistently with the results of the comparative analysis of the amino acid profiles of the GTPVs. The failure to correctly classify GTPV SA2/2017 using conventional sequence comparison methods is likely due to the mishandling of the 21-nucleotide deletion observed in ESA GTPVs. Indeed, the heatmaps produced using distance calculated based on conventional methods such as MCL, K2P and p-distance and the one based on k-mer frequencies after removal of the positions with gaps, also failed to classify GTPV SA2/2017 in the correct cluster. Our study further highlights the power of the alignment-free methods in handling useful information such as gaps. Contrasting with global and local alignment algorithms which work base-by-base, frequency-based algorithms can better deal with complexities caused by mismatches, gaps and sequence inversions often found while comparing sequences [13].

In a previous study, using thirteen GTPVs GPCR sequences, Le Goff et al. (2009) [2], observed three groups: however, due to the limited number of sequences and the absence of sequences from Eastern Asia and Eastern Africa, they were unable to draw any definite conclusion for the classification of GTPVs based on their geographical origins.

Though the heatmap based on the alignment-free method provided a better classification of GTPVs according to their geographical origins, this method identified GTPV Yemen/83 as an African GTPV and GTPV Oman/84 and SPPV Oman/84 as ESA GTPVs.

This misclassification could also suggest that GTPV Yemen/83 is an African GTPV imported to Yemen and GTPV Oman and SPPV Oman are two ESA GTPVs imported to Oman. In the 1980s, when GTPV Yemen/83 was collected, Yemen imported sheep and goats from Eastern Africa, mainly, Somalia, Sudan, Eritrea and Sudan, before the outbreaks [19]. It was suggested that new GTPV strains, imported via Eastern Africa, could have caused these outbreaks [19]. Similarly, Oman imported sheep and goats from Africa, other countries in the Middle East and India, during the same period, with the risk of introducing new GTPVs in the country [19]. The introduction of foreign GTPV strains from Africa to Yemen and/or South Asia to Oman via trade or illegal imports is also consistent with livestock trading routes between African and Asian countries [20].

The spread of foot and mouth disease (FMD) in the Middle East, which better illustrates the association between trade and introducing transboundary diseases, was associated with the importation of animals from Asia and Africa leading to epidemics involving new FMD virus strains [20]. The presence of all three GTPV profiles in the Western Asia further supports the hypothesis of GTPVs importation in the region via livestock trade.

These geographical signatures in the GPCR were highly conserved, with no change observed for each region over 50 years. For instance, the WCA GTPVs included in this study were from outbreaks that occurred between 1961 and 2017, the ESA GTPVs, between
1946 and 2009 and the African GTPVs, between 1955 and 2013, implying a great adaptation of these profiles for each region.

The CaPVs' GPCR gene enabled the discrimination of GTPVs based on their geographical origins. While the GPCR gene encodes for a host factor protein, it remains unclear whether the GPCR profile was driven by the adaptation of GTPV to the affected breeds in each region.

It is also worth noting that some ESA and African GTPVs, including six African GTPVs that were sequenced for this study, were from sheep, showing that these isolates can affect both small ruminant species. In contrast, though only a few GTPVs were available from WCA, we noticed that all of them were exclusively from goats.

**5. Conclusions**

This study showed that, using alignment-free sequence comparison methods, the GPCR gene sequences enable the classification of GTPVs based on geographical origin.

The alignment-free approach based on k-mer frequencies enabled a higher resolution as compared to standard approaches using alignment to calculate distance matrix. The result of the alignment-free approach correlates better with the specific features found in the amino acid sequence comparison.

The approach described here can help support investigations to determine GTPV's origin of introduction whenever it emerges in a disease-free region and determine the origin of contamination with a new strain in GTPV endemic areas.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/microorganisms9040855/s1. Figure S1: The heatmap for the nucleotide's variation of 49 GTPVs based on the MLC distances. The Complete-linkage clustering method was used to re-order the sequences. The vertical side colors indicate the origin of the samples (red for Africa, orange for West and Central Asia, yellow for East and South Asia). Figure S2: The heatmap for the nucleotides k-mers frequencies variations of 49 GTPVs after manual removal of 21 nucleotides from the GPCR genes of WCA and African GTPV. The 21 nucleotides correspond to the gap observed on the multiple sequence alignment in ESA GTPV. The Complete-linkage clustering method was used to re-order the sequences. The vertical side colors indicate the origin of the samples (red for Africa, orange for West and Central Asia, yellow for East and South Asia). Figure S3: Multiple sequence alignment of the complete amino acid sequences of the GPCR of 49 GTPVs from various geographical locations.

**Author Contributions:** Conceptualization, C.E.L. and T.R.C.; methodology, T.R.C., Y.L., C.E.L., R.G. and A.L., software, Y.L.; formal analysis, T.R.C., C.E.L. and T.B.K.S.; investigation, T.R.C., J.K.L., D.D., B.B. and J.E.; resources, A.L., M.S., J.K.L., J.E., B.B. and D.D.; data curation, M.S., J.E. and D.D.; writing—original draft preparation, T.R.C., C.E.L., F.J.B. and R.G.; writing—review and editing, F.J.B., A.D., G.C., R.G., A.L., M.S., J.K.L., T.B.K.S., B.B. and J.E. visualization, Y.L. and C.E.L.; supervision, A.D., G.C., R.G. and C.E.L.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the VETLAB network initiative of the Joint FAO/IAEA Division, funded through the African Renaissance and International Cooperation Fund of South Africa and the Peaceful Uses Initiatives (PUI) funded by Japan and the United States of America.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are openly available in NCBI at https://www.ncbi.nlm.nih.gov/nuccore/, accessed on 15 April 2021.

**Acknowledgments:** The authors would like Jenna E. Achenbach (Battelle, VA, USA) for the help rendered the assessment of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.
References

1. Buller, R.M.; Arif, B.M.; Black, D.N.; Dumbell, K.R.; Esposito, J.J.; Lefkowitz, E.J.; McFadden, G.; Moss, B.; Mercer, A.A.; Moyer, R.W.; et al. Family poxviridae. In *Virus Taxonomy: Classification and Nomenclature of Viruses*. Eighth Report of the International Committee on Taxonomy of Viruses; Mayo, M.A., Maniloff, J., Desselberger, U., Ball, L.A., Eds.; Elsevier Academic Press: San Diego, CA, USA, 2005; pp. 117–133.

2. Le Goff, C.; Lamien, C.E.; Fakhfakh, E.; Chadeyras, A.; Aba-Adulugba, E.; Libeau, G.; Tuppurainen, E.; Wallace, D.B.; Adam, T.; Silber, R.; et al. Capripoxvirus G-protein-coupled chemokine receptor: A host-range gene suitable for virus animal origin discrimination. *J. Gen. Virol.* **2009**, *90*, 1967–1977. [CrossRef] [PubMed]

3. Dutta, T.K.; Roychoudhury, P.; Kawlni, L.; Lalmuanpuia, J.; Dey, A.; Muthuchelvan, D.; Mandakini, R.; Sarkar, A.; Ramakrishnan, M.A.; Subudhi, P.K. An outbreak of Goatpox virus infection in Wild Red Serow (*Capricornis rubidus*) in Mizoram, India. *Transbound. Emerg. Dis.* **2019**, *66*, 181–185. [CrossRef] [PubMed]

4. Kitching, R.P.; Taylor, W.P. Clinical and antigenic relationship between isolates of sheep and goat pox viruses. *Trop. Anim. Health Prod.* **1985**, *17*, 64–74. [CrossRef] [PubMed]

5. Lamien, C.E.; Lelenta, M.; Goger, W.; Silber, R.; Tuppurainen, E.; Matijevic, M.; Luckins, A.G.; Diallo, A. Real time PCR method for simultaneous detection, quantitation and differentiation of capripoxviruses. *J. Virol. Methods* **2011**, *171*, 134–140. [CrossRef] [PubMed]

6. Lamien, C.E.; Le Goff, C.; Silber, R.; Wallace, D.B.; Gulyaz, V.; Tuppurainen, E.; Madani, H.; Caufour, P.; Adam, T.; El Harrak, M.; et al. Use of the Capripoxvirus homologue of Vaccinia virus 30 kDa RNA polymerase subunit (RPO30) gene as a novel diagnostic and genotyping target: Development of a classical PCR method to differentiate Goat poxvirus from Sheep poxvirus. *Vet. Microbiol.* **2011**, *149*, 30–39. [CrossRef] [PubMed]

7. Tuppurainen, E.S.; Venter, E.H.; Shisler, J.L.; Gari, G.; Mekonnen, G.A.; Juleff, N.; Lyons, N.A.; De Clercq, K.; Upton, C.; Bowden, T.R.; et al. Review: Capripoxvirus Diseases: Current Status and Opportunities for Control. *Transbound. Emerg. Dis.* **2017**, *64*, 729–745. [CrossRef] [PubMed]

8. Gelaye, E.; Belay, A.; Ayelet, G.; Jenberie, S.; Yami, M.; Loitsch, A.; Tuppurainen, E.; Grabherr, R.; Diallo, A.; Lamien, C.E. Capripox disease in Ethiopia: Genetic differences between field isolates and vaccine strain, and implications for vaccination failure. *Antivir. Res.* **2015**, *119*, 28–35. [CrossRef] [PubMed]

9. Santhamani, R.; Venkatesan, G.; Minhas, S.K.; Shivachandra, S.B.; Muthuchelvan, D.; Pandey, A.B.; Ramakrishnan, M.A. Detection and characterization of atypical capripoxviruses among small ruminants in India. *Virus Genes* **2015**, *51*, 33–38. [CrossRef] [PubMed]

10. Santhamani, R.; Yogisharadhya, R.; Venkatesan, G.; Shivachandra, S.B.; Pandey, A.B.; Ramakrishnan, M.A. Molecular characterization of Indian sheeppox and goatpox viruses based on RPO30 and GPCR genes. *Virus Genes* **2014**, *49*, 286–291. [CrossRef] [PubMed]

11. Zhou, T.; Jia, H.; Chen, G.; He, X.; Fang, Y.; Wang, X.; Guan, Q.; Zeng, S.; Cui, Q.; Jing, Z. Phylogenetic analysis of Chinese sheeppox and goatpox virus isolates. *Virol. J.* **2012**, *9*, 1–8. [CrossRef] [PubMed]

12. Zielezinski, A.; Vinga, S.; Almeida, J.; Karlowski, W.M. Alignment-free sequence comparison: Benefits, applications, and tools. *Genome Biol.* **2017**, *18*, 1–7. [CrossRef] [PubMed]

13. Vinga, S.; Almeida, J. Alignment-free sequence comparison—A review. *Bioinformatics* **2003**, *19*, 513–523. [CrossRef] [PubMed]

14. Bonham-Carter, O.; Steele, J.; Bastola, D. Alignment-free genetic sequence comparisons: A review of recent approaches by word analysis. *Brief. Bioinform.* **2013**, *15*, 890–905. [CrossRef] [PubMed]

15. Babicki, S.; Arndt, D.; Marcu, A.; Liang, Y.; Grant, J.R.; Maciejewski, A.; Wishart, D.S. Heatmapper: Web-enabled heat mapping for all. *Nucleic Acids Res.* **2016**, *44*, W147–W153. [CrossRef] [PubMed]

16. Kumar, S.; Stecher, G.; Tamura, K. MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. *Mol. Biol. Evol.* **2016**, *33*, 1870–1874. [CrossRef] [PubMed]

17. Drummond, A.J.; Suchard, M.A.; Xie, D.; Rambaut, A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. *Mol. Biol. Evol.* **2012**, *29*, 1969–1973. [CrossRef] [PubMed]

18. Yu, G.; Smith, D.K.; Zhu, H.; Guan, Y.; Lam, T.T. ggtree: An r package for visualization and annotation of phylogenetic trees with their covariates and other associated data. *Methods Ecol. Evol.* **2017**, *8*, 28–36. [CrossRef]

19. Kitching, R.P.; McGrane, J.J.; Taylor, W.P. Capripox in the Yemen Arab Republic and the Sultanate of Oman. *Trop. Anim. Health Prod.* **1986**, *18*, 115–122. [CrossRef] [PubMed]

20. Di Nardo, A.; Knowles, N.J.; Paton, D.J. Combining livestock trade patterns with phylogenetics to help understand the spread of foot and mouth disease in sub-Saharan Africa, the Middle East and Southeast Asia. *Rev. Sci. Tech.* **2011**, *30*, 63. [CrossRef] [PubMed]