All ticks were collected from Wuwei City (102.65°E, 37.94°N), Gansu province, in May 2022. Feeding ticks were carefully extracted from sheep using tweezers while questing ticks were gathered from grassland using the Flag Cloth Law method, which involves dragging a 90 cm × 60 cm white flannelette flag across the surface. It should be noted that the sampling of both questing and feeding ticks was opportunistic and was therefore not standardized with respect to collection area or time.

Following collection, all ticks were immediately stored in tubes with perforated lids and refrigerated at 4°C to minimize tick mortality. Afterward, ticks’ surfaces were cleaned using 75% ethyl alcohol, and identification was carried out with a stereoscopic microscope (RH-2000, HIROX, Japan), based on observed structural characteristics of the ticks.

Three ticks, morphologically identified as Dermacentor at the genus level, were placed in a tube and labeled according to the collecting date and lifestyles for subsequent operations. Each tube was injected with 1 mL of Stroke-Physiological Saline Solution (SPSS, 0.9% NaCl), and 3–4 0.2 mm steel balls were added before the samples were ground at a frequency of 65 Hz for 500 s. Following this, the samples were incubated in a dry water bath for 10 min at a consistent temperature of 55°C, and the supernatant fluid was carefully transferred to a new tube after a rapid centrifugation at 12,000 rpm for 1 min in a standard centrifuge. Nucleic acid extraction was then performed using the Ex-DNA/RNA nucleic acid extraction kit (TIAN LONG, China) according to the manufacturer’s instructions (Shen et al., 2020), and the extracted samples were stored at −20°C for further analysis.

Tick identification was achieved by amplifying DNA using PCR, specifically targeting the mitochondrial gene cytochrome c oxidase I (COI; Lv et al., 2014). The primers used for the molecular biological identification of ticks were COI (Forward primer: 5′-GGAACAATATATTTAATTTTTGG-3′, Reverse primer: 5′-ATCTATCCCTACTGTAAATATATG-3′).

The PCR reactions contained 40 ng of genomic DNA and were conducted in 50 μL reaction volumes including 25 μL 2× DreamTaq PCR Master Mix (Thermo, America), 2 μL of each primer (10 μmol/L), 2 μL of DNA (20 ng/μL), and 19 μL of ddH₂O. Reactions were performed in a thermocycler under the following conditions: an initial preheating step at 95°C for 3 min, followed by 35 cycles of denaturation at 95°C for 30 s; annealing at 45°C for 30 s; extension at 72°C for 1 min; and a final extension at 72°C for 10 min.

The PCR products were subjected to 1.2% agarose gel electrophoresis (AGE) under specific conditions, including a voltage of 220 V, a current of 400 mA, and a duration of 30 min. Following this, the gel was imaged using a gel imager to visually represent the DNA fragments. For precise identification of the tick species, the nucleotide sequence of the PCR products was determined using Sanger sequencing, which was facilitated by AuGCT DNA-SYN Biotechnology Co., Ltd. located in Beijing, China.

Aimed to amplify the sequences of species-specific for the 16S rRNA gene of the ribosome’s large subunit, the methods followed involved the generation of a DNA library using QIAseq FX DNA Library Kit (QIAGEN, Germany), as per the manufacturer’s instructions, before diluting the nucleic acid for subsequent sequencing, so that the concentration could be detected by the Qubit-4 (Thermo, America).

The primary protocol encompassed several steps, including library preparation via the fusion method, amplicon purification, quantification, and pooling. The final products were individually validated using AGE and any samples with concentrations falling below the standard 50 ng/μL, as defined by the Nextera XT Index Kit v2 (Illumina), were deemed unqualified and hence excluded from the pool. All pooled DNA samples were then subjected to paired-end sequencing, utilizing the MiSeq Reagent Kit V3 on the Illumina HiSeq platform (Odendaal et al., 2022). This sequencing process adhered strictly to the manufacturer’s instructions and was characterized by an insert size of 350 bp and a read length of 250 bp.

