In May 2012, Ochotona thibetana Milne and Edwards, 1871, (Lagomorpha, Ochotonidae) was captured with mouse cages in Diqing Tibetan Autonomous Prefecture, Yunnan Province, China, and taken to the laboratory for ectoparasite collection. All specimens of F. diqingensis were collected on O. thibetana and stored in an Eppendorf tube with 95% ethanol at −80°C for subsequent identification and DNA extraction. F. diqingensis and O. thibetana were deposited at Dali University. All specimens were trapped and preserved with the approval of the Animal Ethics Committee of Dali University under the approval number MECDU-201912-20.

The morphological identification is primarily based on “Fauna Sinica insecta Siphonaptera” (14). For specimens to be photographed, we first placed them in distilled water to clean them and to stretch the wrinkled fleas. The specimens were sealed and dried, then photographed and identified in detail for morphological features using an OLYMPUS (SZ2-ILST) micrographic system. Five specimens were sent to Shanghai Winnerbio Technology Co., Ltd. (Shanghai, China). The tissues were cryoprocessed, ground, and mixed with cetyltrimethylammonium bromide (CTAB). CTAB is a cationic active agent that lysed membrane proteins and lipids, as well as depolymerized nuclear proteins. CTAB-nucleic acid complexes are soluble and stable at high salt (>0.7 mM) concentrations. Still, at low salt concentrations (0.1–0.5-mM NaCl), CTAB-nucleic acid complexes precipitate due to reduced solubility, while the majority of the proteins, polysaccharides, etc. remain dissolved in the solution. The CTAB-nucleic acid complex is separated from sugars, proteins, etc., by centrifugation, and nucleic acids can be precipitated by ethanol, while CTAB dissolved in ethanol and removed. DNA concentration and purity were measured using a spectrophotometer and electrophoretically tested for completeness. Sequencing was performed on the Illumina NovoSeq 6,000 platform (Illumina, Inc.). During Illumina sequencing, a sequencing library was constructed using about 1 μg of genomic DNA. According to the manufacturing protocol, the DNA samples were sheared into 400–500-bp fragments using a Covaris M220 Focused Acoustic Shear (Covaris, Inc.). The sheared fragments were used to prepare Illumina sequencing libraries. The prepared library was then used for paired-end Illumina sequencing (2 × 150 bp) on an Illumina NovaSeq 6,000 platform.

MitoZ 2.3 (27) was used to assemble the mitogenome, and BWA v0.7.17 (28) and Samtools v0.1.20 (29) were used to assess the confidence of the assembly data (sequencing depth ≥ 100×). Annotation was performed using Geneious Prime 11.0 (30). Protein-coding genes and rRNA genes were identified by BLAST2.16 (31) and MITOS Web Server (32). The tRNA gene was predicted using tRNAscan SE2.0 (33) and ARWEN (34), and then manually checked and corrected. The annotated mitogenome sequence of F. diqingensis has been deposited to GenBank (accession number: PP083946).

CodonW was used to analyze the codon preference of three species. Skew was calculated as AT skew = (A − T)/(A + T) and GC skew = (G − C)/(G + C). ClustalX (35) and MEGA (36) were used to calculate the conserved and variant sites of tRNA genes; DnaSP v5 (37) was used to analyze the nucleotide diversity, variable site, conserved site, and singleton site of 13 protein-coding genes.

