**Microbiome Shift, Diversity, and Overabundance of Opportunistic Pathogens in Bovine Digital Dermatitis Revealed by 16S rRNA Amplicon Sequencing**

**Hector M. Espiritu <sup>1</sup> , Lovelia L. Mamuad <sup>1</sup> , Seon-ho Kim <sup>1</sup> , Su-jeong Jin <sup>1</sup> , Sang-suk Lee <sup>1</sup> , Seok-won Kwon <sup>2</sup> and Yong-il Cho 1,\***


Received: 27 August 2020; Accepted: 30 September 2020; Published: 3 October 2020

**Simple Summary:** Bovine digital dermatitis (BDD) is a foot infection known as the primary cause of lameness in cattle due to painful lesions, posing serious impacts on the productivity and welfare of affected animals. Members of the bacterial group *Treponema* have long been considered as the main causative agents because previous investigations by bacterial isolation, tissue analyses, and high molecular sequencing have persistently identified this group in BDD. However, other studies indicated that the presence of several bacteria on the lesion due to the slurry environment the cattle foot are exposed to, suggests an interdependent polybacterial nature which could also play a role in disease development and progression. Therefore, we analyzed the diversity and relationship of the diverse microbiome in BDD lesions compared to normal skin from cattle foot by using next-generation high throughput sequencing. Based on the results obtained, we concluded that the shift in microbial composition which leads to richer diversity in BDD, and the overabundance of opportunistic bacterial pathogens could be associated with BDD pathogenesis.

**Abstract:** This study analyzed the diversity and phylogenetic relationship of the microbiome of bovine digital dermatitis (BDD) lesions and normal skin from cattle foot by using 16S rRNA amplicon sequencing. Three BDD samples and a normal skin sample were pre-assessed for analysis. The Illumina Miseq platform was used for sequencing and sequences were assembled and were categorized to operational taxonomic units (OTUs) based on similarity, then the core microbiome was visualized. The phylogeny was inferred using MEGA7 (Molecular evolutionary genetics analysis version 7.0). A total of 129 and 185 OTUs were uniquely observed in normal and in BDD samples, respectively. Of the 47 shared OTUs, 15 species presented increased abundance in BDD. In BDD and normal samples, Spirochetes and Proteobacteria showed the most abundant phyla, respectively, suggesting the close association of observed species in each sample group. The phylogeny revealed the evolutionary relationship of OTUs and the Euclidean distance suggested a high sequence divergence between OTUs. We concluded that a shift in the microbiome leads to richer diversity in BDD lesions, and the overabundance of opportunistic pathogens and its synergistic relationship with commensal bacteria could serve as factors in disease development. The influence of these factors should be thoroughly investigated in future studies to provide deeper insights on the pathogenesis of BDD.

**Keywords:** bovine digital dermatitis; cattle lameness; microbiome; *Treponema* spp

#### **1. Introduction**

Bovine digital dermatitis (BDD) is known as the most important foot infection causing lameness in cattle [1]. This severe lameness, which is the primary clinical manifestation of BDD caused by painful hyperkeratotic lesions [2], poses serious concerns on the welfare of the affected animals [3]. Serious economic impacts have also been implicated with BDD due to significant milk production losses, poor reproductive performance [4,5], and most extremely, premature culling of the affected animals [6].

The first reported case of BDD was in Italy in 1974 by Chelli and Mortellaro [7], and since then, it has been globally reported reaching an endemic state in many countries [8]. Etiological investigations have identified a variety of bacteria in BDD lesions [9], but advancements in sequencing technology have provided essential information on the identity of associated causal agents [8]. In previous investigations, Spirochetes are considered to be the major pathogen in BDD [1,10]. However, studies indicated that the presence of several bacteria on the lesion due to the slurry environment the cattle foot are exposed to [2], suggests a synergistic polybacterial nature [11] which plays a role in disease development and progression [12].

Many reports have attempted to elucidate the association of the microbiome of BDD towards disease development [11,13–15]. Surprisingly, for more than 20 years since it emerged as a major problem in the cattle industry in certain countries, studies on BDD have just recently started in Korea. In a recent report, *Treponema* spp. was established as the dominant pathogen in BDD lesions in Korea [16], but an emphasis on the plethora of bacteria that was observed is still lacking. Although the majority of studies considered that the main pathogen of BDD is *Treponema* spp., its polybacterial nature must also be considered in order to understand its pathogenesis [15]. Therefore, this study analyzed the diversity and phylogenetic relationship of the microbiome of BDD lesions and normal skin from the interdigital space of the cattle using 16S rRNA amplicon sequencing.

#### **2. Materials and Methods**

#### *2.1. Sample Collection*

This study was performed following general ethical principles and with the consent of farm owners. The collection of samples was done in Holstein–Friesian cattle in housed dairy farm, and was conducted during hoof cleaning, trimming and treatment performed by a professional veterinarian with years of expertise in the field. Inspection was done by anatomic pathological observation of grossly visible active BDD lesions on the proximal border of interdigital space characterized by the presence of ulceration, with hyperkeratosis and proliferative growth with hair-like projections as described by Zinicola et al. (2015) [17]. For normal skin sample, healthy animals with no sign of lameness and no history of BDD were inspected. After thorough cleaning of the foot surface, lidocaine (2%) was subcutaneously injected around the lesion and a 5-mm punch biopsy was taken from the center of the lesion for BDD, while for normal skin, biopsy sample was taken where BDD most often take place. Samples were washed thoroughly with buffered phosphate saline (pH 7.4) and delivered to the laboratory with ice. A total of 66 pre-assessed active BDD samples were subjected to detection of *Treponema* spp. by genus-specific PCR as described in the results of our previous study [16]. Since 100% of the active BDD lesions are positive in PCR, three samples were randomly selected for metagenomics sequencing. A normal sample was provided for BDD negative control. Samples were categorized as Group A for normal skin sample, and Group B for BDD-infected samples.