The barcode-based data were analyzed using the Quantitative Insights Into Microbial Ecology version2 (QIIME2) software suite (Caporaso et al., 2010), which provides a comprehensive software environment, data standards, and tool wrappers. The sequence data, presented in FASTQ format, were demultiplexed and consolidated into the same directory before input as QIIME artifacts (.qza) or QIIME visualizations (.qzv). After importing the paired-end reads from the original DNA fragments into the virtual machine, we utilized the DADA2 method to denoise data by truncating both forward and reverse sequences at 230 base pairs each. Following the DADA2 protocol, the sequences were sorted and grouped into Autosave (ASV), which is replaced by Operational Taxonomic Unit (OTU) in the following. The sequences were matched against the Silva database according to the Uchime algorithm (Edgar et al., 2011). Upon creation or acquisition of the classifier artifact, the “qiime taxa barplot” command was used to perform taxonomic analysis of the OTU representative sequences. The bacterial composition of each sample was visualized in a bar figure using the R packages ggplot2 and RColorBrewer (Wickham, 2009; Neuwirth, 2022).

For single-sample diversity analysis or alpha diversity, we determined sufficient sequencing depth based on rarefaction curves for the observed number of OTUs across all samples shown as illustrated in the Good’s coverage. After the sequences of all samples were randomly drawn to a uniform data volume, we calculated the number of unique OTUs, the community richness as determined by the Shannon estimator, and the community evenness as indicated by the Pielou index in each sample using the R package Vegan (Dixon, 2003). The Wilcoxon rank-sum test was employed to assess whether these indices differed significantly and the analysis of similarities (ANOSIM) was used to quantify the variation in bacterial composition resulting from different tick lifestyles using package Vegan.

Following the exclusion of Rickettsiaceae, the abundance of bacteria in samples was recalculated and beta diversity analyses were performed at the family level using the R packages Vegan, phyloseq, and ggplot2 (McMurdie and Holmes, 2013). Weighted and unweighted UniFrac analyses were utilized to examine diversity and to compare different groups, which were subsequently plotted in a principal coordinate analysis (PCoA). The Welch’s t-test was applied to assess differences in bacterial composition between two tick lifestyles. Variations in bacterial composition between groups were visualized using STAMP software (Parks et al., 2014). Lastly, to identify microbial taxa with significant differences between Groups A and B from phylum to species level, we employed Linear Discriminant Analysis (LDA) Effect Size (LEfSe; Segata et al., 2011). A threshold LDA score of >2 was used within 95% confidence intervals.¹

The representative COI gene sequences and pathogen reads obtained were cross-referenced with existing results in GenBank, employing the Basic Local Alignment Search Tool (BLAST) search engine provided by the National Center for Biotechnology Information (NCBI).² Multiple sequence alignment was conducted using the ClustalW algorithm within the MEGA-11 software suite, applying default parameters. The phylogenetic tree was then constructed using the Kimura two-parameter model of Neighbor-Joining method based on MEGA-11, with bootstrap values estimated from 1,000 replicates.

The Spearman rank correlation test was conducted to present the interrelationship between two bacteria in samples using the R packages Hmisc and Igraph (Csardi and Nepusz, 2005), and the networks were constructed to visualize relationships between OTUs at the genus level. The two bacteria were linked by red or blue lines which represented a positive or negative correlativity, respectively, a process facilitated by Gephi software.³ OTUs that were either unassigned or occurred only once were disregarded at the genus level and a cut-off value was established for reads to eliminate pathogens. Microbial occurrences that represented less than 10% of the frequency distribution of reads in a microorganism were discarded, and the correlation coefficient (|r|-value) was set at >0.6 with a 95% confidence level (p < 0.05).

