Identification of tick-borne pathogens using metagenomic analyses in H. longicornis feeding on humans in downtown Beijing

Abstract

On August 14th, 2018, a Beijing resident living in Xicheng District found a female H. longicornis tick attached to the skin at the front of his upper shin. On examination, the patient was afebrile and appeared well. The species of the tick was identified through morphological characteristics and phylogenetic analysis based on cytochrome C oxidase subunit I. This H. longicornis tick was screened for tick-borne pathogens such as viruses, bacteria and parasites. RNA pathogens were screened by PCR and sequencing, while DNA pathogens were screened by metagenomic analyses. It was found that the tick was positive for the DNA sequences of zoonotic and animal pathogens such as A. phagocytophilum, Ehrlichia minasensis and C. burnetii. Considering the good health condition of the patient, we hypothesized that the pathogens originated from the tick specimen itself rather than host blood meal. For the first time, our study reveals the possible risk of transmission of tick-borne pathogens to human beings through tick bit in downtown Beijing. Further research is needed to screen for tick-borne pathogens among unfed ticks collected from central Beijing.

Background

Ticks were the first arthropods to be recognized as vectors that can transmit pathogens to human beings and are second only to mosquitoes as vectors of infectious diseases (Dantas-Torres et al. 2012). More than 30 emerging tick-borne pathogens have been confirmed in mainland China since the 1980s, such as species of the spotted fever group rickettsiae and species in the family Anaplasmataceae (Fang et al. 2015). Humans are occasional host of ticks, and these emerging tick-borne agents have become a real threat to human health in China (Socolovschi et al. 2009).The northwestern part of Beijing, especially Yanqing and Huairou Districts, is dominated by the Jundu Mountains, while the western part, Mentougou and Fangshan Districts, is framed by Xishan or the Western Hills (Cheng et al. 2016). Forest areas and grasslands in the suburban mountains of Beijing provide habitats for livestock, wildlife and ticks (Li et al. 2002; Lu et al. 2013). There are several reports that ticks bit humans in the suburban areas of Beijing (Li et al. 2002; Lu et al. 2013). Until now, the cases of humans bitten by ticks reported in Beijing have all occurred in suburban areas rather than downtown areas. Downtown Beijing includes Chaoyang, Haidian, Fengtai and Shijingshan Districts, especially the Dongcheng and Xicheng Districts, which are in the center of Beijing.

In this study, a case of a person bitten by H. longicornis that occurred in Xicheng District is reported for the first time. Generally, ticks inhabit suburban areas, and ticks are rarely found in downtown Beijing. Tick-borne pathogens, including viruses, bacteria and parasites, were screened in this tick specimen, and our results showed that this tick was positive for the DNA sequences of zoonotic and animal pathogens such as A. phagocytophilum, Ehrlichia minasensis and C. burnetii. This highlights the risk of transmission of tick-borne pathogens from ticks to humans through the bit in the downtown area of Beijing, as H. longicornis is a tick vector of high medical significance.

Results

Case report

A 65-year-old healthy man discovered a tick attached to the skin on the front of his upper shin on August 14th, 2018. He had not traveled outside Xicheng District of Beijing within 7 d, and he had no exposure to farm animals, wild animals or pets. Within the previous week, the patient had not been outside his resident community or surrounding markets/streets. After identified it as a possible tick through visual inspection, this tick was removed from the skin by wiping it with alcohol. The area around the attachment site was raised, solid and soft, with slightly red concentric rash around it. There were no other symptoms, such as fever, chills, headache, fatigue or muscle aches. As the patient’s symptoms were mild, serological examination was not recommended. After a week of observation, the patient did not develop worsening symptoms, and the previous symptoms, such as concentric rash, vanished without medical attention.

Species identification

Morphological characteristics of body dorsal and ventral views, dorsal integument, capitulum and legs of female ticks are used for species identification. The tick was identified as a female Hemaphysalis tick by its unique combination of morphological characters.

Sequence analysis revealed that COI sequence derived from the specimen shared 99% sequence identity with COI for H. longicornis (MG721044.1 and JQ737096.1). In addition, phylogenetic analysis using sequences derived from COI sequences confirmed that the species of the tick specimen was H. longicornis (Fig. 1).

Fig. 1

Phylogenetic analyses of tick species based on COI. Neighbor-joining phylogenetic analysis based on partial tick COI sequences. Bootstrap values are indicated at the nodes. The scale bar indicates the degree of divergence represented by a given length of a branch. The red dot indicates the COI sequence acquired in this study

Molecular detection of tick-associated RNA pathogens through PCR