#### *2.2. Metagenomics Sequencing and Diversity Analysis*

The four randomly pre-assessed samples were submitted to Macrogen (Korea) for high-throughput sequencing. DNA samples extracted from lesion and normal sample biopsies were subjected to quality control by Picogreen method before library construction. By targeting the V3-V4 region of the 16S rRNA gene, libraries were constructed and were purified. Sequencing was carried out using

Illumina Miseq platform (Illumina, San Diego, CA, USA). The base call binary data produced by Real-Time Analysis (RTA) were converted to FASTQ files by bcl2fastq package (Illumina, San Diego, CA, USA) and were filtered using Scythe (v0.994) (https://github.com/vsbuffalo/scythe) and Sickle (https://github.com/najoshi/sickle) programs to remove adapter sequences. The obtained 16S rRNA sequences were binned into Operational Taxonomic Units (OTUs) based on 97% identity using Quantitative Insights Into Microbial Ecology (QIIME) [18]. The microbiome was visualized using Metagenomics Core Microbiome Exploration Tool (MetaCoMET) [19] using a Biological Observation Matrix format (BIOM) [20] generated using Mothur [21].

#### *2.3. Divergence and Phylogenetic Analysis*

The phylogeny of the microbiome was inferred by aligning the obtained 16S rRNA sequences of the operational taxonomic units identified using ClustalW [22], and the Newick tree data were obtained from MEGA7 [23] by maximum likelihood method following the general time-reversible model as the fit model for calculating the rate of nucleotide base substitution. The final dendrogram was constructed using Iroki [24]. The divergence distance between OTUs was calculated by Euclidian distance method based on the rate of nucleotide base substitution between OTUs, and presented as a heatmap matrix generated by Heatmapper [25] for the aligned sequences of OTUs using the Euclidean distance method.

#### **3. Results**

In total, 66 samples pre-assessed from our previous research [16] were used in this study. From these, one normal skin sample (Group A) and three BDD lesions (Group B) were subjected for analysis.

The Chao1 diversity index shows that Group B has richer species diversity compared to Group A (Figure 1A). The size of the microbiome presented in Figure 1B shows a higher observed OTUs in BDD lesions compared to the normal skin sample. Out of 267 OTUs observed, 138 were uniquely observed in Group B, 82 in Group A, while 47 were overlapping in both groups. Fifteen of these shared OTUs increased their abundance from normal to BDD and includes T. pedis (20.93%), a group of unclassified species (12.4%), Treponema denticola (9.8%), T. medium (6.48%), Porphyromonas levii (1.56%), P. somerae (1.22%), and Acholeplasma vituli (0.89%). OTUs absent in Group A with increased abundance ratio in Group B were as follows in decreasing order: Carboxylicivirga mesophila (5.89%), T. lecithinolyticum (5.35%), A. morum (5.03%), Spirochaeta africana (3.74%), Mycoplasma feliminutum (3.65%), A. modicum (2.77%), Pelobacter propionicus (2.04%), H. sueciensis (1.64%), P. uenonis (1.30%), Falcatimonas natans (1.25%), Devosia confluentis (1.17), Christensenella minuta (1.08%), and M. fermentans (1.30%). Additional information can be accessed in Supplementary Table S1. The shift in abundance ratio of each species from normal to BDD (Figure 1C) shows the drop in abundance of dominant bacteria, Psychrobacter fulvigenes and Pseudomonas caeni from the normal sample, and the increased abundance of Treponema spp., and other bacteria in BDD. In addition, the Euclidian distance heat-map based on OTU abundance between samples and representative phyla (Figure 1D) illustrated the high association of Spirochetes with BDD. Moraxellaceae and Pseudomonadaceae, both under phylum Proteobacteria presented a close relationship with normal tissue. Moreover, Bacteroidetes and Firmicutes are considerably associated in both groups, while Tenericutes is associated with BDD only.

*Animals* **2020**, *10*, x FOR PEER REVIEW 4 of 9

**Figure 1.** (**A**) Box plot of Chao1 diversity index for bovine digital dermatitis (BDD) (Group A) and normal (Group B), (**B**) the microbiome presenting the number of operational taxonomic units (OTUs) in Group A and Group B showing the unique and shared OTUs between groups, and (**C**) relative abundance of operational taxonomic units from Group A and Group B (legend: 22 most abundant). (**D**) Euclidian distance heatmap presenting the association between phylum and sample groups. Red denotes closer association while blue denotes a lesser association with normal sample and BDD. **Figure 1.** (**A**) Box plot of Chao1 diversity index for bovine digital dermatitis (BDD) (Group A) and normal (Group B), (**B**) the microbiome presenting the number of operational taxonomic units (OTUs) in Group A and Group B showing the unique and shared OTUs between groups, and (**C**) relative abundance of operational taxonomic units from Group A and Group B (legend: 22 most abundant). (**D**) Euclidian distance heatmap presenting the association between phylum and sample groups. Red denotes closer association while blue denotes a lesser association with normal sample and BDD.

The maximum likelihood tree constructed using the general time-reversible model illustrated the phylogenetic relationship of all the OTUs (Figure. 2A). The tree showed the phylum classification of the OTUs, and the abundance ratio and OTU count (Figure 2B). Firmicutes was the most diverse phylum representing 41.85% and 37.21% of the observed OTUs in BDD and normal skin, respectively. Spirochetes (which is the most abundant) only represented 3.75% of the over-all observed OTUs. In addition, the Euclidean distance heat map matrix showed the estimates of the evolutionary divergence between OTUs in lesions and normal tissue (Figure 2C), and the frequency distribution of the computed distance of the pairwise comparisons was graphed in Figure 2D. Overall, 63.89% of the pairwise comparisons were above the median distance. The maximum likelihood tree constructed using the general time-reversible model illustrated the phylogenetic relationship of all the OTUs (Figure 2A). The tree showed the phylum classification of the OTUs, and the abundance ratio and OTU count (Figure 2B). Firmicutes was the most diverse phylum representing 41.85% and 37.21% of the observed OTUs in BDD and normal skin, respectively. Spirochetes (which is the most abundant) only represented 3.75% of the over-all observed OTUs. In addition, the Euclidean distance heat map matrix showed the estimates of the evolutionary divergence between OTUs in lesions and normal tissue (Figure 2C), and the frequency distribution of the computed distance of the pairwise comparisons was graphed in Figure 2D. Overall, 63.89% of the pairwise comparisons were above the median distance.

*Animals* **2020**, *10*, x FOR PEER REVIEW 5 of 9

**Figure 2.** General time-reversible tree representing all observed OTUs in BDD-infected and normal skin samples color-coded based on phylum classification (**A**). The scale value is 2.57. The blue bar graph denotes the abundance ratio of each species on BDD (**B1**) and normal sample (**B2**), while the green bar graph shows the actual count per species (OTU) for BDD (**B3**) and normal (**B4**) sample. Euclidean distance heat map matrix is based on the nucleotide base substitution rate presenting the divergence between all OTUs in BDD and normal sample (**C**). The graph (**D**) shows the frequency distribution of pairwise comparison between OTUs in each calculated distance from low (blue) to **Figure 2.** General time-reversible tree representing all observed OTUs in BDD-infected and normal skin samples color-coded based on phylum classification (**A**). The scale value is 2.57. The blue bar graph denotes the abundance ratio of each species on BDD (**B1**) and normal sample (**B2**), while the green bar graph shows the actual count per species (OTU) for BDD (**B3**) and normal (**B4**) sample. Euclidean distance heat map matrix is based on the nucleotide base substitution rate presenting the divergence between all OTUs in BDD and normal sample (**C**). The graph (**D**) shows the frequency distribution of pairwise comparison between OTUs in each calculated distance from low (blue) to high (red).

#### high (red). **4. Discussion**

pathogenesis.

**4. Discussion**  This study analyzed the diversity and phylogenetic relationship of the microbiome of BDD lesions compared to normal skin from the interdigital space of the cattle. From our previous study [16], we elucidated the relative abundance of Spirochetes, specifically *Treponema* spp. in BDD lesions and concluded that *Treponema* spp., was the dominant pathogen involved in BDD in Korea. Although the major abundance of *Treponema* spp. in lesions based on several investigations suggest that BDD is polytreponemal, the presence of other bacteria implies that they may also play certain roles in its This study analyzed the diversity and phylogenetic relationship of the microbiome of BDD lesions compared to normal skin from the interdigital space of the cattle. From our previous study [16], we elucidated the relative abundance of Spirochetes, specifically *Treponema* spp. in BDD lesions and concluded that *Treponema* spp., was the dominant pathogen involved in BDD in Korea. Although the major abundance of *Treponema* spp. in lesions based on several investigations suggest that BDD is polytreponemal, the presence of other bacteria implies that they may also play certain roles in its pathogenesis.

The data analyzed exhibited a change in the microbiome in BDD lesions from the normal skin

The data analyzed exhibited a change in the microbiome in BDD lesions from the normal skin sample. This altered microbiome has been previously observed by Krull et al., (2014) [13] and Zinicola et al., (2015) [11] by investigating the microbiome of each stage and layer of BDD lesions, respectively. They observed the replacement of species from the normal sample by other species along with the disease progression and the alteration of the microbiome in each layer of lesions. In the current data, the 82 OTUs from the normal sample were replaced in BDD-infected sample by 138 OTUs. This alteration also led to a richer diversity in BDD as opposed to the findings of Krull et al., (2014) where disease progression leads to a drop in diversity [13]. Although samples subjected in our study is categorized only as active BDD, it is recommended to use several normal samples and BDD samples categorized in varying lesion stages to gain more conclusive data. Another interesting finding of this study is the shared OTUs between normal and BDD found in the overlap comprised of three *Treponema* spp. (*T. pedis*, *T. medium*, *T. denticola*), along with two *Phorphyromonas* spp. as the five most abundant bacteria in BDD. The increased abundance of these species from normal to BDD, along with others could imply that these bacteria are commensals which later progresses as opportunistic pathological agents which initiates disease development when triggered under favorable conditions.

The abundance data previously presented by Mamuad et al., (2020) [16], to some extent, agree with the reports of Moreira et al., (2018) [14], Zinicola et al., (2015) [11], Beninger et al., (2016) [26], Krull et al., (2014) [13], and Nielsen et al., (2016) [12], that BDD is polybacterial, with *Treponema* spp., as the most abundant genus. However, the dominance and diversity of species under this genus still varies between reports from different geographical locations. Moreira et al., (2018) reported that *Treponema pedis* is the most abundant species in BDD in Brazil, which is in accordance with Mamuad et al., (2020). However, in one study in the USA, *T. denticola* and *T. phagedenis* were the most abundant species in active and inactive lesions, respectively [11]. Additionally, the results of investigation of Beninger et al., (2015) [26] and Krull et al., (2014) [13] were in agreement with each other that the most abundant species in BDD is *T. phagedenis*, which was absent in all samples in our data. This suggests that BDD microbiome has geographical and/or sample-to-sample variation.

Given that certain members of genus *Treponema* are recognized as the abundant bacteria, and considered as major contributors in the development and progression of the disease, other bacterial genera have also been reported in BDD from other countries as also in this study such as *Porphyromonas* [14], *Campylobacter*, *Acholeplasma* [11], *Peptoniphilus* and *Romboutsia* [27]. In this study, a major abundance of *Carboxylicivirga mesophila* and *Treponema lecithinolyticum* was observed as the fourth and fifth most abundant identified bacteria in BDD, respectively, which were not observed from previous studies. Further investigation is required to verify if these bacteria have actual involvement with the disease, since *C. mesophila* was first isolated from tidal flat sediment in Korea [28], while *T. lecithinolyticum* was commonly found in human oral microbiome [29].

The distance heat map between phyla and samples based on OTU count and taxonomy showed the close relationship of the observed species under phylum Spirochetes and Proteobacteria in BDD and normal skin, respectively. Phylum Bacteroidetes and Firmicutes can be associated in both groups, while Tenericutes is associated only with BDD which agrees with Nielsen et al., (2016) [12]. In a study by Bay et al., (2018) Firmicutes was the most abundant phylum in other polybacterial foot infections in bovine models such as interdigital hyperplasia, interdigital phlegmon, sole ulcer, toe necrosis, and white line disease, with Spirochetes being the fifth most abundant phylum overall, after Bacteroidetes, Actinobacteria, and Proteobacteria [27]. This suggests that there are variations on the major pathogens between these diseases.

The GTR tree shows the phylogenetic relationship of all OTUs present in both BDD-infected and normal skin samples, classified under 10 phyla, with Firmicutes representing the highest diversity. This large number of Firmicutes was supported by the findings of Yano et al., (2010) [30] and Santos et al., (2012) [31] in both normal and BDD-infected samples, regardless of its relative abundance. In this study, Spirochetes with the highest abundance in BDD have lower diversity both before and during infection. Compared with Firmicutes with highest diversity in normal and BDD, there was no increase in its

abundance before and during infection. This suggests that diversity richness of a certain phylum may be irrelative with the abundance or its pathogenic involvement in BDD. The microbiome of BDD in this study was verified to be highly diverse, thus, we hypothesize that synergism between overabundant opportunistic pathogens and the diversity of commensals makes the disease more complicated. The pairwise distances of OTUs based on the computed rate of nucleotide substitution show higher frequencies for above-median distances between species found in all samples. This supports the idea that bacterial diversity is evolutionarily diverged.

#### **5. Conclusions**

We concluded that a shift in the microbiome leads to richer diversity in BDD lesions, and the overabundance of opportunistic pathogens and its possible synergistic relationship between less abundant commensal bacteria could serve as factors in disease development and progression. Spirochetes is the most abundant phylum associated with BDD in other previous studies, and in this study, we deliberated that the abundance of species on each of the observed phylum varied between reports, suggesting either geographical or sample-to-sample variation. The influence of the overabundance of opportunistic species and the synergistic interaction of the plethora of commensal bacteria should be thoroughly investigated in future studies by including additional samples categorized to varying degree of the severity of infection, not just in BDD, but also in other lameness-related foot diseases, to provide deeper insights on the pathogenesis and microbiome relationship of these debilitating diseases.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-2615/10/10/1798/s1, Table S1: Abundance ratio of operational taxonomic units in normal and BDD samples.

**Author Contributions:** Conceptualization, Y.-i.C.; methodology, H.M.E., L.L.M., S.-j.J., S.-h.K., and S.-w.K. software, H.M.E. and L.L.M.; validation, Y.-i.C.; formal analysis, H.M.E. and L.L.M.; investigation, H.M.E. and L.L.M.; resources, Y.-i.C. and S.-w.K.; data curation, H.M.E. and L.L.M.; writing—original draft preparation, H.M.E. and L.L.M.; writing—review and editing, Y.-i.C., and S.-s.L.; visualization, H.M.E. and L.L.M.; supervision, Y.-i.C., and S.-s.L.; project administration, S.-j.J.; funding acquisition, Y.-i.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Research Foundation of Korea (NRF-2018R1C1B6004589).

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Molecular Detection and Phylogeny of** *Anaplasma phagocytophilum* **and Related Variants in Small Ruminants from Turkey**

**Münir Akta¸s \* , Sezayi Özübek and Mehmet Can Uluçe¸sme**

Department of Parasitology, Faculty of Veterinary Medicine, University of Firat, 23119 Elazig, Turkey; sozubek@firat.edu.tr (S.Ö.); mculucesme@firat.edu.tr (M.C.U.)

**\*** Correspondence: maktas@firat.edu.tr; Tel.: +90-(424)-237-0000

**Simple Summary:** We explored the existence of *Anaplasma phagocytophilum* and related variant in samples of goats and sheep obtained from Antalya and Mersin provinces, representative of Mediterranean region of Turkey. Based on *16S rRNA* and *groEL* genes of *A. phagocytophilum* and related variants, we examined blood samples by polymerase chain reaction (PCR) followed by sequencing. The results showed that the prevalence of *A. phagocytophilum* and *A. phagocytophilum*-like 1 infection was 1.4% and 26.5%, respectively. Sequencing confirmed molecular data and showed the presence of *A. phagocytophilum* and *A. phagocytophilum*-like-1 variant in the sampled animals.

**Abstract:** *Anaplasma phagocytophilum* causes tick-borne fever in small ruminants. Recently, novel *Anaplasma* variants related to *A. phagocytophilum* have been reported in ruminants from Tunisia, Italy, South Korea, Japan, and China. Based on *16S rRNA* and *groEL* genes and sequencing, we screened the frequency of *A. phagocytophilum* and related variants in 433 apparently healthy small ruminants in Turkey. *Anaplasma* spp. overall infection rates were 27.9% (121/433 analyzed samples). The frequency of *A. phagocytophilum* and *A. phagocytophilum*-like 1 infections was 1.4% and 26.5%, respectively. No *A. phagocytophilum*-like 2 was detected in the tested animals. The prevalence of *Anaplasma* spp. was comparable in species, and no significant difference was detected between sheep and goats, whereas the prevalence significantly increased with tick infestation. Sequencing confirmed PCR-RFLP data and showed the presence of *A. phagocytophilum* and *A. phagocytophilum*-like-1 variant in the sampled animals. Phylogeny-based on *16S rRNA* gene revealed the *A. phagocytophilum*-like 1 in a separate clade together with the previous isolates detected in small ruminants and ticks. In this work, *A. phagocytophilum*-like 1 has been detected for the first time in sheep and goats from Turkey. This finding revealed that the variant should be considered in the diagnosis of caprine and ovine anaplasmosis.

**Keywords:** tick-borne fever; *Anaplasma phagocytophilum*-like 1; PCR-RFLP; small ruminant

### **1. Introduction**

*Anaplasma phagocytophilum* is the agent of tick-borne fever (TBF) or pasture fever, a disease affecting some species of domestic ruminants (cattle, sheep, goats). The bacterium is a pathogenic species for livestock such as ruminants as well as humans in temperate and tropical countries [1–4]. *Anaplasma phagocytophilum* is transmitted by *Ixodes* spp. and infects host neutrophils and monocytes, where reproduction occurs [1,5]. *Anaplasma phagocytophilum* infection is known as pasture fever and characterized by fever, anorexia, lateral recumbency, dullness, and loss of milk yield in affected hosts [2,4,6].

Increased attention to *A. phagocytophilum* reveals new information about the genetic diversity of the pathogen. Recently two *Anaplasma* variants related to *A. phagocytophilum* have been documented in cattle, sheep, goats, and ticks [7–9]. In Japan, *A. phagocytophilum*-like 1 has been detected in deer and *Hemaphysalis longicornis* [10], cattle [11], *Ixodes* spp. [12], and *Haemaphysalis megaspinosa* [13]. *A. phagocytophilum*-like 2 has been identified in *Hyalomma*

**Citation:** Akta¸s, M.; Özübek, S.; Uluçe¸sme, M.C. Molecular Detection and Phylogeny of *Anaplasma phagocytophilum* and Related Variants in Small Ruminants from Turkey. *Animals* **2021**, *11*, 814. https:// doi.org/10.3390/ani11030814

Academic Editor: Robert W. Li

Received: 25 February 2021 Accepted: 11 March 2021 Published: 14 March 2021

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

**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/).