The ticks were categorized into two groups based on their mode of capture: Group A consisted of feeding ticks, while Group B comprised questing ticks. A total of 150 adult ticks were identified as Dermacentor nuttalli (Acari: Ixodidae) according to statistics, confirmed through morphological identifications and a PCR assay of the species-specific COI region, which a high identity range of 92.66–99.75%. The phylogeny is shown in Figure 1.

Prior to the 16S metagenomics sequencing, our laboratory conducted nested PCR assays using specific primers, through which five pathogens, namely Rickettsia, Anaplasma, Borrelia burgdorferi, Babesia, and Bartonella, were detected. The primer sequences are outlined in Table 1, displayed from 5′ to 3′ for both forward and reverse sequences.

Of all the tested nucleic acid samples, the prevalence rates for Rickettsia and Anaplasma stood at 32.67 and 4.67%, respectively, surpassing those of the other three pathogens. Notably, Borrelia burgdorferi, Babesia, and Bartonella were not detected using conventional PCR assays. A subset of the samples (4.67%) demonstrated detected with Rickettsia and Anaplasma simultaneously. Results of the positive rate are presented in Supplementary Table 1.

Using the Illumina HiSeq platform and after eliminating low-quality sequences, 5,480,349 reads were obtained in the data analysis. The bacterial genera were determined by comparing the taxonomic profiles at the genus level against the Silva database. Approximately 8 × 10⁶ paired-end V3-V4 16S reads were procured from all samples, with the read count varying from 7,332 to 417,497 and Supplementary Table 2 included 16S NGS sample productions. After pre-processed (merged, quality filtered, and removal of singletons and chimeras) and post-processed, a cumulative total of 5 million sequences were obtained for the samples and their replicates. The average number of sequences was 109,607 and a range from 42,877 to 257,450 as illustrated in Figure 2A. The rarefaction curves, which level off well before 10,000 reads (Figure 2B), indicate that the number of reads is sufficient to compile a reliable list of bacterial genera. The raw data had uploaded to NCBI with the Bio-Project Accession Number PRJNA1015185.

At the family level, Rickettsiaceae, Coxiellaceae, and Enterobacteriaceae represented an average of 84.21, 9.42, and 3.84%, respectively. The top 20 taxonomics for each sample are illustrated in Figure 3A. Although Rickettsiaceae is predominant in most samples, the Coxiellaceae and the Fancisellaceae are significantly prevalent in samples number 25 (33,450 reads, 68.18%) and 50 (65,231 reads, 48.29%) respectively. As shown in Figure 3B, Rickettsiaceae comprises 2,330,700 (81.11%) and 2,241,292 (87.71%) of the reads in Groups A and B, respectively. Notably, Coxiellaceae and Enterobacteriaceae in Group A (339,077 reads, 11.80%; 160,804 reads, 5.60%) have approximately double and triple the reads in Group B (172,405 reads, 6.75%, 47,821 reads, 1.87%) respectively. In contrast, the Francisellaceae family, although less prominent, shows a slight distinction between the groups, with Group B (66,723 reads, 2.61%) having more reads and proportions than Group A (1,187 reads, 0.04%).

And at the genus level, Rickettsia is the dominant genus across all samples, comprising between 25.36 and 98.95% of the total community (average: 82.63%), followed by Coxiella and Enterobacter, which account for averages of 9.42, and 3.49%, respectively. Figure 3C provides a visual representation of the relative abundance of the top 25 species at the genus level. The figure illustrates that the genus-level classifications of Rickettsia and Ac37b collectively constitute the family-level classification of Rickettsiaceae. The comparison of the top 25 species between groups is shown in Figure 3D.

Regarding alpha diversity, there were no significant differences between groups for either Shannon Index (W = 359, p = 0.375) or Simpson Index (W = 353, p = 0.441), despite the range of the two Index being wider for Group A compared to Group B, as presented in Supplementary Figure 1. The Wilcoxon rank test of other indexes also had no statistical differences which are shown in Supplementary Table 3. There were no significant differences in bacterial microbiome compositions between parasitic and questing ticks according to ANOSIM (pseudo-R = 0.019, p = 0.175), as illustrated previously in Figure 2F.

