Assessment of bacterial diversity of Rhipicephalus microplus ticks from two livestock agroecosystems in Antioquia, Colombia

Rhipicephalus microplus is recognized as a tick species highly prevalent in cattle, with a wide pantropical distribution that seems to continue spreading geographically. However, its role as a biological vector has been scarcely studied in the livestock context. In this study, a 16S rRNA next-generation sequencing analysis was used to determine bacterial diversity in salivary glands and gut of R. microplus from two contrasting livestock agroecosystems in Antioquia, Colombia. Both the culture-independent approach (CI) and the culture-dependent (CD) approach were complementarily adopted in this study. A total of 341 unique OTUs were assigned, the richness showed to be higher in the Northern than in the Middle Magdalena region, and a high diversity was found at the phylum and genus levels in the samples obtained. With the CI approach, Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria were the most common phylum of bacteria regardless of the organ, or geographic origin of the specimens analyzed. While the relative abundance of bacteria at a phylum level with the CD approach varied between analyzed samples, the data obtained suggest that a high diversity of species of bacteria occurs in R. microplus from both livestock agroecosystems. Bacterial genera such as Anaplasma, Coxiella, and Ehrlichia, recognized for their implications in tick-borne diseases, were also detected, together with endosymbionts such as Lysinibacillus, previously reported as a potential tool for biological control. This information is useful to deepen the knowledge about microbial diversity regarding the relations between endosymbionts and pathogens and could facilitate the future development of epidemiological surveillance in livestock systems.


Ethical approval
This study was approved by the Health Research Ethics Committee of the Universidad Pontificia Bolivariana (Record No. 5 of April 24, 2017). It was also granted an environmental license issued by the Colombian Government, through the National Environmental Licensing Authority (Autoridad Nacional de Licencias Ambientales-ANLA, Resolution ANLA 744 of July 26, 2016).

Sampling sites
In the Northern region, animals grazed in predominantly Kikuyu grass surrounding pastures (Pennisetum clandestinum). The main breeds were Holstein, Jersey, and Jersey x Holstein cattle (Bos taurus). Farms were located between 2298 and 2621 meters above sea level.
In the middle Magdalena region, the main breeds were Gyr, Jersey, and Bos indicus x Bos taurus cattle grazed in the surrounding jungle/forests, with predominantly Brachiaria (Brachiaria sp.) and Jaragua grass (Hyparrhenia rufa) pastures. Farms were between 182 and 497 meters above sea level.
All farms were in the department of Antioquia, and their longitude and latitude coordinates were acquired using a global positioning system (GPS Trimble 1 Juno 1 3B, Westminster, USA). The study was conducted in five sites selected from different agroecosystems areas. The geographical distribution of the study areas is presented in Fig 1.