Boreus elegans Carpenter, 1935 (38) was used as the outgroup. Bayesian inference (BI) (39) tree and Maximum likelihood (ML) (40) tree were constructed based on two concatenated data matrices (PCG: 13 protein-coding genes; PCGRNA: 13 protein-coding genes and 2 rRNA gene sequences) of 23 species (Supplementary Table S1). Sequence alignment was initially performed using MAFFT (41) software and then optimized using MACSE v 2.06 (42) software. MACSE uses the classic “Needleman–Wunsch” algorithm to correctly identify pseudogenization events and maintain the ancestral codon structure. Since comparing PCG sequences directly using MACSE was too slow, the PCG sequences were first compared using MAFFT and then optimized using MACSE. Significant gaps and ambiguous sites were omitted using Gblocks (43). The best-fit nucleotide substitution model for the tree was determined using ModelFinder (44). The best-fit nucleotide substitution model for constructing ML and BI trees is shown in Supplementary Table S2. IQ-TREE (45) was used to construct the ML tree. The number of bootstrap repetitions in the maximum likelihood analysis was set to 1,000, and the support rate of the node branch uses the maximum likelihood bootstrap proportion (BSP). When BSP ≥ 70%, the node has a high support rate. MrBayes 3.2.6 (39) was used to construct the BI tree. The BI tree underwent 3,000,000 generations, with trees sampled every 1,000, each with four independent Markov Chain Monte Carlo (MCMC). When Bayesian posterior probabilities (BPP) ≥ 95%, the posterior probability had a higher support rate. Figtree (46) was used to visualize and edit the evolutionary tree.

Based on the reconstructed phylogenetic tree, the divergence time of 23 species of 15 genera in 7 flea families was estimated using Bayesian Evolutionary Analysis Sampling Trees (BEAST 2.0) (47), with data from 13 PCGs. Due to limited flea fossil records and available sequence data, fossil calibrations (amber specimens, see Introduction section) were employed to confirm true Siphonaptera, which convey verifiable knowledge of the existence of specific lineages of Siphonaptera over time. (1) The genus Pulex Linnaeus, 1758. Although this flea fossil was placed into the Miocene, the age range of Dominican amber was controversial, with estimates running from 20 to 15 Mya (Iturralde–Vinent and MacPhee, 1996) to 45–30 Mya (Schlee, 1990). It was set to 35–15 Mya to reflect this uncertainty (1, 48). (2) Ctenophthalmidae: 40–35 Mya was used for calibration points, and the root node was not calibrated. First, the best model for accurate phylogenetic estimation was selected using ModeFinder (44). The most suitable replacement model, GTR, was set in BEAUti, along with the adoption of a strict molecular clock model. The tree prior model used the Yule process, and the prior distribution was set to normal to estimate the divergence time. The Markov chain (MCMC chain) ran for 40,000,000 generations, with trees sampled every 1,000 generations. The convergence of the parameter estimates was performed using Tracer. The ESS value of the parameters was >200. The top 10% of the tree was discarded as burn-in. Trees were edited and visualized using tvBOT (49).

From a careful examination of specimens of F. diqingensis, the following characteristics were used for identification: The sternum IX (rear arm) is not developed and narrow. It bears six bristles. There is a row of seven bristles above it, and the spacing between the two bristles is very short. The top of the immovable process is up to the midpoint of the movable process, which is higher than F. spadix Jordan and Rothschild, 1921. The immovable process is conical, and its upper anterior angle is mostly a right angle. The concave degree of the upper part of the trailing edge of the movable process is usually deeper. There are three sub-thorns in the middle and lower part of the trailing edge of the movable process. There are three spiniforms in the middle and lower part of the trailing edge of the movable protrusion. The shape of slermum VII (rear arm) differs from F. spadix. The tail boundary of spermatheca is unclear (Figure 1).

The coding region of F. diqingensis mitogenome (GenBank accession number: PP083946) is 16,153 bp in size (the non-coding region is not complete, which will not be discussed in this article) (Figure 2). The nucleotide composition and skew are shown in Table 1. Among the 37 genes, 4 PCGs (nad1, nad4, nad4L, and nad5), 8 tRNA genes (trnY, trnC, trnQ, trnV, trnL₁ (tag), trnP, trnH, and trnF), and 2 rRNA genes (rrnL and rrnS) were encoded on the N-chain, while the remaining 23 genes (nad2, cox1, cox2, atp8, atp6, cox3, nad3, nad6, cob, trnI, trnM, trnW, trnL₂ (taa), trnK, trnD, trnG, trnA, trnR, trnN, trnS₁ (tct), trnE, trnS₂ (tga), and trnT) were encoded on the J-chain. The AT-skew of the F. diqingensis mitogenomes was −0.027, and the GC-skew was −0.228. There are a total of 13 intergenic spacer regions (106 bp), with the largest intergenic spacer region being 32 bp, followed by 21 bp. There are 11 overlap regions (28 bp), with the largest overlap region being 7 bp and the smallest overlap region being 1 bp.