RNA viruses belonging to Flavivirus, Phlebovirus, Nairovirus and Hantavirus can cause serious human diseases (Honig et al. 2004; Lv et al. 2018; Matsuno et al. 2015; Wójcik-Fatla et al. 2011). Viruses belonging to the first three have been demonstrated to spread to people through ticks, while Hantavirus can be spread through rodents (Honig et al. 2004; Lv et al. 2018; Matsuno et al. 2015; Wójcik-Fatla et al. 2011). In this study, RNA was extracted from the H. longicornis specimen and was screened for the presence of Flavivirus, Phlebovirus, Nairobivirus and Hantavirus. Results showed that it was negative for these viruses.

Investigation of tick-associated DNA pathogens via metagenomic analyses

Metagenomic sequencing and analyses were utilized to screen the DNA of pathogens in the tick specimen. Krona analyses showed that several bacterial pathogens belonging to Alphaproteobacteria and Gammaproteobacteria were detected in the H. longicornis specimen. A. phagocytophilum, Ehrlichia minasensis and Rickettsia amblymmatis belonging to the Alphaproteobacteria detected in the tick specimen are shown in Fig. 2. The counts of mappable reads are indicators of the abundance of a species, although the two are not exactly the same. A. phagocytophilum accounted for 51% of mappable reads for the bacterial species within Alphaproteobacteria, while E. minasensis accounted for 24%. C. burnetii and Coxiella-like endosymbionts belonging to Gammaproteobacteria detected in the tick specimen are shown in Fig. 3. C. burnetii accounted for 16% of the counts of mappable reads for the bacterial species within Gammaproteobacteria. The relative abundance of DNA segment reads of bacterial pathogens is shown in Fig. 4. A. phagocytophilum, E. minasensis and C. burnetti were among the most abundant DNA pathogens within the tick specimen. DNA sequences of tick-borne parasites such as Babesia spp. and Theileria spp. were not detected in this tick specimen.

Fig. 2

Hierarchical classification and distribution of various species belonging to Alphaproteobacteria detected in the tick specimen using KRONA extension. The red dots indicate pathogens that can infect humans

Fig. 3

Hierarchical classification and distribution of various species belonging to Gammaproteobacteria detected in tick specimens using KRONA extension. The red dots indicate pathogens that can infect humans

Fig. 4

The relative abundance of various bacterial species detected in the tick specimen. The red dots indicate pathogens that can infect humans

Discussion

H. longicornis ticks, commonly associated with livestock (e.g., sheep, goats, horses and cattle), and wild animals such as deer and birds, occasionally bite humans (Choe et al. 2011). A previous investigation showed that the dominant tick species in forest areas of Beijing was H. longicornis, followed by H. concinna, D. silvarum, Ixodes persulcatus and Rhipicephalus sanguineus (Li et al. 2002). Several cases in which humans were bitten by ticks have been reported in forested areas in Beijing, such as Huairou and Yanqing districts, and these patients all had a travel history to the suburbs of Beijing (Liu et al. 2020). Fortunately, none of the patients in the above cases experienced serious symptoms and did not require any further treatment. In fact, the number of reports of tick bites is very limited. Although there is the possibility that some cases were not reported, the frequency of tick bites is probably not very high in Beijing. In this study, the patient’s recent 7 days of outdoor exposure before the tick bite was restricted to his resident community and the surrounding markets/streets. H. longicornis tick was found attached to the patient’s skin, and there was little blood within the tick. It is speculated that the length of blood-feeding time of the tick on the patient was not long and that the tric bit patient in the lawns of his resident community or adjacent areas in Xicheng District. This research showed a risk of being bitten by H. longicornis tick in the downtown area of Beijing for the first time.

Severe fever with thrombocytopenia syndrome bunyavirus (SFTSV) is a currently emerging pathogen, and H. longicornis has been proven to be the definite reservoir and vector responsible for SFTSV transmission to humans (Zhuang et al. 2018). In this study, bacterial pathogens, including A. phagocytophilum, E. minasensis, R. amblymmatis, C. burnetii and Coxiella-like endosymbionts, were detected in H. longicornis. A. phagocytophilum and C. burnetii were among the most abundant zoonotic pathogens within the tick specimen (Kim et al. 2003; Lee et al. 2004). Fortunately, the patient’s mild symptoms vanished soon without medical care and he did not develop a clinical illness over the following 3 months. Several key factors affect the infection by tick-borne pathogens through tick bites: the activity and quantity of pathogens, the infection ability of pathogens, the immunity capability of patients, and the length of blood-feeding time of patients. It is suspected that a short length of blood-feeding time and a high level of immunity of healthy hosts prevent the infection of these pathogens. However, there have been several reports about tick-borne pathogens such as A. phagocytophilum in humans with tick bites (Pascoe et al. 2019; Jahfari et al. 2016). The probability of infection with a tick-borne pathogen other than Lyme spirochetes after a tick bite is approximately 2.4% in the Netherlands (Jahfari et al. 2016). This study revealed the possible risk of transmission of tick-borne pathogens to human beings through tick bites in downtown Beijing.

