Population genetics and genetic variation of Ctenocephalides felis and Pulex irritans in China by analysis of nuclear and mitochondrial genes

Background Fleas are the most economically significant blood-feeding ectoparasites worldwide. Ctenocephalides felis and Pulex irritans can parasitize various animals closely related to humans and are of high veterinary significance. Methods In this study, 82 samples were collected from 7 provinces of China. Through studying the nuclear genes ITS1 and EF-1α and two different mitochondrial genes cox1 and cox2, the population genetics and genetic variation of C. felis and P. irritans in China were further investigated. Results The intraspecies differences between C. felis and P. irritans ranged from 0 to 3.9%. The interspecific variance in the EF-1α, cox1, and cox2 sequences was 8.2–18.3%, while the ITS1 sequence was 50.1–52.2%. High genetic diversity was observed in both C. felis and P. irritans, and the nucleotide diversity of cox1 was higher than that of cox2. Moderate gene flow was detected in the C. felis and P. irritans populations. Both species possessed many haplotypes, but the haplotype distribution was uneven. Fu's Fs and Tajima's D tests showed that C. felis and P. irritans experienced a bottleneck effect in Guangxi Zhuang Autonomous Region and Henan province. Evolutionary analysis suggested that C. felis may have two geographical lineages in China, while no multiple lineages of P.irritans were found. Conclusions Using sequence comparison and the construction of phylogenetic trees, we found a moderate amount of gene flow in the C. felis and P. irritans populations. Both species possessed many haplotypes, but the distribution of haplotypes varied among the provinces. Fu’s Fs and Tajima’s D tests indicated that both species had experienced a bottleneck effect in Guangxi and Henan provinces. Evolutionary analysis suggested that C. felis may have two geographical lineages in China, while no multiple lineages of P.irritans were found. This study will help better understand fleas' population genetics and evolutionary biology. Graphical Abstract