However, after removing Rickettsiaceae computationally, the composition unveils noteworthy variations between groups (Figure 2C). It is apparent that Figure 2D explicitly presents the dominant role of Coxiellaceae within two groups (A: 339,077 reads, 62.46%; B: 172,405 reads, 54.88%), and subsequently, a significant disparity in the presence of Enterobacteriaceae (A: 160,804 reads, 29.62%; B: 47,821 reads, 15.22%) is discernible in Figure 4. Interestingly, while the relative abundance of Francisellaceae in Group B (21.24%) surpasses that in Group A (0.22%), no statistical difference in microbiome composition at the family level is noted (p = 0.55, Wilcoxon rank-sum test) as depicted in Figure 2E. A comparative analysis of the Ixodes microbiome diversity reveals greater species evenness in Group B (0.628) relative to Group A (0.520) based on the Simpson index, while the Shannon index suggests that the species richness in Group A (1.044) is lower than in Group B (1.363) as shown in Supplementary Table 4.

Principal coordinate analysis identifies the species basically similar between the groups. The Unweighted UniFrac PCoA further reveals that PC1 (26.44%) surpasses PC2 (9.34%), and in the Weighted UniFrac PCoA plot, PC2 (36.13%) is inferior to PC1 (46.65%), as represented in Figures 5A,B. A comparison of the relative abundance of the top 20 bacterial genera using the Welch’s t-test (two sides), identifies no significant differences at the family level within the 95% confidence intervals, as shown in Figure 6A. While Enterobacteriaceae is marginally more prevalent in Group A than Group B, the mean proportion of Francisellaceae and Coxiellaceae is comparatively diminished. However, substantial differences become apparent at the genus level as shown in Figure 6B, with the result of the relative abundance of Friedmanniella and Bordetella being elevated in parasitic ticks as compared to questing ticks. Figure 7A demonstrates differences among specimens at p < 0.05 as determined by LEfSe analysis, and a slight discrepancy is detectable within 95% confidence intervals from phylum to species in Figure 7B, with 11 taxa in Group B and 10 taxa in Group A detected at higher relative abundance.

For the selected representative reads from Rickettsia, Anaplasma, Francisella, and Coxiella, the phylogenetic tree, as shown in Figure 8, manifests sequences alignment implemented with MEGA 11 software, utilizing the Kimura two-parameter model of Neighbor-Joining method. Table 2 presents the average intraspecies pairwise genetic distances for the four Rickettsia and 14 reference sequences retrieved from NCBI, which range from 0.045 to 0.173 and 0.000 to 0.012, respectively. The table also indicates the pairwise divergence (Divergence = 1 − percent identity) between stains, denoted above the diagonal in Table 2.

The network, as depicted in Figure 9A, illustrates the correlational relationship among bacterial genera based on co-occurrence patterns, with varying line colors representing different relationships between two genera. The network consists of 36 nodes and 133 links, with only two instances of negative correlation noted between Francisella and Phyllobacterium (r = −0.65), and between Francisella and Arthrobacter (r = −0.63). There are 22 pairs exhibiting strong positive correlation (r > 0.8), with the top three pairs being Sphingobacteriales-uncultured and Phyllobacterium (r = 0.88), Blastocatellia-uncultured and Propionibacteriaceae-uncultured (r = 0.92), and Enterobacter and Escherichia-Shigella (r = 0.98). The detailed results for all bacteria are provided in Supplementary Table 5. In order to examine the correlation between bacterial genera across groups, we plotted network diagrams for Group A (Figure 9B) and Group B (Figure 9C). And there are 390 correlational relationships among bacteria in Group A and more than 350 in Group B, which are shown in Supplementary Tables 6, 7.