Previous studies have shown that Ehrlichia spp. such as E. chaffeensis have been detected in H. longicornis ticks (Lee and Chae 2010). In this study, E. minasensis was detected from H. longicornis specimens in Beijing. To our knowledge, this is the first report of E. minasensis in China (Cabezas-Cruz et al. 2016; Cicculli et al. 2019). Based on the results of this study, it is speculated that the prevalence of E. minasensis in the local livestock of Beijing is possible, but further investigation is needed.

• R. amblyommii is a pathogen belonging to the Rickettsiae spotted fever group (SFG) (Merhej et al. 2014; Blanton LS et al. 2014). R. amblyommii was detected in A. americanum ticks in the United States of America, and then R. amblyommii was detected in Am. pseudoconcolor ticks in Brazil (Goddard and Norment 1986; Silva et al. 2018). To our knowledge, R. amblyommii has been detected in H. longicornis in China for the first time in this case. The presence of R. amblyommii increases the complexity status of spotted fever rickettsiosis in China as R. raoultii, R. sibirica, Candidatus R. longicornii and R. jingxinensis belonging to SFG have been previously proven to be present in China (Liu et al. 2020; Yu et al. 1993).

Conclusion

We certified that this tick was positive for the DNA sequences of zoonotic and animal pathogens such as A. phagocytophilum, Ehrlichia minasensis and C. burnetii by using metagenomic analyses. The tick, identified as H. longicornis, fed on a human in the downtown area of Beijing.

Materials and methods

Specimen identification

The tick was kept alive and sent to the Institute of Animal Inspection and Quarantine (IAIQ), Chinese Academy of Inspection and Quarantine (CAIQ), for species identification. It was killed in hot water (80 °C) according to a previously described method and dried on filter paper (Soares et al. 2013). Digital images of the specimen were taken with a stereomicroscope (Discovery V20, Zeiss, Oberkochen, Germany). The tick was identified to species by morphology according to standard morphological characteristics, including body dorsal and ventral views, dorsal integument, capitulum and legs, and corroborated by deriving sequences for mitochondrial cytochrome C oxidase subunit I (COI) (Murrell et al. 2001). Digital images of the tick specimen are provided in Fig 5. Partial tissue of the tick was used for ribonucleic acid (RNA) and deoxyribonucleic acid (DNA) extraction. The remaining tissue of the tick was preserved at − 80 °C.

Fig. 5

Morphological characteristics of the tick (females) by stereomicroscopy. A dorsal view (20×); B ventral view (20×)

DNA and RNA extractions from the tick

The tick was washed three times with distilled water and dried on filter paper. Then it was dissected with disposable scalpels: one quarter utilized for DNA extraction and another quarter utilized for RNA extraction. The remaining half was preserved at − 80 °C. DNA and RNA were extracted according to a previously described method (Lv et al. 2018).

Molecular identification of the tick species

Sequences of COI were amplified by PCR using the primer pair LCO1490/HCO2198 (Folmer et al. 1994). This universal primer pair was designed to amplify a 710 bp fragment of COI from 11 invertebrate phyla, including ticks. PCR amplification and sequence alignment were analyzed as previously described (Lv et al. 2014). COI sequence has been submitted to GenBank, and the accession number is MZ452024.

Detection of tick-associated RNA viruses through PCR

Reverse transcription was carried out with M-MLV reverse transcriptase (Promega, WI, USA) in a reaction volume of 40 μL, which included 18 μL extracted RNA, 8 μL 5× RT buffer, 2 μL dNTPs (10 mM), 3 μL DTT (0.1 M), 4 μL M-MLV reverse transcriptase (200 U/μL), 1 μL RNasin (40 U/μL), 2 μL 10× hexanucleotide mix, and 2 μL molecular grade H₂O.

Detection of Flavivirus RNA was conducted using hemi-nested PCR targeting the RNA-dependent RNA polymerase gene, as described previously (Lv et al. 2018). In brief, 5 μL cDNA was utilized as a template for each reaction. Reaction master mix was prepared according to the manufacturer’s protocol, and PCR conditions described in previous studies were used (Lv et al. 2018).

Detection of Hantavirus RNA was conducted using nested PCR targeting the L-segment, as described previously (Wójcik-Fatla et al. 2011). In brief, 5 μL cDNA was utilized as a template for each reaction. Reaction master mix was prepared according to the manufacturer’s protocol, and PCR conditions described in previous studies were used.

Detection of Phlebovirus (L-segment) was performed by PCR using a mixture of ppL1/ppL2 primers (40 pmol each) with the reaction conditions previously described (Matsuno et al. 2015).

Detection of Nairovirus (L-segment) was performed by PCR using 6942+ and 7385- primers with the reaction conditions previously described (Honig et al. 2004).