*asiaticum* [14], sheep and goats from China [15]. Recently those *Anaplasma* variants have been documented in ruminants from Tunisia [7,8], South Korea [16], and Italy [17].

Various *Anaplasma* species including *A. phagocytophilum* have been documented in ruminants and ticks in Turkey [5,6,18–21]. However, until now no data on *A. phagocytophilum* variants is available in Turkey. In the current study, 16S rRNA, groEL (heat shock protein) PCR and sequencing were performed to identity *A. phagocytophilum* and *A. phagocytophilum*like variants in small ruminants from sampling sites in Antalya and Mersin provinces, where the representative Mediterranean area of Turkey.

#### **2. Materials and Methods**

#### *2.1. Study Region and Sample Collection*

This survey was conducted in small ruminants farmed in three districts (Alanya, Akseki, Manavgat) from Antalya (latitude 36◦ 530 N, longitude 30◦ 420 E) and two districts (Anamur, Bozyazı) from Mersin (latitude 36◦ 470 N, longitude 34◦ 370 E) provinces of Turkey (Figure 1). This area has a Mediterranean climate, with hot humid summers and warm rainy winters. The goats and sheep are kept in closed areas in villages near to the coast during the winter months, and they are taken to the plateaus in the Taurus Mountains in the early spring and grazed in the pastures here until autumn. *Animals* **2021**, *11*, x 3 of 13

