Sequencing of the V3–V4 region of the 16S rRNA gene and data analysis

The 16S rRNA gene sequences were amplified by polymerase chain reaction using a set of primers (forward primer 5′-GCACCTAAYTGGGYDTAAAGNG-3′ and reverse primer 5′-TACNVGGGTATCTAATCC-3′). The resulting libraries were sequenced using the TruSeq Nano DNA LT library preparation kit (Illumina Inc., USA) on an Illumina Miseq platform (2 × 300 bp paired-end reads) following the manufacturer's instructions. Then, the library quality was assessed on an Agilent Bioanalyzer 2100 system (Agilent Technologies Inc., USA). The V3–V4 amplified region (200–450 bp) was demultiplexed and quality controlled using QIIME 2 [15]. In short, sequences were trimmed according to the following criteria: 1) Phred quality score ≥ 20; 2) sequence length > 200 bp; and 3) no ambiguous sequences. Next, the reads were clustered into OTUs performed using an identity threshold of 97% and the GreenGene database [16]. After clustering, the sequences were assigned taxonomy with the UCLUST clustering algorithm and following the open-reference OTUs picking protocol in QIIME 2 [17]. The OTUs (chimeric and spurious) with a relative abundance of < 0.001% of the total read counts were removed. All resulting sequences are available through the NCBI BioProject (ID PRJNA599265). For data analyses, the Chao1 index for richness, Shannon and Simpson diversity indices were estimated using the QIIME 2 software and visualized by the R software v 3.1.2 [18]. The linear discriminant analysis (LDA) effect size (LEfSe) method was used to determine the microbial differences between sexes in the collected ticks. In addition, analyses of similarities (ANOSIM) and permutational multivariate analysis of variance (PERMANOVA) were performed to compare the microbial composition among the groups using weighted and unweighted UniFrac metrics. Additionally, principal coordinate analysis (PCoA) plot was used for data visualization. The results were considered statistically significant for p value ≤ 0.05.