The total length of the 13 protein-coding genes was 11,121 bp, the longest was nad5 (1,716 bp), and the shortest was atp8 (162 bp). All 13 protein-coding genes use the typical ATN as the start codon (ATG: 7, ATA: 3, ATT: 2, and ATC: 1), and the stop codons are TAA, TAG, and incomplete T (Table 1). The 13 protein-coding genes of F. diqingensis encode a total of 3,707 amino acids. The most frequently used code was UUA (RSCU = 4.46), and the lowest frequency used code was CCG (RSCU = 0.1) (Supplementary Table S4). The average length of 22 tRNAs was 65.3 ± 2.59 bp, and the size of a single gene was 62–70 bp. The secondary structure of tRNA is shown in Supplementary Figure S1. Among the 22 tRNA genes, except for trnS₁ (tct) lacking D-arm, all other tRNA genes formed typical cloverleaf structures. In addition to the typical Waston-Crick pairing (A-T, G-C), a total of 23 mismatches (G-T: 9, T-G: 7, T–T: 5, and A-A: 2) occurred in the mitogenome folding process of F. diqingensis, of which trnH mismatch was the most occurred (4 times) (Supplementary Figure S1). The total length of the two rRNAs was 2093 bp, both located on the N chain. The AT content of rrnL was 81.8%, and the AT content of rrnS was 80.1%.

F. diqingensis, F. spadix, and P. custodis showed high AT (75.5, 78.8, and 76.8%) contents. The AT-skew of F. diqingensis was −0.027, and the GC-skew was −0.228. The AT-skew of F. spadix was −0.036 and the GC-skew was −0.214. The AT-skew of P. custodis was −0.008, and the GC-skew was −0.259. The AT content of the first, second, and third codons of protein-coding genes of three species of flea was analyzed, and it was found that the AT content in the first and second codons was much smaller than that in the third codon. It shows that the third codon of the three species in Amphipsyllinae prefers to use A and T bases (Table 2). The protein-coding genes of three species use the typical ATN as the start codon, with ATG being used most frequently (Supplementary Table S5). Except for nad2 and nad4 with incomplete T or TA as the stop codon, the remaining 11 protein-coding use the typical TAN as the stop codon, with TAA used most frequently (Supplementary Table S5).

The tRNA secondary structures of F. diqingensis, F. spadix, and P. custodis in Amphipsyllinae have little difference among different species, and all of them form typical cloverleaf structure (except trnS₁ (tct)). The 22 tRNAs had 1,458 sites, of which 1,312 were conserved, accounting for 90.0% of the total sites. There were 146 variable sites, accounting for 10.0% of the total sites. There were 14 mutation sites of trnC and 1 mutation site of trnQ (Supplementary Table S6). The majority of these mutation sites were at the junction of arms and rings (Supplementary Figure S2).

The nucleotide diversity index (π) of 13 protein-coding genes in Amphipsyllinae was 0.3842. There were 3,991 variable sites (V), 7,147 conserved sites (C), and 3,890 singleton sites (S).