**Figure 1.** Map of Turkish provinces, indicating the localities studied in the study. (**A**) Geographical position of the provinces of Antalya and Mersin in Turkey. (**B**) Position of localities sampled in the provinces of Antalya and Mersin. The sample size was calculated using the online tool Sample Size Calculator **Figure 1.** Map of Turkish provinces, indicating the localities studied in the study. (**A**) Geographical position of the provinces of Antalya and Mersin in Turkey. (**B**) Position of localities sampled in the provinces of Antalya and Mersin.

Elazig Veterinary Control Institute (number: 2018/02).

*2.2. DNA Extraction and Amplification of 16S rRNA Gene*

(www.calculator.net/sample-sizecalculator.html), for a confidence level (CL) of 95%, an error margin of 5%. According to this, during April–July 2019, a total of 433 apparently

ples were drawn from the punctured jugular vein into anticoagulated (K3-EDTA) vacutainer tubes and stored at −20 °C freezer until DNA extraction. The goats and sheep were also checked for tick infestations, and a total of 1475 ticks were removed. The collected ticks were preserved in 70% ethanol in Eppendorf tubes. They were identified using taxonomic keys [22]. The animals were grouped into categories according to species (goat and sheep) and the presence of ticks (yes/no). This study secured the approval of the

DNA was isolated from 200 µL volumes of whole blood using a DNAeasy Blood Minikit according to the vendor's recommendations. Genomic DNA from blood of clinically infected cattle with *A. phagocytophilum* [6] was used for positive control. *Anaplasma* 