Ctenocephalides felis and the human flea Pulex irritans are characterized by various hosts, including humans, carnivores, rodents, and ungulates [8,9] and exhibit high genetic diversity and complex genetic structure. C. felis and P. irritans are closely related to human and animal life and are, therefore, of high veterinary significance [8,10]. Understanding fleas' vector transmission and disease epidemiology will help develop strategies to prevent and control them [11,12].
With the development of genetic technologies, which complement, to some extent, the limitations of traditional morphology, molecular biology methods have been widely used in taxonomy, systematics, and population genetics, including those of fleas [13,14]. The genetic diversity of fleas has also been studied using nuclear markers, such as histone H3, EF-1α, ITS1, and ITS2 [13,[15][16][17]. Mitochondria caught the attention of evolutionary biologists in the 1960s because of their small size, high abundance in cells, and simple mode of inheritance [18,19]. The rapid mutation rates of mtDNA compared with nDNA facilitates the analysis of recent divergences within and between species [20]. The mtDNA markers cytochrome c oxidase subunits 1 and 2, named cox1 and cox2, are frequently used to evaluate genetic diversity and identify cryptic species and the population structure of invertebrates [10,[21][22][23]. Van der Mescht et al. suggested in 2015 that host specificity may influence the level of intraspecies genetic differentiation [24]. Hornok et al. investigated this hypothesis in 2018 and found large differences in mitochondrial sequences in some synanthropic flea species, such as C. felis and P. irritans [25].
In 2019, Zurita et al. studied the classification, origin, evolution, and phylogeny of P. irritans from two populations in Spain and Argentina, using the internal transcribed spacer (ITS)1 and ITS2 and partial cytochrome c oxidase subunit 1 (cox1) and cytochrome b (cytb) genes. The study found a considerable degree of molecular differentiation between the two populations and found two clear geographic genetic lineages, indicating the presence of two cryptic species [17].
Fleas have previously been studied in parts of China [29,32], but on the Chinese mainland, there have been no studies into the intraspecific genetic diversity, evolutionary relationships, or transmission patterns of C. felis and P. irritans. In this study, two nDNA markers, ITS1 and EF-1α, were combined with two different mtDNA markers, cox1 and cox2, to (i) provide detailed information that expands the knowledge of the genetic diversity of C. felis and P. irritans in various provinces of China; (ii) provide a basis for the re-evaluation of the population structure and genetic differentiation of C. felis and P. irritans in China; and (iii) infer phylogenetic relationships from the geographical distribution of different isolates.

Parasites and DNA extraction
A total of 82 fleas were obtained from different regions of China (Fig. 1). The fleas were observed and photographed with a stereomicroscope, and the features of the shape of the frons, mane, comb, and other parts were used for preliminary identification [10,13,25,[33][34][35]. All fleas were supplied in ≥ 70% (v/v) ethanol and stored at − 20 °C. The whole specimen was used to extract DNA. According to the specification; total genomic DNA was isolated using sodium dodecyl sulfate/proteinase K treatment, followed by spin-column purification (Wizard ® SV Genomic DNA Purification System; Promega Madison, WI, USA) from the manufacturer. The sample codes, host, locality, and GenBank accession numbers are listed in Table 1.
All PCR reactions (25 μl) contained 1 µl of DNA sample, 0.5 μl of each primer, 12.5 μl of 2 × Premix Emer-aldAmp Max HS PCR Master Mix (TaKaRa, Dalian, China), and 10.5 μl of ddH 2 O. All PCRs were run on a thermocycler (Bio-Rad, Hercules, CA, USA). The results were validated by electrophoresis on 1% agarose gels. PCR products were bi-directionally sequenced by Sangon Biotech Co., Ltd. (Shanghai, China). BLAST analysis was performed on the sequencing results of four genes of all samples to further determine the species of flea collected.

Population differentiation
All ITS1, EF-1α, cox1, and cox2 sequences were compared using Clustal X 1.83 software and then manually cut. MegAlign software was used to analyze the sequence differences among four genes [38]. The DnaSP 5.0 program was used to determine the haplotypes, nucleotide diversity (Pi), and haplotype diversity (Hd) of each gene [39]. DnaSP was also used to calculate the correlation among geographic locations, genetic differentiation index (Fst), and the corresponding gene flow (Nm) between populations [12]. For the cox1, cox2, and concatenated sequences, DnaSP 5.0 was used to calculate Tajima's D [40] and Fu's Fs [41] statistics to study the population history of C. felis in six regions and P. irritans in two regions. DnaSP 5.0 and PopART 1.7 were used to create a statistically parsimonious network to infer haplotype relationships [42].

Sequences analysis and reconstruction of phylogenetic relationships
The EF-1α sequence and the tandem sequence of cox1 + cox2 of C. felis and P. irritans from China were compared and analyzed. For the EF-1α sequence, Megabothris calcarifer was used as the outgroup, and for cox1 + cox2, Ceratophyllus wui was used as the outgroup.
Clustal X 1.83 was used to align the sequences of these four genes separately [38], and then the EF-1α sequences and the concatenated sequences of cox1 and cox2 were compared using Gblock 0.91 [43]. Using MEGA X, the phylogenetic trees were constructed using the maximum likelihood method [44]. Relative support for the tree topology was obtained by bootstrapping using 1000 replicates.
For C. felis, the Hd and Pi were 0.834 and 0.18944, respectively, according to cox1 sequence analysis. The Hd was 0.819, and the Pi was 0.02544, according to cox2 sequencing analysis. In tandem sequence analysis, the Hd was 0.901, and the Pi was 0.09145 (Table 2). In the overall mean data, the cox1 gene showed higher Hd than the cox2 gene, while for the Pi, cox1 showed higher genetic diversity than cox2. For P. irritans, Hd and Pi were 0.917 and 0.07117, respectively, according to cox1 sequence analysis. In cox2 sequence analysis, Hd and Pi were 0.889 and 0.45925, respectively. The results of tandem sequence analysis indicated that the Hd and Pi were 0.9763 and 0.29144, respectively ( Table 2). The Hd of the cox1 gene was higher than that of the cox2 gene, while the genetic diversity of Pi of the cox1 gene was lower than that of the cox2 gene.
According to the known classification criteria, Fst > 0.25, is defined as great differentiation; 0.15 < Fst < 0.25 is defined as moderate differentiation, and Fst < 0.15 is defined as a trivial difference [45]. According to Slatkin [46], Nm > 1 indicates high gene flow; 0.25 < Nm < 0.99 is intermediate gene flow; and Nm < 0.25 is classified as low gene flow. For C. felis, Fst, and Nm among the six regions are shown in Table 3

Population expansion
In the existing C. felis, the Fu's F neutrality test of cox1, cox2, and the serial sequence were positive, except for the cox2 sequence in Guangxi Zhuang Autonomous Region, and cox1, cox2, and the serial sequence in Jiangsu province. In Tajima's D neutrality test, cox1 from Jiangsu province, cox2 from Guangdong province, and cox1, cox2, and their serial sequences from Hunan province were all positive, while cox2 from Guangxi Zhuang Autonomous Region had no Tajima's D neutrality test value, and the rest were negative. There were no statistically significant differences except for cox1 in Henan province and cox1 and the serial sequence in Guangxi Zhuang Autonomous Region. The significance of the neutrality test in Guangxi Zhuang Autonomous Region is a sign of population balance. cox2 in Henan province and cox1 and the tandem sequence in Guangxi Zhuang Autonomous Region were positive, indicating that these two regions' populations shrank dramatically because of equilibrium selection and bottleneck effects ( Table 2).
For P. irritans, Fu's F neutrality test was positive except for cox1, cox2, and the serial sequence from Sichuan province. In the neutrality test of Tajima's D, cox2 from Sichuan province was positive, and the rest were negative. Only the cox1 test values of these two regions had statistical significance, among which the population neutral test value of Henan province was significant, a sign of population balance.

Network analysis
The complete network of haplotypes of C. felis and P. irritans was constructed (Figs. 2 and 3). In C. felis, 17 cox1 haplotypes (A1-17) were distributed in a star shape around the A1 haplotype, which contained 21 individuals, including 4 from Hunan province, 5 from Guangdong province, 5 from Jiangsu province, and 7 from Shanxi province. Henan province had the largest haplotypes, with ten haplotypes (A2, A5, A6, A7, A8, A9, A10, A11, A12, and A13), while Guangxi Zhuang Autonomous Region had the fewest haplotypes (A3 and A4). The A4 haplotype was found only in Guangxi Zhuang Autonomous Region, the A6-A13 haplotypes were found only in Henan province, the A14 and A15 haplotypes were found only in Hunan province, the A16 haplotype was found only in Jiangsu province, and the A17 haplotype was found only in Shanxi province (Fig. 2, Table 4). Compared with cox1, cox2 haplotypes (B1-B8) were fewer and less widely distributed (Fig. 2). Guangdong province had the most haplotypes (B1-B4), while Guangxi Zhuang Autonomous Region had only one haplotype, B5. The B2 haplotype was found only in Guangdong province, the B6 haplotype was found only in Henan province, and the B8 haplotype was found only in Jiangsu province. Among the 21 haplotypes generated by the sequence data of the cox1 + cox2 series, AB11 was at the center of the remaining 20 haplotypes, forming a star pattern (Fig. 2). There were ten haplotypes (AB2, AB10, AB14, AB15, AB16, AB17, AB18, AB19, AB20, and AB21) in Henan province and two haplotypes (AB5 and AB64) in Guangxi Zhuang Autonomous Region, results which were the same as those of the cox1 haplotype. Haplotypes AB3 and AB4 were found only in Guangdong province, the AB5 and AB6 haplotypes were found only in Guangxi Zhuang Autonomous Region, the AB14-21 haplotypes were found only in Henan province, the AB7 and AB8 haplotypes were found only in Hunan province, the AB9, AB11, and AB12 haplotypes were found only in Jiangsu province, and the AB13 haplotype was found only in Shanxi province. In the populations of C. felis, the most Table 4 Haplotype distribution of C. felis and P.irritans from China based on cox1, cox2, and concatenated cox1 + cox2 genes corresponding to the seven sampling locations. The number in parentheses refers to the relative frequency of the given haplotype common haplotypes were A1 of cox1, B1 of cox2, and AB1 of cox1 + cox2 (Fig. 2).
In the cox1 network of P. irritans, 16 haplotypes (C2-17) formed a star shape with the C1 haplotype as the center. C1 was the only haplotype that appeared in both Sichuan and Henan provinces. The C2-12 haplotypes were found only in Henan province, and the C13-17 haplotypes were found only in Sichuan province (Fig. 3, Table 4). Among the 13 haplotypes of cox2 (D1-13), the D2 and D3 haplotypes were found only in Sichuan province, and the D4-13 haplotypes were found only in Henan province, except D1 was found in both places (Fig. 3, Table 4). Twenty haplotypes (CD1-20) were generated from the cox1 + cox2 series data. The CD1-7 haplotypes were found only in Sichuan province, and the CD8-20 haplotypes were found only in Henan province (Fig. 3, Table 4). In the populations of P. irritans, the most common haplotypes were C1 of cox1, D1 of cox2, and CD1 of cox1 + cox2 (Fig. 3).

Phylogenetic analysis
To study the phylogenetic relationships among C. felis and P. irritans isolates from different regions in China, we analyzed sequences of EF-1α and the tandem sequences The phylogenetic tree constructed based on the EF-1α is shown in Fig. 4, and the phylogenetic tree of the tandem sequences of the cox1 and cox2 genes is shown in Fig. 5. The C. felis and P. irritans samples in these two trees were separated with high support. In the evolutionary tree constructed based on EF-1α, C. felis was not divided into two distinct branches. However, in the evolutionary tree constructed from the cox1 and cox2 series, C. felis could be divided into two distinct clades; almost all the samples from Guangxi Zhuang Autonomous Region clustered into one branch, and the rest clustered into another. In all of the evolutionary trees, P. irritans clustered into a single branch.

Discussion
Flea is one of the important blood-sucking ectoparasites in veterinary medicine, which can also seriously affect human health. With the development of molecular Phylogenetic relationships between C.felis and P.irritans based on cox1 + cox2 sequences genetics, the understanding of flea genetic diversity and population dynamics has become the premise for elucidating flea's basic biology and population characteristics, which is of great significance for the detection and control of fleas. This study showed high genetic diversity based on cox1 + cox2 sequences of C. felis and P. irritan from China. The overall genetic diversity of C. felis from China was higher than that of Malaysia [12], Thailand, and Fiji and Seychelles [10]. The nucleotide diversity of cox1 was higher than that of cox2 in C. felis and P. irritans. This finding is consistent with those of Lawrence et al. [10], who found higher levels of nucleotide diversity in the cox1 sequence in the C. felis population of Australia and contradicts those of Azrizal-Wahid et al. [12]. The neutral test results of Fu's Fs and Tajima's D showed that C. felis and P. irritans in Guangxi Zhuang Autonomous Region and Henan province experienced equilibrium selection and bottleneck effect. Due to long-term interactions between the parasite and the host, certain aspects of the parasite tend to closely track the biological characteristics of the host [47]. Therefore, it stands to reason that the population history of the host may influence the population history of the flea.
C. felis and P. irritan showed several common haplotypes (A1; A4; B1; B4; AB1; AB5; C1; D1; CD1) and many rare haplotypes (Figs. 2 and 3, Table 4). Multiple cox1 gene haplotypes were detected in C. felis from Australia, Hong Kong, New Zealand, and Malaysia [12,[29][30][31], and only one cox2 gene haplotype was found in C. felis from Australia [30]. However, in our study, the number of cox1 and cox2 gene haplotypes in C. felis was higher than those in the above study. The high genetic diversity of cox1 and cox2 genes in C. felis and P. irritan in China may be responsible for many haplotypes. In this study, the most common haplotypes of C. felis were A1 haplotype of cox1 sequence, B1 haplotype of cox2 sequence, and AB1 haplotype of cox1 + cox2 sequence, while the most common haplotypes of P. irritan were C1 haplotype of cox1 sequence, D1 haplotype of cox2 sequence, and CD1 haplotype of cox1 + cox2 sequence. This finding suggests that A1, B1, AB1, C1, D1, and CD1 may be the oldest haplotypes because of their higher frequency and wider distribution [48]. At the same time, the Fst values of C. felis and P. irritan were higher (Fst = 0.34897, Fst = 0.37433), indicating that the genetic differentiation of C. felis and P. irritan was greater. Ctenocephalides felis and P. irritan had moderate gene flow (Nm = 0.93, Nm = 0.84), which may result from restricted gene flow between different regions. China is a vast country with complex terrain. Flea transmission between regions may be limited despite well-developed road transportation and frequent human exchanges, resulting in a low distribution of flea haplotypes and limited gene flow.
To further understand the population outcome and genetic differentiation of fleas, C. felis and P. irritan in this study were compared with previous results. Among them, C. felis from Vietnam and Philippines matched their cox1 gene with haplotype A4 from Guangxi Zhuang Autonomous Region in this study [32]. The haplotype H1 produced by the cox2 gene from Malaysian C. felis studied by Azrizal-Wahid et al. [12] was consistent with the haplotype B5 from Guangxi Zhuang Autonomous Region in this study. For P. irritans, haplotypes C1 and C10 matched haplotypes H6 and H4 from Spain, respectively [17]. Other haplotypes in this study were not recorded in other countries, indicating that the mutations in our study may not have been dispersed or mutations in other areas have not entered our study area. From the perspective of geographical relationship, southeast Asia is relatively close to Guangxi Zhuang Autonomous Region of China, which may be the reason for the haplotype matching of C. felis.
Phylogenetic analysis based on cox1 + cox2 sequences showed that C. felis could be divided into two different branches, and C. felis from Guangxi Zhuang Autonomous Region had undergone genetic differentiation. Our study identified two clades of C. felis from China through the cox1 and cox2 genes of mtDNA, indicating two geographical lineages. Only one branch in the evolutionary tree was constructed by EF-1α and cox1 + cox2 sequences. There were two geographical lineages of P.irritans in Spain and Argentina in previous studies, with two possible hidden species [17]. In our study, no multiple lineages of P.irritans were found, which may be due to the lack of hidden species of P.irritans in China or the limitation of sample size.
Our study showed that C. felis populations were geographically isolated, and high levels of genetic diversity and moderate levels of gene flow indicated that C. felis and P.irritans had undergone genetic differentiation and lineages recombination. C. felis may have two genetic lineages, and the Guangxi Zhuang Autonomous Region has undergone genetic differentiation. In contrast, no multiple lineages have been found in human fleas. These results complement previous studies and better understand Chinese fleas population genetics and evolutionary biology.

Conclusions
In this study, 59 C. felis samples and 23 P. irritans samples from seven provinces in China were analyzed for population genetic variation and structure using the nuclear genes ITS1 and EF-1α and the mitochondrial genes cox1 and cox2. The results showed that C. felis populations