Expanding the Universe of Hemoplasmas: Multi-Locus Sequencing Reveals Putative Novel Hemoplasmas in Lowland Tapirs (Tapirus terrestris), the Largest Land Mammals in Brazil

The lowland tapir (Tapirus terrestris) is the largest land mammal in Brazil and classified as a vulnerable species, according to the assessment of the risk of extinction. The present study aimed at investigating the occurrence and genetic diversity of hemoplasmas in free-ranging T. terrestris from the Brazilian Pantanal and Cerrado biomes. Blood samples were collected from 94 living and eight road-killed tapirs, totalizing 125 samples Conventional PCR targeting four different genes (16S rRNA, 23S rRNA, RNAse P, and dnaK) were performed, and the obtained sequences were submitted for phylogenetic, genotype diversity, and distance analyses. The association between hemoplasma positivity and possible risk variables (age, gender, and origin) was assessed. Out of 122 analyzed samples, 41 (41/122; 33.61% CI: 25.84–42.38%) were positive in the 16S rRNA-based PCR assay for hemoplasmas. Positivity for hemoplasmas did not differ between tapirs’ gender and age. Tapirs from Pantanal were 5.64 times more likely to present positive results for hemoplasmas when compared to tapirs sampled in Cerrado. BLASTn, phylogenetic, genotype diversity, and distance analyses performed herein showed that the sampled lowland tapirs might be infected by two genetically distinct hemoplasmas, namely ‘Candidatus Mycoplasma haematoterrestris’ and ‘Candidatus Mycoplasma haematotapirus’. While the former was positioned into “Mycoplasma haemofelis group” and closely related to ‘Candidatus Mycoplasma haematoparvum, the latter was positioned into “Mycoplasma suis group” and closely related to ‘Candidatus Mycoplasma haematobos’. The impact of both putative novel species on tapir health status should be investigated.


Introduction
The lowland tapir (Tapirus terrestris) is the largest land mammal in Brazil and the species is widely distributed throughout the country. Although the lowland tapir can still be found in four Brazilian biomes (Atlantic Forest, Pantanal wetlands, Amazon, and Cerrado), this species is classified as vulnerable according to extinction risk assessments [1,2].
The emergence of infectious diseases in tapir populations may occur as a consequence of habitat loss and increased contact between animals. In this instance, preventive medicine should be applied for both in situ and ex situ conservation attempts [3]. In the last decades, studies focusing on direct detection of different pathogens, such as Trypanosoma deer (Blastocerus dichotomus), crowned solitary eagle (Buteogallus coronatus), and hyacinth macaw (Anodorhynchus hyacinthinus). The impacts caused by wildfires are related as direct injuries and death, or indirect injuries caused by habitat loss and resource depletion [35].