The sample size was calculated using the online tool Sample Size Calculator (www. calculator.net/sample-sizecalculator.html, accessed on 1 February 2019), for a confidence level (CL) of 95%, an error margin of 5%. According to this, during April–July 2019, a total of 433 apparently healthy small ruminants (296 goats, 137 sheep) were included in the survey. Blood samples were drawn from the punctured jugular vein into anticoagulated (K3-EDTA) vacutainer tubes and stored at −20 ◦C freezer until DNA extraction. The goats and sheep were also checked for tick infestations, and a total of 1475 ticks were removed. The collected ticks were preserved in 70% ethanol in Eppendorf tubes. They were identified using taxonomic keys [22]. The animals were grouped into categories according to species (goat and sheep) and the presence of ticks (yes/no). This study secured the approval of the Elazig Veterinary Control Institute (number: 2018/02).

#### *2.2. DNA Extraction and Amplification of 16S rRNA Gene*

DNA was isolated from 200 µL volumes of whole blood using a DNAeasy Blood Minikit according to the vendor's recommendations. Genomic DNA from blood of clinically infected cattle with *A. phagocytophilum* [6] was used for positive control. *Anaplasma phagocytophilum*-like variants DNAs, received from Alberto Alberti (University of Sassari, Sassari, Italy) were used as positive controls.