The relative synonymous codon usage of three species of Amphipsyllinae was analyzed and found that the codon usage preference patterns of F. diqingensis, F. spadix and P. custodis were very similar (Supplementary Table S4, Supplementary Figure S3). F. diqingensis, F. spadix, and P. custodis had 26, 27, and 27 high-frequency codons (RSCU >1), respectively. It can be seen from RSCU that there is a specific preference for codon usage, but the reason for its preference is not apparent. We use parity rule 2 (PR2), effective number of codon (ENC), neutrality curve, and corresponding analysis (COA) to analyze the reasons that affect codon usage preference. PR2 showed that mutation, selection, and other factors affected codon usage through qualitative analysis and found that protein-coding genes of three species of fleas were not located at the central site (0.5, 0.5). Among them, F. diqingensis and F. spadix have a T/C preference, and P. custodis has an A/C preference (Figure 3). The ENC values of protein-coding genes of 3 species in Amphipsyllinae ranged from 30.11 to 46.54. The majority of them were below the standard curve (the atp8 gene is not included because CodonW calculates the NC value based on the principle of ENC (50) calculation, and when the gene does not contain amino acids with eight synonymous genes, the ENC value of this sequence is not calculated). Among them, 13 protein-coding genes Nc ≤ 35 have a total of 15 genes, F. diqingensis has 4 genes, F. spadix has 9 genes, and P. custodis has 2 genes (Supplementary Figure S4). In the neutrality curve, GC12 of three species has no significant linear relationship with GC3s (the coefficient of determination R² is used mainly to quantify the strength of the linear relationship between two variables. The closer the R² value is to 1, the stronger the linear relationship between GC12 and GC3s; the closer the R² value is to 0, the weaker the linear relationship between GC12 and GC3s). Therefore, the absolute value of the regression coefficient cannot represent the degree of mutation pressure on protein-coding genes in the species (Figure 4). In the COA analysis, three species of fleas are scattered together. There is no situation in which one species is clustered alone, indicating that F. diqingensis, F. spadix, and P. custodis have similar codon usage preferences (Figure 5).

We obtained four trees from two concatenated data matrices (PCG and PCGRNA) and two tree-construction methods (ML and BI) (Figure 6). Four trees are highly similar; the difference is that the taxonomic status of Hystrichopsylla weida qinlingensis Zhang et al. (1984), Dorcadia ioffi Smit (1953), and Aviostivalius aklossi bispiniformis Li and Wang (1958), is ambiguous. In Figure 6A, PCG: The 23 flea species are divided into two major lineages, forming three clades: (Pulicidae + Pygiopsyllidae) forms one clade, while ((Ceratophyllidae + Leptopsyllidae) + ((Vermipsyllidae + Hystrichopsyllidae) + Ctenophthalmidae)) forms two clades that are more closely related. In Figure 6B, PCGRNA: The 23 flea species are divided into two major lineages, forming three clades: Pygiopsyllidae form one clade, while ((Ceratophyllidae + Leptopsyllidae) + ((Vermipsyllidae + Hystrichopsyllidae) + Ctenophthalmidae)), and Pulicidae forms two clades that are more closely related. The families Leptopsyllidae, Ctenophthalmidae, and the subfamily Ceratophyllinae are paraphyletic groups. The families Pulicidae and Ceratophyllidae form a monophyletic group. The families Pulicidae and Pygiopsyllidae are closely related to the base of the phylogenetic tree. In Amphipsyllinae, F. diqingensis and F. spadix were clustered into one clade (BPP = 1, BSP = 100). P. custodis and Macrostylophora euteles Jordan et Rothschild, 1,911 were clustered into one clade (BPP = 1, BSP = 100). The families Ctenophthalmidae, Vermipsyllidae (one species), and Hystrichopsyllidae (one species) had low node support.

According to the estimated divergence time, Siphonaptera originated in the Early Cretaceous around 119.3 Mya (the divergence times and confidence intervals for nodes/clades are shown in Supplementary Table S7). Clade 1 (Pygiopsyllidae and Pulicidae) begins to diverge from clades 2 (Ctenophthalmidae; Vermipsyllidae; and Hystrichopsyllidae) and 3 (Leptopsyllidae and Ceratophyllidae) at 101.21 Mya. Clade 2 (Ctenophthalmidae; Vermipsyllidae; and Hystrichopsyllidae) and Clade 3 (Leptopsyllidae and Ceratophyllidae) began to diverge at 86.26 Mya. The divergence time between Pygiopsyllidae and Pulicidae is approximately 95.36 Mya, and the relationship between Pygiopsyllidae and Pulicidae is closer. Amphipsyllinae diverged approximately between 73.37 and 40.32 Mya. The Ceratophyllidae family diverged most recently at 73.37 Mya, while F. spadix and F. diqingensis diverged at approximately 40.32 Mya. The majority of the intraordinal divergence into Siphonaptera’s extant lineages occurred after the K-Pg boundary (Figure 7).