Sampling
This study was approved by the Ethics Committee for Animal Experimentation of FCAV/UNESP (Faculty of Agricultural and Veterinary Sciences of the São Paulo State University) under protocol number 4558/20. The "Instituto Chico Mendes de Conservação da Biodiversidade (ICMBIO)" provided the required annual permits for the capture and immobilization of tapirs and collection of biological samples (SISBIO# 14,603). All protocols for the capture, anesthesia, handling, and sampling of tapirs have been reviewed and approved by the Veterinary Advisors of the Association of Zoos and Aquariums (AZA)-Tapir Taxon Advisory Group (TAG), and the Veterinary Committee of the IUCN SSC Tapir Specialist Group (TSG).
From 2013 to 2018, blood samples from free-ranging T. terrestris were collected from 94 living and 8 road-killed individuals, totalizing 125 samples. From this total, 78 ( Sampling was performed during tapir anesthesia for the installation of GPS collars by professionals from the "Iniciativa Nacional para a Conservação da Anta Brasileira (INCAB-IPÊ)" (Lowland Tapir Conservation Initiative (LTCI-IPÊ)). Some living animals (n = 20) were sampled more than once. Moreover, eight samples (8/125; 6.40% CI: 3.28-12.12%) belonging to eight road-killed animals on highways BR-267 and MS-040 in the Cerrado biome were collected during necropsy procedures.
The sampled living animals included females (

DNA Extraction and PCR Protocols for Mammals' Endogenous Genes
DNA extraction was performed using a commercial kit (InstaGene™ Matrix, Biorad ® , Hercules, CA, USA) and following the manufacturers' instructions. In order to ensure successful DNA extraction, a conventional PCR for the mammal-endogenous gene glyceraldehyde-3-phosphate dehydrogenase (gapdh) targeting a 450 bp fragment [36] was performed in all samples. Samples that were not successfully amplified by the gapdhbased PCR protocol were submitted to a second protocol targeting a 227 bp fragment of the irpb (interphotoreceptor retinoid-binding protein) gene [37]. Samples that did not yield amplicons in either of the used PCR protocols were excluded from the subsequent analysis.

Conventional Polymerase Chain Reaction (cPCR) Assays for Hemoplasmas Based on the 16S rRNA, 23S rRNA, RNAse P and dnaK Gene Fragments
Samples were screened for the presence of hemoplasmas DNA using a semi-nested PCR protocol targeting a fragment of approximately 1107 bp of the 16S rRNA gene [38,39]. Positive samples were submitted to PCR protocols targeting the 23S rRNA (approximately 800 bp) [18], RNAse P (approximately 165 bp) [16], and dnaK (approximately 544 bp) genes [40] for additional molecular characterization. Thermal conditions and PCR reagent concentrations were slightly modified from the originally published protocols for RNAse P Microorganisms 2022, 10, 614 4 of 22 and dnaK genes (Table 1). DNA obtained from a naturally infected sheep with Mycoplasma ovis [18] and ultrapure RNAse and DNAse-free water (Promega, Madison, WI, USA) were used as positive and negative controls, respectively, in all PCR assays for hemoplasmas. The products obtained in PCR assays were separated by electrophoresis on a 1% agarose gel stained with ethidium bromide (Life Technologies™, Carlsbad, CA, USA) at 100 V/150 mA for 50 min. The gels were imaged under ultraviolet light (ChemiDoc MP Imaging System, Bio Rad™, Hercules, CA, USA) using the Image Lab Software v4.1 (Biorad, Hercules, CA, USA).

Genetic Diversity Assessment
The genetic diversity assessment was performed using only the obtained 16S rRNA sequences. In order to calculate nucleotide diversity (π), polymorphism level (haplotype diversity-[dh]), number of haplotypes (h), and the average number of nucleotide differences (K) among the sequences obtained, the DnaSP program v5 (http://www.ub.edu/ dnasp/, accessed on 21 December 2021) [48] was used. The Genotype network was constructed in PopART (http://popart.otago.ac.nz, accessed on 11 December 2021), using the TCS inference method [49]. GPS coordinates collected for each sampling location were used to develop a map representing the genotype distribution using PopART software. Additionally, a distance-based analysis was performed using SplitsTree v4.14.6 (University of Tubingen, Tubingen, Germany) [50] to investigate the genetic relationship among hemoplasma genotypes detected in the present study and those previously deposited in GenBank. The pairwise distance matrix from the alignment of the 16S rRNA sequences detected in tapirs in the current study was calculated using the p-distance model and included Transitions + Transversions substitutions with uniform rates. Analysis was performed using MEGAX software (https://www.megasoftware.net, accessed on 17 December 2021) [51,52]. Data were transferred to a Microsoft Excel 2016 sheet to construct the heat map according to the rates obtained on the pairwise distance matrix.

Statistical Analysis
The chi-square test was used to determine associations between variables (gender, location, and age) and the outcomes (positive or negative PCR results for hemoplasmas). Odds ratio (OR), 95% confidence interval, and p-values were calculated for each variable. Results considered significantly different when p < 0.05. Data were compiled and analyzed in Epi Info™ software (v7.1.5, Centers for Disease Control and Prevention, Atlanta, GA, USA).

PCR for Mammals' Endogenous Genes
Out of 125 blood samples tested for the presence of the gapdh gene fragment, 9 (9/125; 07.20% CI: 03.83-13.12%) did not amplify fragments of the expected size and were submitted to the irpb-based PCR protocol. Out of 9 analyzed samples, 3 (3/125; 02.40%; CI: 00.82-06.82%) did not successfully amplify fragments of the expected size and were excluded from the following analyses. All excluded samples were obtained from road-killed tapirs.

PCR Assays for Hemoplasmas
From 122 analyzed samples, 41 ( .12%) sub-adults. Some animals (n = 3) presented positive samples in more than one sampling, totalizing seven positive samples (mean value: 2.3 positive samples/tapir) from the same animals. The two samples from road-killed animals that yielded positive results for the partial 16S rRNA gene were obtained from a sub-adult female and an adult male. In total, 35.29% (36/102; CI: 26.71-44.95%) tapirs (living or road-killed) presented positive results for hemotropic Mycoplasma sp. Results regarding positive samples, sampling dates, age/gender/location from tapirs, and GenBank accession numbers were summarized in Table 2.  14%) submitted to dnaK-based PCR protocol presented expected band sizes in agarose gel electrophoresis and one was successfully sequenced and deposited in the Gen-Bank database (OM339521) ( Table 2). Only one sample (ID: MIA-P) presented sequences for all four targeted genes. The BLASTn analysis results for each sequence obtained by amplifying 23S rRNA, RNAse P, and dnaK genes in the present study are shown in Tables 3-5.

Phylogenetic Inference
Phylogenetic trees were constructed for four partial gene fragments by Bayesian inference using 10 7 generations of MCMC (Monte Carlo Markov Chains) with two independent runs and 10% of burn-in (Figures 1-4). For the partial 16S rRNA gene analysis (Figure 1), a total size alignment of 909 bp was constructed using 99 homolog sequences and one outgroup (AB042061).
The best-fit model for this gene analysis was determined as F81+G. For the partial 23S rRNA gene analysis (Figure 2), an alignment with a total size of 720 bp was constructed using 41 homolog sequences and one outgroup (NR103037). The best-fit model for this gene analysis was determined as F81+G. For the partial RNAse P gene analysis (Figure 3), a total size alignment of 230 bp was constructed using 23 homolog sequences and one outgroup (U64878). The best-fit model for this gene analysis was determined as F81+I+G. For the partial dnaK gene analysis (Figure 4) a total size alignment of 622 bp was constructed using 19 homolog sequences and one outgroup (KJ690086). The best-fit model for this gene analysis was determined as F81+G.
Regarding the 16S rRNA phylogenetic analysis, our sequences were divided into two different major clades. The first clade was also sub-divided into three minor clades and positioned in the "Mycoplasma suis group", whereas the second clade comprised of only three sequences and was distantly positioned in the "Mycoplasma haemofelis group". These results suggest the occurrence of two distinct species occurring in sampled tapirs. In order to facilitate the understanding from now on these two clades will be identified as Ca1 and Ca2, respectively. Sequences that fitted in Ca1 on partial 16S rRNA-based phylogeny followed a similar pattern in other genes phylogenies, except for the RNAse P gene. In the phylogenetic tree based on the RNAse P gene, both obtained sequences were positioned in the "M. suis group", albeit each one was positioned in a different group on 16S rRNA phylogeny. Moreover, one sample (ID: AA-P) that was positioned in Ca1 on 16S rRNA-based phylogeny was found in Ca2 in the phylogenetic inference based on the partial 23S rRNA.

Genetic Diversity Assessment
Regarding the 16S rRNA genotype analysis, 22 different genotypes were identified among 30 sequences ( Figure 5). Values of nucleotide diversity (π), polymorphism level (haplotype diversity-[dh]), number of haplotypes (h), average number of nucleotide differences (K) among the sequences obtained, and the number of mutations between sequences from Ca1 and Ca2 are described in Table 6. 16S rRNA genotypes and sampling location of each obtained sequence are shown in Supplementary Table S1. The map representing genotype distributions along the sampling area reinforced the geographical distance between genotypes identified in Pantanal and Cerrado areas (Figures 6-8).
lowed a similar pattern in other genes phylogenies, except for the RNAse P gene. In the phylogenetic tree based on the RNAse P gene, both obtained sequences were positioned in the "M. suis group", albeit each one was positioned in a different group on 16S rRNA phylogeny. Moreover, one sample (ID: AA-P) that was positioned in Ca1 on 16S rRNAbased phylogeny was found in Ca2 in the phylogenetic inference based on the partial 23S rRNA.

Genetic Diversity Assessment
Regarding the 16S rRNA genotype analysis, 22 different genotypes were identified among 30 sequences ( Figure 5). Values of nucleotide diversity (π), polymorphism level (haplotype diversity -[dh]), number of haplotypes (h), average number of nucleotide differences (K) among the sequences obtained, and the number of mutations between sequences from Ca1 and Ca2 are described in Table 6. 16S rRNA genotypes and sampling location of each obtained sequence are shown in Supplementary Table S1. The map representing genotype distributions along the sampling area reinforced the geographical distance between genotypes identified in Pantanal and Cerrado areas (Figures 6-8).

Figure 5.
Genotype diversity among 16S rRNA gene sequences detected herein. Analysis was made using DnaSP6. Inference and graphic representation were made by TCS network method on PopART software. Genotypes in blue were obtained from samples from tapirs in Pantanal regions meanwhile genotypes in green were obtained from samples from tapirs in Cerrado regions. Table 6. Values obtained regarding genotype diversity by DnaSP software and based on partial hemoplasma 16S rRNA sequences detected in tapirs from the present study. Figure 5. Genotype diversity among 16S rRNA gene sequences detected herein. Analysis was made using DnaSP6. Inference and graphic representation were made by TCS network method on PopART software. Genotypes in blue were obtained from samples from tapirs in Pantanal regions meanwhile genotypes in green were obtained from samples from tapirs in Cerrado regions.  Figure 6. Representation of the genotype distribution along all sampling areas. Map was constructed using PopART software based on GPS coordinates data of each sampling.    Figure 6. Representation of the genotype distribution along all sampling areas. Map was constructed using PopART software based on GPS coordinates data of each sampling.

Distance Analysis by SplitsTree
The distance analysis performed by the neighbor-joining method and displayed by SplitsTree v4.14.6 demonstrated that Ca1 and Ca2 sequences were distinctly disposed among other hemotropic Mycoplasma species and Candidatus (Figure 9).

Distance Analysis by SplitsTree
The distance analysis performed by the neighbor-joining method and displayed by SplitsTree v4.14.6 demonstrated that Ca1 and Ca2 sequences were distinctly disposed among other hemotropic Mycoplasma species and Candidatus (Figure 9).

Distance Analysis by SplitsTree
The distance analysis performed by the neighbor-joining method and displayed by SplitsTree v4.14.6 demonstrated that Ca1 and Ca2 sequences were distinctly disposed among other hemotropic Mycoplasma species and Candidatus (Figure 9).

Distance Matrix Analysis
The heat map with the distance matrix from all partial 16S rRNA sequences detected herein corroborated with other achieved results from the present study and demonstrated a marked difference between sequences from Ca1 and Ca2 ( Figure 10). The minimum and maximum divergence percentage rates found for all the sequences that fit in Ca1 were 0.0% and 2.92%, respectively. When comparing to the phylogenetic closest sequences (FJ004275, KM275255, MN543630, AY297712, MH094850, EF416569, GQ129112, KF366443, GQ129114, AY532390, AY383241), the minimum and maximum divergence found between Ca1 members and those sequences was 1.52% and 3.27%, respectively. When comparing the Ca1 members and Ca2 members, the minimum and maximum divergence found was 15.23% and 16.22%, respectively. The minimum and maximum divergence percentage rates found for all members from Ca2 were 0.0%. When comparing to the phylogenetic closest sequences (FN421445, MG948628, MG948631, MH388480, MW463059, EF424082, EF616467, EF460765, HFU95297, GQ129115, KM275247, KP715860), the minimum and maximum divergence found between Ca2 members and those sequences was 5.12% and 7.23%, respectively. The distance matrix is available in Supplementary Table S2.

Distance Matrix Analysis
The heat map with the distance matrix from all partial 16S rRNA sequences detected herein corroborated with other achieved results from the present study and demonstrated a marked difference between sequences from Ca1 and Ca2 ( Figure 10). The minimum and maximum divergence percentage rates found for all the sequences that fit in Ca1 were 0.0% and 2.92%, respectively. When comparing to the phylogenetic closest sequences (FJ004275, KM275255, MN543630, AY297712, MH094850, EF416569, GQ129112, KF366443, GQ129114, AY532390, AY383241), the minimum and maximum divergence found between Ca1 members and those sequences was 1.52% and 3.27%, respectively. When comparing the Ca1 members and Ca2 members, the minimum and maximum divergence found was 15.23% and 16.22%, respectively. The minimum and maximum divergence percentage rates found for all members from Ca2 were 0.0%. When comparing to the phylogenetic closest sequences (FN421445, MG948628, MG948631, MH388480, MW463059, EF424082, EF616467, EF460765, HFU95297, GQ129115, KM275247, KP715860), the minimum and maximum divergence found between Ca2 members and those sequences was 5.12% and 7.23%, respectively. The distance matrix is available in Supplementary Table  S2.

Statistical Analysis
Statistical analysis results including OR (odds ratio) and p-value are summarized in Table 7. Regarding the analysis between location and positivity for hemoplasmas, the number of positive tapirs sampled in Pantanal wetlands was statistically higher when compared to those sampled in the Cerrado region (p-value = 0.0001). The other analyzed variables (gender and age) did not present statistically significant differences for this outcome (p-value > 0,05). The OR value demonstrated that tapirs from Pantanal were 5.64 times more likely to present positive results for the PCR protocol targeting the partial Mycoplasma spp. 16S rRNA gene used in the present study when compared to tapirs sampled in Cerrado.

Statistical Analysis
Statistical analysis results including OR (odds ratio) and p-value are summarized in Table 7. Regarding the analysis between location and positivity for hemoplasmas, the number of positive tapirs sampled in Pantanal wetlands was statistically higher when compared to those sampled in the Cerrado region (p-value = 0.0001). The other analyzed variables (gender and age) did not present statistically significant differences for this outcome (p-value > 0.05). The OR value demonstrated that tapirs from Pantanal were 5.64 times more likely to present positive results for the PCR protocol targeting the partial Mycoplasma spp. 16S rRNA gene used in the present study when compared to tapirs sampled in Cerrado.

Discussion
We described herein the occurrence of at least two genetically distinct hemoplasma species occurring in free-ranging tapirs. Reports of Candidatus and novel species of these bacteria in Brazilian wild fauna are becoming common due to the growing use of molecular and genetic analysis tools. Based mainly on the partial 16S rRNA and 23S rRNA genes amplification, novel species of hemotropic Mycoplasma spp. infecting capybaras [22], hairy dwarf porcupines [23], opossums [24], and coatis [25] have been proposed. When analyzing the topologies found by partial 16S rRNA gene phylogeny from these reports, sequences from each novel Candidatus species were positioned separately in clades inside the "Mycoplasma haemofelis group", with satisfactory bootstrap values for species separation. In the present work, we detected different sequences from tapirs that fit in both the "Mycoplasma haemofelis group" (Ca2) and "Mycoplasma suis group" (Ca1) by partial 16S rRNA, 23S rRNA, and dnaK genes-based phylogeny, which was unprecedented for putative novel species of hemoplasmas reported in wild hosts from Brazil.
The 16S rRNA is considered the "gold standard" target gene for PCR assays aiming at detecting and identifying hemoplasmas and different pairs of primers have been used for its purpose [15]. Although this gene may be considered highly conserved, the 16S rRNA gene from bacterial species may exhibit considerable variations even on supposed conserved regions [53], which allow species differentiation. In the present study, all partial 16S rRNA sequences obtained for the Ca1 genotypes presented a similarity range of 97.20-98.20% with 'Candidatus Mycoplasma haematoparvum' detected in dogs from Italy and the USA (MH094850, AY383241) by BLASTn analysis. Meanwhile, the Ca2 genotypes demonstrated a similarity range of 94.99% with 'Candidatus Mycoplasma haematobos' detected in cattle from Cuba (MG948631). The p-distance values obtained by the distance matrix analysis showed that Ca1 members presented divergence rates of 0.0-2.92%, whereas the same comparison using the closest phylogenetic group and members of Ca1 presented rates of 1.52-3.27%. When using the 16S rRNA gene, a similarity rate of at least <97% is expected for bacteria belonging to different species. However, divergence rates <3% do not necessarily indicate that sequences belong to the same species [54]. Once phylogenetic analyses based on four molecular markers (16S rRNA, 23S rRNA, RNAse P, and dnaK) and 16S-based distance analyses strongly supported the differentiation between Ca1 and the closest hemoplasma species ('Ca. M. haematoparvum'), it is likely that Ca1 represents a genetically distinct group of hemoplasmas.
The use of different genes for the molecular characterization of hemoplasmas is reported. Amplification of partial 23S rRNA gene was reported in phylogenetic studies of hemoplasmas identified in captive cervids from southern Brazil [55] and for the phylogenetic study of Mycoplasma ovis infecting sheep from the same region [18]. Recently, this gene has also been used for the characterization of novel Candidatus species of hemoplasmas from wild animals in Brazil [22][23][24][25]. In the present study, partial 23S rRNA sequences that fit in Ca1 presented similarity percentages ranging from 89.59-91.50% with a sequence of 'Candidatus Mycoplasma haematominutum' from the UK (HE613254). Meanwhile, 23S rRNA sequences that fit in Ca2 presented similarity ranges of 90.20-90.57% with a sequence of Mycoplasma haemofelis from the USA (NR103993). In fact, putative novel hemoplasma species may present lower similarity rates in BLASTn analysis for 23S rRNA sequences when compared to 16S rRNA sequences [22,55,56]. These differences may be explained by the fact that, although the 23S rRNA is considered as phylogenetically conserved as 16S rRNA, it presents a higher degree of sequence variability [57].
The RNAse P (rnpB) gene codifies the RNA subunit of endoribonuclease P with a length of approximately 400 bp [58]. Although only a few numbers of studies used phylogenetic trees of hemoplasmas using fragments of the RNAse P gene, our phylogeny analysis using this gene agreed with these studies, showing a separation of species among the "Mycoplasma haemofelis group" and "Mycoplasma suis/'Ca. M. haematominutum' group" [16,58]. However, both sequences obtained in the present study (OM317758-OM317759) were positioned in the "Mycoplasma haemofelis group", albeit these same two samples were positioned separately in the "Mycoplasma haemofelis group" and "Mycoplasma suis group" in all other analyzed targeted genes. Although it is expected that phylogenies using RNAse P fragments of hemoplasmas may present similar topologies to those using 16S rRNA [58], the fragments obtained herein are very short (~100 bp), precluding robust species differentiation.
Amplification and subsequent phylogenetic analysis of the dnaK gene have been demonstrated to be a useful tool for the separation of hemoplasmas from other Mollicute species [59] and also for the separation of subspecies of hemoplasmas [60]. This gene is responsible for coding a chaperon protein (heat shock protein 70) and is considered a great genetic marker for species differentiation, since it may present more variable regions within the sequences when compared to ribosomal fragments [59]. Unfortunately, only one sequence from the present study was successfully sequenced using the dnaK (OM339521) PCR protocol, and comparison with other sequences obtained from tapirs was not possible. However, this single sequence followed a similar topology pattern of 16S rRNA and 23S rRNA phylogenies, since it fitted on Ca1 in the "Mycoplasma suis group", reinforcing the phylogenetic position found for this putative novel hemoplasma species.
A high diversity of genotypes (n = 22) was found in the present work by analysis of 30 partial 16S rRNA sequences of hemoplasmas obtained from tapirs' samples. Richness on genotype diversity was also reported for hemoplasmas detected in bats from Brazil [56]. Although the genotypes found in tapirs from the Pantanal diverged from those found in tapirs from the Cerrado, a genetic proximity between these genotypes was demonstrated on the genotype network analysis. The occurrence of more than one genotype in the same region infecting wild hosts has already been reported in bats from Brazil [56]. Interestingly, one animal that was sampled at three different times (IDs: WE-P-1, WE-P-2, and WE-P-3) presented three different genotypes according to the time of sampling. Moreover, one animal (ID: AA-P), that was sampled only once, was positioned in Ca1 on 16S rRNA-based phylogeny and in Ca2 in the phylogenetic inference based on the partial 23S rRNA. These data suggest that tapirs may be susceptible to co-infections with different hemoplasma genotypes or species, with these co-infections occurring simultaneously or at different times.
Some animals that were sampled more than once presented positive samples in more than one sampling (IDs: MA-P-1, MA-P-2, JE-P-1, JE-P-2, WEP-1, WEP-2, WE-P3). Meanwhile, some animals that also were sampled more than once presented positive results for one sample only (IDs: TD-P-1, VA-P-1, ANO-C-2, CNA-C-2, SO-C-2, FFO-P-2, JO-P-2, DO-P-2, TD-P-3, SAO-P-2). These results may suggest either that tapirs may be maintained as chronically infected hosts or that hemoplasma bacteremia may be too low, precluding molecular detection. Chronic infection by hemoplasmas has been commonly reported in domestic species, such as pigs [61], cattle [62], and cats [63,64]. Although it is not possible to extend these findings for wildlife, chronic infection by hemoplasmas is appointed as a cause of unthriftiness in newborn piglets [65] and lower calf birth weight in cattle [62] which may raise a red flag in the context of species conservation.
Herein, gender and age were not associated with hemoplasma infection in the sampled tapirs. Although tapirs from Cerrado presented more health abnormalities when compared to tapirs from other biomes [65], we found that animals sampled in the Pantanal were 5.64 times more prone to be infected by hemoplasmas when compared to tapirs sampled in the Cerrado. In cats, population density may be associated as a risk factor for hemoplasma infection [66]. Some Pantanal sites and habitats were already reported as able to sustain high densities of tapir populations when these animals are not exposed to anthropic actions [67].
Tapirs from both Cerrado and Pantanal regions were found parasitized by Amblyomma spp. and Rhipicephalus microplus ticks. Although information regarding the transmission of hemoplasmas by ticks and other vectors is still lacking, R. microplus ticks were appointed as capable of transmitting 'Ca. M. haematobos' for egg and larval stages and also to transmit this hemoplasma species for mice [68]. Future studies aiming at investigating the role of ticks in the transmission of hemoplasmas among tapirs are needed.
Mycoplasma haemofelis and 'Ca. M. haematominutum' DNA was detected in saliva and salivary glands of infected cats, suggesting that social interactions, such as aggression, may be related to hemoplasma infections [69]. Lowland tapirs were considered to be as solitary and show tolerance to other individuals not influenced by kinship [70]. More studies are necessary to elucidate the transmission routes for hemoplasmas among free-ranging lowland tapirs in Brazil.
Besides the assessment, for the first time, of the occurrence of hemotropic Mycoplasma among tapirs in two distinct Brazilian biomes, the present work showed high genetic diversity, and, at least, two genetically distinct hemoplasma species infecting these mammals. Accordingly, the obtained results reinforce the need for multi-locus and large-scale sequencing aiming at unraveling accurately the genetic diversity of hemoplasmas in wild animals.

Conclusions
At least two genetically distinct species of hemotropic Mycoplasma spp. occurs in free-ranging T. terrestris from the Pantanal and Cerrado regions in Brazil. The occurrence of hemoplasmas did not differ according to the gender or age of the sampled tapirs. Animals sampled in the Pantanal may be at a higher risk of becoming infected by hemoplasma when compared to those in the Cerrado. We propose that the two genetically divergent species found infecting tapirs from the present study represent putative novel Candidatus species and the names 'Candidatus Mycoplasma haematoterrestris' and 'Candidatus Mycoplasma haematotapirus' are proposed for the species found in Ca1 and Ca2, respectively.