To detect *A. phagocytophilum* and *A. phagocytophilum*-like variants, a nested 16S rRNA PCR was carried out described by Kawahara et al. [10]. The PCR reaction conditions were made according to the previously described studies [10,21]. The nested amplicons were examined by 1.5% agarose gel electrophoresis and visualized using the gel Documentation System (Vilber Lourmat, Marne La Vallee Ceedex, France).

#### *2.3. Restriction Fragment Length Polymorphism (RFLP)*

*Xcm*I and *Bsa*I restriction enzymes allow the specific discrimination amongst *A. phagocytophilum* and related variants [8,17]. For differentiation of *A. phagocytophilum* and related variants, the nested amplicons obtained in this study were digested with the *Xcm*I and *Bsa*I restriction enzymes as previously described [8,17].

#### *2.4. GroEL PCR*

To confirm the results of the RFLP assay, the positive samples were screened by a groEL nested PCR for the amplification of *A. phagocytophilum* [23]. The semi-nested PCR reported by Ybañez et al. [24] with the primers EEGro1F/AnaGroe712R and AnaGroe240F was utilized for amplifying of *A. phagocytophilum*-like 1 *groEL* gene. Oligonucleotide primers used in this study were presented in Table 1.



#### *2.5. Sequencing and Phylogenetic Analyses*