Sample collection and R. microplus identification
A total of 25 partially-fed female ticks (larger than 3.5 mm) were collected from 15 cattle between August and September 2017, following a protocol for direct collection of ticks from each bovine [15] and transported alive in sterile tubes. Each specimen was analyzed to confirm the growth stage, sex, and taxonomic status as R. microplus, based on morphological keys [16,17]. The microdissection of each type of organ, salivary glands or gut, was performed with a sterile scalpel blade (N˚. 11), as previously described [18,19]. Using a stereomicroscope needle (Leica EZ4, Heerbrugg, Switzerland), the dorsal cuticle flap was lifted, and each organ was removed with fine forceps and. Five female ticks from each farm were dissected, and each salivary gland, gut, or exoskeleton obtained were grouped in pools by organ (S1 Table).
Genomic DNA was extracted from all partially-fed female ticks carcass material (Exoskeleton-material) of five arthropods pooled by the municipality to obtain molecular confirmation of their taxonomic status as R. microplus. This process was completed by using the mammalian tissue genomic DNA Purification Protocol from the GeneJET™ Genomic DNA Purification Kit (Thermo Scientific™, CA, USA). During to microdissection and DNA extraction process, ticks were manipulated under sterile conditions in a Class II biological safety cabinet (BSC) using cycles of ultraviolet (UV) light between uses to avoid the contamination.
Two different molecular markers were amplified following separate PCR protocols. The fragment mitochondrial DNA 16S rRNA was amplified using 16S-F (5'-TTAAATTGCTGTRGTATT-3') and 16S-R1 (5'-CCGGTCTGAACTCASAWC-3') primers [20], whereas the fragment 18S rRNA was amplified using usingNS3 (5'-GCAAGTCTGGTGCCAGCAGCC-3') and NS4 (5'-CCGGTCTGAACTCASAWC-3') [21] primers. DNA extraction and PCR master mix preparation and amplification were performed in separate rooms to avoid contamination. The amplification protocols previously proposed [20,21] were optimized using BIOLASE™ DNA Polymerase (Bioline, London, UK) following manufacturer instructions. This optimization resulted in annealing temperatures of 50˚C and 60˚C for mitochondrial DNA 16S rRNA and 18S rRNA primers, respectively. Blank controls (nuclease-free water) were included during the PCR, and PCR products were verified through agarose gel electrophoresis. Bidirectional sequencing (Macrogen Corp., MD, USA) was performed to obtain partial sequences of each molecular marker. Reference sequences of 18S rRNA and mitochondrial DNA 16S rRNA of Ixodidae ticks were downloaded from NCBI to conduct the molecular phylogenetic analysis, which included multiple sequence alignments with ClustalW 2.1 using Geneious Prime 1 2019.0.4 (https://www.geneious.com). A molecular evolution model was selected using the Bayesian criterion to select and to conduct a molecular phylogenetic analysis by the Maximum Likelihood (ML) method using MEGAX 10.0.5 [22]. The dendrogram was estimated by a heuristic search and a Bootstrap resampling with 1000 pseudo replicas (S1 Fig).

Bacterial diversity assessment of R. microplus
Each pool of salivary glands (n = 5) or gut (n = 5) was homogenized and divided into two aliquots, one for the culture-independent (CI) approach and the other for the culture-dependent (CD) approach.

Culture-independent (CI) approach
Each pool of salivary glands or gut was immediately placed in 2 ml sterile microcentrifuge tubes containing 180 μl buffer lysis, and 20 μl Proteinase K. Genomic DNA was extracted using PureLink 1 Genomic DNA Mini Kit (Thermo Fisher Scientific, USA) following the manufacturer instructions. The DNA quality was checked by spectrophotometry in a Nanodrop 1 2000 (Thermo Fisher Scientific, USA).
PCR amplification was achieved using Bakt_341F (CCTACGGGNGGCWGCAG) and Bakt_805R (GACTACHVGGGTATCTAATCC) universal primers while targeting the V3 through V4 regions of the 16S rRNA gene [23]. PCR products were processed for library preparation with Nextera XT library preparation kit, PhiX control DNA, and sequenced with paired-end 2x250 reads on the MiSeq System (Illumina, USA) at Macrogen Inc. (Seoul, Rep. of Korea).

Culture-dependent (CD) approach
A second aliquot sample of each pool of salivary glands or gut was resuspended in vials containing 200 ul of sterile 0.9% NaCl and homogenized by vortexing to achieve organ disruption and allowing the greatest release of bacteria; then, each sample was brought to a final volume of 1000 μl and vortexing again. Subsequently, 10 μl of the suspension were directly seeded in each different aerobic culture media (blood agar, chocolate agar, Mueller Hinton agar, Hecktoen agar, and SS agar). Such volume was chosen to obtain independent colonies in the culture media. Finally, the Petri dishes were incubated at 35˚C for 48 hours.
After incubation for 48 hours, the bacteria growth as colonies on the different solid media used for the culture of each organ (intestine or salivary glands) was verified. After this, washing the surface of the solid media with a final volume of 50 ml of sterile 0.9% NaCl, ensuring all the colonies were recovered [24]. Finally, the material obtained was centrifuged at 16000 RPM, and the pellet was processed for DNA extraction, PCR amplification, and sequencing targeting the V3 through V4 regions of the 16S rRNA gene, following the same culture-independent protocol, as explained in the previous section.

Sequence processing and data analysis
Raw sequences were analyzed with Mothur v.1.39 [25], according to the MiSeq Standard Operating procedure -SOP [26]. Briefly, merged reads were filtered to a maximum length of 465 nt without ambiguous bases. Chimeras were removed with the VSEARCH algorithm [27]. OTUs with 97% of identity were built with the Opti clustering method included in Mothur v.1.39. Then, they were classified using the Ribosomal Database Project Trainset v16 MicrobiomeAnalyst [28]. Data were filtered to remove low count OTUs (parameters: minimum count 2; prevalence in samples: 10%). These were scaled using the cumulative sum scaling algorithm. Afterward, MicrobiomeAnalyst [28] web server was employed to calculate diversity indexes and core microbiome. For the former, regions were compared in terms of evenness (Shannon and Simpsons) and richness (observed OTUs). Methodological approaches, organs, and differences were tested using a T-test. For the latter, the following parameters were included: taxonomic level (genus), sample prevalence (20%), and relative abundance (0.01). Linear discriminant analysis (LDA) effective size [29] was performed to identify biomarkers through Galaxy web application [30] with the following parameters: Kruskal-Wallis test between classes: 0.05; Wilcoxon test between subclasses: 0.05 and Threshold on the logarithmic LDA score: 2. On the other hand, beta diversity was defined as the dissimilarity between communities, estimated using Principal coordinate Analysis (PcoA) ordination plot employing the Bray-Curtis dissimilarity index to evaluate diversity between each organ and methodological approaches, Permutational Multivariate Analysis of Variance (PERMANOVA) was conducted to test the association of microbial composition with covariates of interest in this study. It is important to mention that it could have possible background noise in this data, due to the absence of blank extraction controls during the DNA extraction and bioinformatics workflows.

Results and discussion
This study describes an initial exploration of the bacterial diversity of the tick species R. microplus in salivary glands and gut. Samples were collected in the Northern and Middle Magdalena, two ecologically contrasting livestock regions of Antioquia, Colombia. It also describes the bacteria identified using a direct approach with the massive sequencing of the V3 through V4 regions of the 16S rRNA gene versus the data obtained in samples subjected to a culturedependent method.

Ticks identification
The specimens analyzed here were identified as R. microplus, based on morphological characteristics confirmed by molecular and phylogenetic analyses inferred from mitochondrial DNA 16S rRNA (S1 Fig). This result confirmed the R. microplus identity of the specimens from Antioquia included in this work. Besides, it showed the close link between the specimens analyzed in this work and the clade B [31], usually restricted to R. microplus ticks found in America, Africa, Malaysia, and Japan [32,33] (S1 Fig). Regarding the phylogenetic analysis based on the 18S rRNA gene, partial sequences were grouped at the genus level as Rhipicephalus.

Bacterial identification using 16S rRNA sequencing
As S1 Table shows, there were differences between the raw reads data of the sequenced libraries obtained with the culture-independent (CI) approach and the culture-dependent (CD) approach. The total number of raw reads was 1,929,897, combining the data of both CD and CI. After applying quality filters, 636,764 reads were assigned to 9,223 raw OTUs (Operational Taxonomic Units), out of which 8,024 were unique OTUs considering 97% of identity.

Bacterial alpha diversity
Shannon's and Simpson's indices were estimated to evaluate diversity between organs, regions, and culture-independent and culture-dependent methodological approaches implemented in this study (S4 Fig). Only statistical differences were found between CI and CD approaches at the phylum level (p < 0.001) and genus level (p < 0.001, p = 0.019 respectively) (S4 Fig). Shannon's and Simpson's indices at the phylum and genus level showed a lower diversity in the CD samples than in the CI samples (S4 Fig). Considering that the data obtained at the phylum level and genus level did not show differences between the organs, it was decided that further analyses would consider the totality of OTUs identified at the genus category regardless of the tick's organ analyzed.
Regardless of the organ or geographic origin of the specimens, the four most abundant phyla detected with the CI approach were Proteobacteria, Bacteroidetes, Firmicutes, and Actinobacteria. These added to approximately 95% of the classified reads (Fig 2A), while the relative abundance of bacteria at the phylum level with the CD approach varied between the samples analyzed. In the Northern region, a dominance of Proteobacteria was detected, whereas, in the Middle Magdalena region, Firmicutes was the most abundant phylum ( Fig  2B). The differences in the bacterial diversity profiles to Northern and Middle Magdalena region were found. These variations can also depend on the host, as some authors have reported for goats and cattle [11], and lizards and mice [34].
Our findings are partially in agreement with a previous study on saliva or whole-body samples obtained from partially or fully engorged adult R. microplus females, in which Proteobacteria, Firmicutes, and Actinobacteria were also the predominant phyla; nevertheless, a relative abundance of the phylum Proteobacteria [12] higher than 75% was found.
In total, 112 bacterial genera were detected. Statistically significant differences were found in richness when CI and CD were compared (p = 0.0017). Only six genera detected with the CD approach were not seen with the CI approach. Those genera were identified as unclassified Alcaligenaceae, Alcaligenes, unclassified Bacillales, unclassified bacteria, Serratia, and unclassified Staphylococcaceae. Also, no differences were found when the bacterial genera richness between organs salivary glands or gut was compared in this study (p = 0.97) (S3 Fig).
Although the Northern region showed a higher bacterial richness when it was compared to the Middle Magdalena region, there were no statistically significant differences (p = 0.06) (S3 Fig). Xu et al. (2015) had reported that microbial diversity in R. microplus gut might be influenced by the geographic location; however, their results were based on PCR and Denaturing Gradient Gel Electrophoresis without statistical analyses [11].
In this study, unclassified Alcaligenaceae, Alcaligenes, unclassified Bacillales, unclassified bacteria, Serratia, and unclassified Staphylococcaceae were identified using the culture-dependent approach. A previous study in ticks only reported possible contamination with Several authors have reported significant variation in bacterial abundance based on the degree of tick engorgement [12,36]. Therefore, a limit to the results in this study whit partiallyfed female ticks is that it cannot be extrapolated to different nutritional status, sex or developmental stages, and the influence of such aspects to bacterial diversity associated with R. microplus should be subject to future studies.

Bacterial beta diversity
This analysis displayed statistical differences in the bacteria community between methodological approaches (CI) and (CD) (PERMANOVA: F-value = 3.4417; R2 = 0.16837; p < 0.001) and between salivary glands and gut (PERMANOVA: F-value = 2.3472; R2 = 0.12132; p = 0.004) the R2 values indicating that the covariates explained a large portion of variation between bacterial communities.

The bacterial core of R. microplus
The bacterial core community, regardless of the origin of the samples and the type of organ, showed the genera Staphylococcus and Lysinibacillus, and the family Enterobacteriaceae with the highest relative abundance in the CI approach (Fig 4A). Lactobacillus, Prevotella, Gardenella, and Coxiella were the four genera with the highest relative abundance in the CD approach (Fig 4B). On the other hand, Staphylococcus, Coxiella, and Prevotella were the three genera with the highest relative abundance in the samples from Middle Magdalena, regardless of the Comparative analysis of bacterial core community by 16S rRNA sequences obtained to samples of R. microplus using the CI and CD approaches. (A) Heatmap displays the genus-level core of the relative abundance of genera based on salivary glands and guts from R. microplus ticks from the Northern and the Middle Magdalena regions of Antioquia, data obtained using the CI approach. (B) Heatmap displays the genus-level core of the relative abundance of genera based on salivary glands and gut from R. microplus ticks from the Northern and the Middle Magdalena regions of Antioquia, data obtained using the CD approach; (%) The proportion and the color spectrum represent the abundance of each genus, dark blue represents minimum relative abundance, and dark red represents maximum relative abundance. (C) Biomarker taxa identified in R. microplus based on data comparison between the culture-independent (green) and the culture-dependent approaches (red), and non-significant (yellow).
https://doi.org/10.1371/journal.pone.0234005.g004 approach or the type of organ analyzed (Fig 5A). Stenotrophomonas and Pseudomonas were the two genera with the highest relative abundance in the Northern region (Fig 5B).
A previous study reported Coxiella and Bacillus as the most prominent bacterial genera present in seabird ticks [37]. Staphylococcus has been previously reported in R. microplus [11] and almost all tick species. This observation might be explained by the fact that it is associated with the microbiota of the skin of mammalian hosts [38]. The genus Pseudomonas is considered ubiquitous in the environment and is frequently detected as one of the main genera in several tick species such as I. scapularis [36], R. turanicus, R. sanguineus [39], A. americanum [40] and in samples of salivary glands and gut of A. maculatum [41].
The increased abundance that environmentally acquired microbes such as Bacillus, Pseudomonas, and Enterobacteriaceae is associated with decreased pathogen colonization in ticks Ixodes [42].
In this study with the CD approach, the genus Lysinibacillus (order Bacillales) was identified. This genus has been reported in ticks Dermacentor [43], in R. microplus [5] and I. scapularis [44]. Furthermore, the entomopathogenic potential of Lysinibacillus species have been used to biological control of mosquito larvae [45]. This highlights the need to determine the role of tick endosymbionts.

Biomarker analysis
The analysis of biomarkers in the CI approach identified several OTUs belonging to the phylum Bacteroidetes and Firmicutes. Regarding the former, classes Flavobacteriia and Bacteroidia were over-represented, and the latter, classes Clostridia and Negativicutes were also overrepresented (Fig 4C, OTUs highlighted in green). With the CD approach, the genus Lysinibacillus (order Bacillales) was identified as a biomarker (Fig 4C, OTUs highlighted in red). It should be noted that with the CI approach, and the genus Coxiella was detected as a biomarker. This finding is not unexpected given that some members of this genus have demanding culture conditions or reproduce in an obligate intracellular way (Fig 4C). Biomarkers found in the CD samples must be interpreted with caution since optimization and culture media may affect their reproducibility, as it has been previously reported [11].
In the Middle Magdalena region, an over-representation of OTUs of the Firmicutes phylum and the Bacilli classes were detected (Fig 5C, OTUs highlighted in red). Proteobacteria phylum and the classes Gammaproteobacteria Betaproteobacteria were predominant in the Northern region ( Fig 5C, OTUs highlighted in green). The microbial profiles might serve as biomarkers of infection frequency and transmission in ticks [46]. Therefore the bacterial diversity profiles observed in this study might provide an initial understanding of potential pathogen transmission by R. microplus in these areas.

Analysis using sequences of the V3 through V4 regions of the16S rRNA gene
In this study, the partial sequence of 16S rRNA corresponding to OTU0017 suggests the phylogenetic relationship of this OTU with different species of the genus Anaplasma; however, the definition of this relationship was not achieved at a species level. A similar result was observed with the sequence corresponding to the OTU0049 and the phylogeny analysis with reference sequences of the genus Ehrlichia (Fig 6A). A previous study that noted the relevance of Anaplasma as tick-borne bacteria, also suggested that midgut and salivary glands are specific barriers to efficient tick transmission [47]. Another study in Colombia found various species of Anaplasma and Ehrlichia circulating in different ticks species, including R. microplus from Cordoba (department to the North of Colombia) [48]. Regarding the phylogeny analysis of Coxiella sp., the OTU0019 was closely related to Coxiella-like endosymbionts (Fig 6B), previously found in R. decoloratus and R. evertsi, which have been reported as evolutionarily related but distinct to C. burnetii [49]. Because Coxiella species have been previously reported in ticks of the genera Amblyomma, Dermacentor, Ixodes, Ornitodorus and Rhipicephalus [39, 50,51], it might be useful to use a multilocus sequence typing strategy in future investigations, in order to achieve the Coxiella identification in R. microplus [52]. genus-level core of the relative abundance of genera based on salivary glands and guts from R. microplus ticks from the Middle Magdalena region, including data from both CI and CD approaches. (B) Heatmap displays the genus-level core of the relative abundance of genera based on salivary glands and guts from R. microplus ticks from the Northern region, including data with both CI and CD approaches (%). The proportion and the color spectrum represent the abundance of each genus, dark blue represents minimum relative abundance, and dark red represents maximum relative abundance. (C) Biomarker taxa identified R. microplus specimens from the Northern (green) and the Middle Magdalena (red) regions and non-significant (yellow). https://doi.org/10.1371/journal.pone.0234005.g005 In contrast to earlier findings, Rickettsia, which had been reported as one of the most frequent bacteria genera detected in R. microplus [12], Ixodes scapularis [53], Amblyomma americanum [54] and Dermacentor occidentalis, was not found in this study [55]. This is an important issue for future research, given the possibility of competition between Coxiella and Rickettsia, as it had been previously suggested in A. americanum [54]. Bacterial competition was presented for I. scapularis since it was found that the frequencies of the Rickettsiaceae family were inversely related to those of Pseudomonadaceae or Enterobacteriaceae [56]. Besides, error rates during sequencing and the effect of taxonomic resolution could interfere to some degree, hindering thus an accurate identification of bacteria genera [57]. A previous study proposed that the bioinformatic analysis may also result biased by the capability and resolution of the similarity algorithm employed [58].
Although in this study, the differences between cattle breeds were not evaluated, their influence on the microbiome of the R. microplus ticks of the two areas should be considered in future studies. Furthermore, employ other bacteria culture methodologies to permit the research of other genera that can be obtained through anaerobic bacteria culture and in the eukaryotic cell culture allowing the isolation and growth of obligate intracellular bacteria that can be transmitted by ticks such as Coxiella, Ehrlichia, Anaplasma, Rickettsia, and others.

Conclusions
This work makes an initial description of bacterial diversity in salivary glands and gut of R. microplus adult female ticks from Northern and the Middle Magdalena regions from Antioquia, Colombia. This knowledge is essential because this specie considered the most predominant tick in livestock systems. Through a complementary use of molecular procedures and traditional microbiological methods, a great diversity at the phylum and genus level was found in the R. microplus specimens analyzed. The culture-independent approach used in this study resulted in the detection of genera such as Coxiella, Anaplasma, and Ehrlichia, described as agents with public health importance. In addition to other endosymbionts, such as Lysinibacillus, previously reported as a potential tool for biological control, were detected for the first time in R. microplus. However, more studies are needed to improve this knowledge, their potential implications for tick-borne diseases, and endosymbionts related to R. microplus from South America ticks' populations.