Sequence and phylogenetic analyses of PCR-amplified DNA segments

DNA amplicons of the correct size were sequenced by BGI Tech Solutions Co., LTD (Liuhe, Beijing). DNA sequence was assembled and edited in MEGA X (Kumar et al. 2018). Sequence alignments were conducted as described (Lv et al. 2018). Sequence alignments were conducted using ClustalW within MEGA V.7.0 (Kumar et al. 2018) using default parameters (open gap penalty =10.0, extend gap penalty = 5.0) before subsequently being checked by visual inspection. Bootstrapping (1000 replicates) was utilized to estimate node support. Pairwise deletion was used for gaps/missing data. Based on the K2P distances, phylogenetic trees were constructed with the combined datasets of all major tick genera using the neighbor-joining method. For COI analysis, all codon positions and noncodon sites were combined and tested.

Metagenomic sequencing and analyses

Both quality and quantity of the extracted DNA were assessed by measuring the absorbance at wavelengths of 260 nm and 280 nm. One microgram of high-quality DNA was used for subsequent library construction. A sequencing library was generated using the NEBNext® Ultra™ DNA Library Prep Kit from Illumina (NEB, USA) based on the kit instructions. The Illumina HiSeq platform was used to sequence the library preparations. The specific processing steps were as follows: a) remove the reads that contained low-quality bases (default quality threshold value ≤38) above a certain portion (default length of 40 bp); b) remove the reads in which the N base reached a certain percentage (default length of 10 bp); c) remove reads that shared the overlap above a certain portion with Adapter (default length of 15 bp). Considering the possibility of host pollution in sample, the clean data needed to be blasted to the host database, which defaults using Bowtie 2.2.4 software (Langmead and Salzberg 2012) to filter the reads that are of host origin. The parameters are as follows according to the previous description: --end-to-end, −-sensitive, −I 200, −X 400 (Karlsson et al. 2012).

The samples combined and then assembled using SOAP denovo (Luo et al. 2012)/MEGAHIT (Li et al. 2015) software as previously described (Brum et al. 2015; Nielsen et al. 2014; Qin et al. 2014).

CD-HIT software (Fu et al. 2012) was utilized to cluster sequences (parameters -c 0.95, −G 0, −aS 0.9, −g 1, −d 0), and MetaGeneMark (Qin et al. 2012) was used to predict and analyze the ORF (open reading frame). Clean data of the sample were mapped to the initial gene catalog using Bowtie 2.2.4, and the number of reads to which genes mapped in the sample with parameter settings as --end-to-end, −-sensitive, −I 200, and −X 400. Genes with ≤2 reads in the sample were filtered, and gene catalog (Unigenes) was obtained and eventually used for subsequent analysis. Based on the number of mapped reads and the length of genes, the abundance information of each gene in the sample was statistically analyzed. Abundance of a species in one sample equals the sum of gene abundance annotated for the species; gene number of a species in a sample equals the number of genes whose abundance is nonzero (Li et al. 2014; Fu et al. 2012; Qin et al. 2012; Le Chatelier et al. 2013). The basic information statistics, core-pan gene analysis, correlation analysis of the sample and Venn figure analysis of the number of genes were all based on the abundance of each gene in the gene catalog.

Taxonomy prediction

DIAMOND software (V.0.9.9, https://github.com/bbuchfink/diamond/) was used to blast unigenes to the sequences of bacteria, fungi, archaea and viruses, which were all extracted from the NR database (Version: 2018–01–02, https://www.ncbi.nlm.nih.gov/) of NCBI with the parameter settings blastp-e 1e-5 (Buchfink et al. 2015).

To ensure the species annotation information of the sequences, we chose the result for which the e value ≤ the smallest e value * 10 to take the LCA algorithm, which was applied to the system classification of the MEGAN software (Huson et al. 2011).

Data containing number of genes and abundance information of the sample in each taxonomy hierarchy (kingdom, phylum, class, order, family, genus and species) were obtained based on the LCA annotation result and the gene abundance result.

Krona analysis, the exhibition of generation situation of relative abundance, the exhibition of abundance cluster heat map, PCA (R ade4 package, V.2.15.3) (Avershina et al. 2013) and NMDS (R vegan package, V.2.15.3) (Magali Noval Rivas et al. 2013) decrease-dimension analysis were based on the abundance table of each taxonomic hierarchy. Differences between groups were analyzed according to the references (White et al. 2009). LEfSe analysis was conducted by LEfSe software (the default LDA score was 3) (Segata et al. 2011). Finally, random forest (RandoForest) (Williams 2014) was used to construct a random forest model. Important species were screened out by MeanDecreaseAccuracy and MeanDecreaseGin, cross-validated with each model (default 10 times) and plotted with the ROC curve.