*Anaplasma phagocytophilum* (*n* = 6) and *A. phagocytophilum*-like 1 (*n* = 10) positive PCR amplicons were purified using the QIAquick PCR Purification Kit (Qiagen, Hilden, Germany). The purified amplicons were sent to BM Labosis (Ankara, Turkey) for Sanger sequencing to determine DNA sequences of the *16S rRNA* gene. Multiple alignments were performed with the CLUSTAL Omega ver. 1.2.1 (https://www.ebi.ac.uk, accessed on 1 February 2019). The representative sequences have been submitted to the GenBank (MT881655 and MT881656 for *16S rRNA* gene of *A. phagocytophilum*-like 1 and *A. phagocytophilum*, respectively). The sequence alignment was performed using MUSCLE in Geneious prime [25].

Phylogenetic analyses of the 16S rRNA sequences obtained in this work and the other sequences submitted to GenBank were carried out. The maximum likelihood analysis (ML) carried out in Mega X [26] was utilized to determine the phylogenetic relationship of the *Anaplasma* spp. To sequence evolution, best-fit model was assessed as TN93+G+I by using the jModel test v.0.1.1 [27]. Reliability of internal branches of the tree was evaluated by the bootstrapping method with 1000 iterations [28].

#### *2.6. Statistical Analysis*

Association of the presence of *Anaplasma* spp. with host species and presence of tick was performed with Epi Info 6.01 (CDC, Atlanta), using the χ2 test and Fisher's exact test.

#### **3. Results**

#### *3.1. Tick Infestation*

Of the 433 small ruminants examined, 190 (43.9%) were infested with at least one tick species. A total of 1475 adult ticks (449 females, 1026 males) were collected from goats (1409/1475, 95.5%) and sheep (66/1475, 4.5%). Six tick species were identified among all collected ticks. *Rhipicephalus bursa* (1269/1475, 86%) was the dominant tick species, followed by *R. turanicus* (98/1475, 6.6%), *Dermacentor marginatus* (94/1475, 6.4%), *Hyalomma marginatum* (8/1475, 0.5%), *R. sanguineus* s.l. (5/1475, 0.3%), and *Ixodes ricinus* (0.06%, only one specimen). The goats were infested with all the identified tick species, whereas sheep were infested with *R. bursa* and *R. turanicus.*

#### *3.2. Prevalence and Distribution of Anaplasma spp.*

The prevalence of *A. phagocytophilum* and related variants in sampled goats and sheep is presented in Table 2. Overall, 121/433 (27.9%) samples collected in studied regions tested positive for *Anaplasma* spp. by 16S rRNA PCR. The infection rate in goats and sheep was determined as 28% and 27.7%, respectively. RFLP revealed the prevalence of *A. phagocytophilum* and *A. phagocytophilum*-like 1 as 1.4% and 26.5%, respectively. No PCR amplicons derived from goats and sheep were digested by the *Bsa*I enzyme, confirming the absence of Chinese variant (*A. phagocytophilum*-like 2). Of the 121 positive samples with 16S rRNA PCR, 110 (95.6%) were positive with groEL nested PCR. Six of them (6/110, 5.4%) were positive for *A. phagocytophilum* and 104 (94.5%) were positive for *A. phagocytophilum*like 1 (Table 2).

Association of the frequency of *A. phagocytophilum* and *A. phagocytophilum*-like 1 variant in small ruminants with species and tick infestation is documented in Table 3. The prevalence of *Anaplasma* spp. was comparable in species, and no difference was detected between infection rates in sheep and goats (*p* = 0.9603). However, the prevalence significantly increased with tick infestation in small ruminants (*p* = 0.0003) (Table 3).


**Table 2.** Samples origin, 16S rRNA PCR, RFLP and groEL PCR.

**Table 3.** Association of the frequency (16S rRNA PCR) of *Anaplasma phagocytophilum* and related variants in small ruminants with species and tick infestation.


#### *3.3. Molecular and Phylogenetic Analyses*

To validate the RFLP results and identify genetic variants of *A. phagocytophilum*like 1, randomly selected 10 representative samples were sequenced. The sequences shared 100% identity to each other. Therefore, one representative sequence for *A. phagocytophilum*-like 1 was submitted to the NCBI GenBank database, and deposited with accession number MT881655. This finding indicated that one variant was identified, and named as Aplike1OvineCaprine in this work. BlastN analysis demonstrated that the Aplike1OvineCaprine variant indicated high similarity (99–100%) to those *Anaplasma* isolates deposited in the GenBank as uncultured *Anaplasma* sp. and *A. phagocytophilum*. Moreover, the Aplike1OvineCaprine variant was 100% identical to those of *A. phagocytophilum*-like 1 detected in sheep (Aplike1Ov1, KX702978) and goat (Aplike1GGo2, KM285227) from Tunisia, and cattle from Turkey (Aplike1Bv, MT338494) (Table 4). The *A. phagocytophilum* Akseki11 Sheep Turkey isolate obtained in this study shared 99.3–99.6% identity isolated from *Niviventer confucianus* (*A. phagocytophilum* ZJ-HGA strain, DQ458805) and human (*A. phagocytophilum* HZ strain, NR\_074113), respectively.

Phylogenetic analysis using the *16S rRNA* gene showed that our variant (Aplike1Ovine-Caprine) clustered a distinct group with those of *A. phagocytophilum*-like 1 previously published sequences reported in sheep, goats, cattle, deer, and *Haemaphysalis ginghaiensis* (Figure 2).



of 1137 nucleotide indicates the substitution of A by G allowing for the distinction between *A. phagocytophilum* and related variants (like 1 and 2) by *Xcm*I enzyme [8].

*haiensis* (Figure 2).

Phylogenetic analysis using the *16S rRNA* gene showed that our variant (Aplike1OvineCaprine) clustered a distinct group with those of *A. phagocytophilum*-like 1 previously published sequences reported in sheep, goats, cattle, deer, and *Haemaphysalis ging-*

**Figure 2.** Maximum likelihood phylogenetic tree was inferred with partial sequences (598–599 bp) of the 16S rRNA gene of *Anaplasma* sp. related to *A. phagocytophilum* isolated from sheep in Turkey (highlighted in red) with other *Anaplasma* spp. retrieved from GenBank. Numbers at the nodes refer percentage occurrence in 1000 the bootstrap replication. The new sequences of *A. phagocytophilum*-like 1 and *A. phagocytophilum* from this study were highlighted in red. *Ehrlichia ruminantium* was used as an outgroup. **Figure 2.** Maximum likelihood phylogenetic tree was inferred with partial sequences (598–599 bp) of the 16S rRNA gene of *Anaplasma* sp. related to *A. phagocytophilum* isolated from sheep in Turkey (highlighted in red) with other *Anaplasma* spp. retrieved from GenBank. Numbers at the nodes refer percentage occurrence in 1000 the bootstrap replication. The new sequences of *A. phagocytophilum*-like 1 and *A. phagocytophilum* from this study were highlighted in red. *Ehrlichia ruminantium* was used as an outgroup.

#### **4. Discussion**

*Anaplasma phagocytophilum* causes tick-borne fever in small ruminants and granulocytic anaplasmosis in horses and dogs [1,2]. It is an emerging tick-borne pathogen for humans as well [3]. The genetic diversity of *A. phagocytophilum* is much greater than expected. Indeed, recent studies have revealed the existence of two distinct *Anaplasma* species or variants related to *A. phagocytophilum,* one in Japan and the other in China [11–14,24]. Then, these pathogens were designated as *A. phagocytophilum*-like 1 and *A. phagocytophilum*-like 2 variants [7,8]. More recently, both genotypes have been documented in ruminants and *R. turanicus* in Tunisia [7,8,17,29], cattle in South Korea [16], and small ruminants in Italy [17]. In the present study, a survey was carried out to detect and identify *A. phagocytophilum* and *A. phagocytophilum*-like variants in small ruminants from the Mediterranean region of Turkey. Our findings provide molecular evidence for the presence of *A. phagocytophilum* and *A. phagocytophilum*-like 1 in sampled sheep and goats. In the previous studies carried out in Turkey, *A. phagocytophilum* has been reported in small ruminants [18,21,30]. However, this is the first time that *A. phagocytophilum*-like 1 variant in sheep and goats have been reported in the country.

Contrary to *A. phagocytophilum*, it has been suggested that both Japanese and Chinese variants do not cause clinical infection in ruminants [8,17]. In this study, a high prevalence for *A. phagocytophilum*-like 1 variant was determined (26.5%), but no clinical infection for tick-borne fever was observed in sheep and goats during sample collection. This result is consistent with the previous suggestions that *A. phagocytophilum*-like variants are considered non-pathogenic for ruminants [8,16,17]. The prevalence of *A. phagocytophilum*like 1 (26.5%) in small ruminants obtained in this study was higher than that observed in Tunisian sheep (7%) and goats (13.1%) [8], however, it was lower than that observed in other studies conducted in Mediterranean small ruminants (122/203, 60%) from Tunisia and Italy [17].

It has been previously suggested that serological cross-reactions occur between *A. phagocytophilum* and other *Anaplasma* species [31,32]. The same situation may be true in some circumstances for molecular markers, for example a pair of primers (SSAP2f/SSAP2r) based on the *16S rRNA* gene of *A. phagocytophilum* were designed for the specific amplification [10]. However, it has been shown that these primers also detect distinct *Anaplasma* variants related to *A. phagocytophilum* [7,8,24]. In this work, the frequency of pathogenic *A. phagocytophilum* was 1.4%, which is not consistent with the previous studies in Turkey that reported values of 66.7% in Central Anatolia [30] and 19.7% in Eastern Anatolia [21]. The high infection rates obtained in the previous studies may be due to the selected primers for the amplification of *A. phagocytophilum*. EE1/EE2 and SSAP2f/SSAP2r primers have been selected to detect *A. phagocytophilum* in the studies conducted in Central Anatolia [30] and Eastern Anatolia [21], respectively. However, the EE1/EE2 primers are universal for the detection of all *Anaplasma* spp. including *A. phagocytophilum*-like variants [33]. It has been also reported that the SSAP2f/SSAP2r can amplify not only *A. phagocytophilum*, but also *A. phagocytophilum*-like variants [7,8,24]. This study provides molecular data for the circulation of *A. phagocytophilum* and *A. phagocytophilum*-like 1 Turkish small ruminants. Therefore, cross-reactivity between *A. phagocytophilum* and related variant should be considered in interpreting the findings of surveys to be carried out in the area, where *A. phagocytophilum* and *A. phagocytophilum*-like variant co-exist.

As several domestic and wild mammals are hosts or reservoirs for *A. phagocytophilum* [1,2], abundance and intensity of the tick vector, *I. ricinus* in Europe including Turkey are considered a major factor affecting the distribution of the pathogen in a specific area. It is well known that there is no *I. ricinus* in the Eastern and Central Anatolian regions of Turkey [34]. It has been reported that *A. phagocytophilum* is transmitted by *Ixodes* spp. (*I. persulcatus*, *I. scapularis* and *I. ricinus*) in some parts of the world including in Europe [1,35]. In Turkey, *I. ricinus* collected from humans were positive for *A. phagocytophilum* [5]. So far, data on the transmission of *A. phagocytophilum*-like variants by ticks are lacking. A recent study reported that *R. turanicus* was common in sampled sheep and goats in Tunisia, and

one *R. turanicus* tick feeding on the goat was found to be infected with *A. phagocytophilum*like 2 [28]. In the present study, potential vectors of *A. phagocytophilum*-like 1 was not studied, but we found that the sampled sheep and goats were commonly infested with *R. bursa* (86%), *R. turanicus* (6.6%), *D. marginatus* (6.4%), and very rarely *I. ricinus* (0.06%, only one specimen). Our finding also showed a correlation between *Anaplasma* positivity and the presence of ticks (*p* = 0.0003), compatible with the finding that the prevalence of *A. phagocytophilum*-like 1 was higher in goats infested by ticks than in not infested [7]. Based on the abundance of *Rhipicepahlus* and *Dermacentor* ticks and the very rarity of *I. ricinus*, we can assume that *Rhipicephalus* and *Dermacentor* may play an important role in the transmission of *A. phagocytophilum*-like 1 rather than *I. ricinus*. This assumption is supported by the previous findings that a high prevalence of *A. phagocytophilum*-like variants have been reported in ruminants in the higher semi-arid area of Tunisia, where *I. ricinus* is not present [8]. However, more detailed studies are needed to validate this assumption and to establish what tick species may play a role in the transmission of *A. phagocytophilum*-like 1 in Turkey.

Our sequencing validated RFLP findings, and showed that the sampled small ruminants were found to be infected with *A. phagocytophilum*-like 1. Phylogenetic analysis indicated two main separate branches. The Aplike1OvineCaprine (MT881655) variant obtained in this study, as well as those previously reported from sheep, goats, cattle, and ticks, formed a monophylogenic clade distinct from *A. phagocytophilum* and *A. phagocytophilum*-like 2, and other members of *Anaplasma* spp. [7,8,24].

#### **5. Conclusions**

This work provides molecular data for the circulation of *A. phagocytophilum-*like 1 for the first time in Turkey. The novel strain is widespread in small ruminants in the Mediterranean area of Turkey with an overall prevalence of 26.5%. This finding revealed that the variant should be considered in the diagnosis of caprine and ovine anaplasmosis.

**Author Contributions:** Conceptualization, project administration, methodology, validation, formal analysis, investigation, writing—review and editing, M.A.; methodology, validation, formal analysis, S.Ö.; validation, formal analysis, M.C.U. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by founding from the Scientific and Technological Research Council of Turkey (TUBITAK) (project no. 118O871).

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by Firat University Animal Experiment Local Ethics Committee (2018/100, 109. 11. 2018).

**Data Availability Statement:** Data available in a publicly accessible repository.

**Acknowledgments:** The authors are grateful to Alberto Alberti (University of Sassari, Sassari, Italy) for positive control DNAs from *A. phagocytophilum*-like variants. We are also thankful to the farmers for their cooperation and the veterinary clinician's Ali Asar for helping with sample collection.

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

#### **References**


*Article*
