Adult specimens were sampled during 2012–2022 from Yunnan, Chongqing, Guangxi, Jilin and Hubei Province (Table S1). The specimens were preserved in 95% alcohol and at –20°C before the DNA extraction. Genomic DNA was extracted using the DNeasy Blood and Tissue kit (QIAGEN, Hilden, Germany) from the thoracic muscle tissue. The DNA concentration was measured using a nucleic acid protein analyzer (Thermo Scientific, Waltham, MA, USA).

Whole genomes were sequenced on the Illumina HiSeq 2500 and Illumina X plus. Raw reads were checked by FastQC 0.11.9 (Brown et al. 2017) and lo w-quality reads were filtered using Trimmomatic 0.32 (Bolger et al. 2014) with an average quality threshold of Q15. The mitogenomes were assembled by IDBA-UD 1.1.3 (Peng et al. 2012). Annotations were conducted using MitoZ 2.4 (Meng et al. 2019) with the “invertebrate mitochondrial code 5” as genetic code and “Arthropoda” as clade. Then, the sequence was checked by manual proofreading according to related species. The circular maps of mitogenomes were drawn using OGDRAW (Greiner et al. 2019), and the base composition and codon usage were analyzed using MEGA 7.0 (Kumar et al. 2016). The calculation formulas of base composition asymmetry were AT-skew = (A−T)/(A+T) and GC-skew = (G−C)/(G+C) (Perna and Kocher 1995). The relative synonymous codon usage (RSCU) of PCGs were calculated using MEGA 7.0 (Kumar et al. 2016). Sequences of 13 PCGs of all Culex (including downloaded from Genbank and the seven newly sequenced mitogenomes) were individually aligned via the L-INS-i algorithm using MAFFT 7.313 (Katoh and Standley 2013) with the “invertebrate mitochondrial” as code table, “codon” as alignment mode and “auto” as strategy. Then, nucleotide diversity (Pi) and non-synonymous/synonymous substitution ratios (Ka/Ks) were calculated using DnaSP 6 (Rozas et al. 2017) with the “first site” as protein coding regions and the “mtDNA Drosophila” as genetic code. The secondary structures of tRNAs were predicted using the MITOS Web Server (http://mitos.bioinf.uni-leipzig.de/index.py) and checked through manual proofreading (Bernt et al. 2013), and the secondary structures of rrnS and rrnL were predicted using RNA Structure (http://rna.urmc.rochester.edu/RNAstructureWeb).

There were 50 specimens of Culicidae in the ingroup and two specimens of Simuliidae were used as outgroup representatives to reconstruct the phylogeny in this study (Table S2). Sequences of PCGs and rRNAs were aligned using MAFFT 7.313 (Katoh and Standley 2013). Each rRNA gene alignment was conducted using the G-INS-i algorithm by MAFFT 7.313 (Katoh and Standley 2013) with the “invertebrate mitochondrial” as code table, “normal” as alignment mode and “auto” as strategy. The ambiguous regions in the alignments of PCGs and rRNAs were removed using Gblocks (Talavera and Castresana 2007). Maximum Likelihood (ML) trees were inferred using IQ-TREE v.1.6.10 (Nguyen et al. 2015) with 1000 ultrafast bootstraps. Bayesian Inference (BI) analyses were conducted using MrBayes v.3.2.2 (Ronquist and Huelsenbeck 2003) under the best substitution models selected by ModelFinder 2.1.10 (Kalyaanamoorthy et al. 2017). Finally, the phylogenetic trees were visualized using FigTree 1.4.4 (Rambaut 2009).

Divergence times in Culicidae were estimated using MCMCTree v4.9e (Yang 2007). The topology obtained through BI analysis of the PCGAA dataset was incorporated as the input file. Three fossil calibrations were selected: Toxorhynchites mexicanus (Zavortink and Poinar, 2000) (16 million years ago Mya) as the earliest fossil record of Toxorhynchites; Cx. winchesteri Cockerell, 1908 (34 Mya) as the earliest fossil of Culex; Priscoculex burmanicus, Zavortink & Brown, 2009 (89 Mya) as the earliest fossil of Anophelinae (https://mosquito-taxonomic inventory.myspecies.info). The independent rate clock model (clock = 2) was used for further calculation. Two independent MCMC chains were conducted for 20 million generations with 25% discarded as burn-in. The convergence of the MCMC chains was evaluated using Tracer v1.7.2 (Rambaut et al. 2018), ensuring that the ESS values exceeded 200.

The complete mitogenomes of the seven newly sequenced species (Fig. S1) range from 15,564 to 15,971 bp in length, exhibiting the typical gene content for the genus Culex: 13 PCGs, 22 tRNAs, two rRNAs, and a control region (CR). Gene strandedness were conserved, with 23 genes (nine PCGs, 14 tRNAs) encoded on the heavy strand (J-strand) and 14 genes (four PCGs, eight tRNAs, and two rRNAs) on the light strand (N-strand) (Tables S3–S9).

There are 34 complete mitogenomes of Culex including those published on GenBank (https://www.ncbi.nlm.nih.gov/genbank), ranging from 15,564 bp (Cx. bitaeniorhynchus) to 16,052 bp (Cx. chidesteri) in length. Across these mitogenomes, the content of C is higher than G and the content of A is higher than T, with A + T content ranging from 77.98% to 79.57% (Fig. 1, Table S10). The length of complete PCGs ranges from 11,199 bp (Cx. bitaeniorhynchus) to 11,232 bp (Cx. huangae). AT content in PCGs ranges from 78.17% (Cx. lygrus) to 76.11% (Cx. fergusoni) (Fig. 1, Table S11).

For complete mitogenomes, AT-skew values range from 0.00 to 0.01, and GC-skew values range from –0.19 to –0.15. In PCGs specifically, AT-skew ranges from –0.17 to –0.15, and GC-skew ranges from 0.03 to 0.07 (Fig. 1, Tables S10, S11), a pattern consistent with previous reports in mosquitoes (Krzywinski et al. 2011; Martinez-Villegas et al. 2019). Notably, base composition analysis at the three codon positions reveals that variation is most pronounced at the third position (Table S12). Due to codon degeneracy, this nucleotide variation doesn’t alter the encoded amino acid sequences (Table S13; Reetz et al. 2008).

This pattern indicates that purifying selection acts primarily at the protein level, while synonymous substitutions accumulate at the nucleotide level – a feature that makes mitochondrial PCGs informative for phylogenetic studies across different taxonomic scales (Zardoya and Meyer 1996; Simon et al. 2006).

Start and stop codon usage across the 34 Culex mitogenomes is largely conserved, with some gene- and species-specific exceptions (Tables S14, S15). While most PCGs initiate with the canonical ATN start codon, two notable deviations are observed: the ND1 gene initiates with TTG in 11 species, and the ATP6 gene in Cx. vishnui also utilizes TTG. The recurrent use of the non-canonical TTG start codon for ND1 across multiple lineages suggests either a shared ancestral feature or a conserved mechanism of translational initiation within the genus (similar non-canonical starts have been reported in other insect mitogenomes) (Cameron 2014). Furthermore, COX1 consistently employs TCG as the start codon across all 34 species. Stop codon usage also varies. Most PCGs (ATP6, ATP8, ND1, ND2, ND4L, ND5, ND6) use the complete TAA stop codons in all species examined, with the exceptions of ND2 in Cx. harrisoni (TAG) and ND6 in Cx. vishnui (TA-tRNA). Incomplete stop codons (T-- or TA-), completed via post-transcriptional polyadenylation, are observed in COX1, COX2, COX3, CYTB, ND3, and ND4, with the specific type varying by species.

Relative synonymous codon usages (RSCUs) in the seven newly sequenced mitogenomes are presented in Figure S2 (Tables S16–S22). The most frequently used codons are UUU (Phe), UAU (Tyr), AUU (Ile), AUA (Met), AAU (Asn) and AAA (Lys). Notably, CAA (Gln) is not used in these seven mitogenomes. This codon bias reflects the high A + T content of the genomes, favoring codons rich in A and T.

Nucleotide diversity (Pi) and the rate of nucleotide substitution (Ka/Ks) are calculated for the 13 PCGs (Fig. 2). Pi values range from 0.041 (ND4L) to 0.080 (CYTB). The Ka/Ks ratios for all 13 PCGs are less than 1, indicating that these genes evolve under purifying selection (Hurst 2002; Hurst 2009; Mori and Matsunami 2018). Among them, COX1 exhibits the lowest Ka/Ks ratio (0.170), suggesting the strongest purifying selection and the slowest evolutionary rate, while ATP8 shows the highest ratio (0.170), indicating a relatively faster evolutionary rate.

The strong purifying selection on COX1 reinforces its suitability as a reliable marker for deep phylogenetic relationships (Folmer et al. 1994; Hebert et al. 2003), while the elevated evolutionary rate of ATP8 may provide greater resolution for separating closely related species (Gissi et al. 2010).

The CR, located between tRNA^Val and tRNA^Ile in all seven newly sequenced mitogenomes, is the most variable region of the Culex mitogenome. Across the 34 species, CR length ranges from 580 to 1195 bp (Table S23), with A + T content varying from 87.89% to 91.56%. AT-skew ranges from –0.13 to 0.02, and GC-skew from –0.46 to –0.24. Two long tandem repeats (>200 bp) are identified in the CR of Cx. bicornutus, and their predicted secondary structures are shown in Figure 3. Additionally, numerous microsatellite-like sequences, such as (TA)n repeats, are present in the CR of all 34 species (Table S24). These sequences may facilitate the completion of DNA replication and transcription, and may serve as useful markers for geographical studies (Tian et al. 2023; Zhao et al. 2017). The extensive length polymorphism and presence of variable microsatellite motifs in the CR suggest this region may contain sufficient population-level variation to resolve phylogeographic structure within widely distributed Culex species complexes (Simon et al. 2006).

Secondary structures of all 22 tRNAs are predicted and compared for the seven newly sequenced mitogenomes (Fig. 4). Most tRNAs are highly conserved. The greatest difference is observed in trnD (11 variable sites), while trnQ, trnH, trnI, trnL1, trnK, trnM, trnS1, and trnM each exhibit only one variable site. Almost all tRNAs fold into typical cloverleaf structures, except trnS2, in which the dihydrouridine (DHU) arm forms a simple loop. This feature is widely documented in metazoan mitogenomes (Wolstenholme 1992).

The rrnL is located between trnL1 and trnV, with length ranging from 1,279 to 1,360 bp across the 34 species. A + T content varies from 82.49% to 83.76%, and G + C content from 16.24% to 17.51%. The content of T is higher than A, and the content of G is higher than C (Table S25). The predicted secondary structures of rrnL contain four canonical domains (I–II, IV–V) and 46 helices (Fig. 5), with domains III and VI absent, a typical structure in arthropods (Cannone et al. 2002).

The rrnS is located between trnV and the CR, with lengths ranging from 783 to 819 bp and A + T content from 80.46% to 81.89%. The content of T is higher than A, and the content of G was higher than C (Table S26). The predicted secondary structure of rrnS contains three domains and 34 helices (Fig. 6).

There ar e five different types of datasets used in the phylogenetic analysis (Fig. 7, Figs S3–S11). The PCG123RNA dataset (13,222 sites) incorporates all codon positions of 13 PCGs and two rRNAs; the PCG123 dataset (11,160 sites) contains three codon positions for all 13 PCGs; the PCG12RNA dataset (9,502 sites) includes only the first and second codon positions of 13 PCGs along with two rRNAs; the PCG12 dataset (7,440 sites) comprises only the first two codon positions from all 13 PCGs, while the dataset PCGAA contains 3,720 sites including all amino acids of 13 PCGs.

There are 10 phylogenetic trees reconstructed based on BI and ML analyses in this study (Figs 7, S3–S11). The phylogenetic results demonstrate that both the subfamily Anopheline and Culicinae are monophyletic.

Across all phylogenetic topologies examined, the genus Culex is consistently recovered non-monophyletic, with Lutzia fuscana invariably neste within it. With the exception of the ML analysis using the PCG123 dataset, all topologies support the monophyly of Cx. (Lophoceraomyia) and Cx. (Culiciomyia). In contrast, Cx. (Neoculex) and Cx. (Culex) are recovered as paraphyletic. Notably, the BI and ML analyses based on the PCGAA dataset reveal that subgenera Cx. (Oculeomyia) and Cx. (Culiciomyia) are nested within Cx. (Culex), which is consistent with ML analyses of all datasets except PCG12RNA. Based on the PCGAA dataset, both BI and ML phylogenetic analyses consistently support the monophyly of Cx. (Lophoceraomyia) and its position as the sister group to a clade containing Lutzia, Cx. (Culex), and several other subgenera nested within Cx. (Culex). This result is consistent with the ML analyses based on the PCG12, PCG123, and PCG123RNA datasets. The BI analysis of PCG12 dataset and both BI and ML analyses of PCG12RNA demonstrate that Cx. (Oculeomyia) is the sister group to Cx. (Culex) with Cx. (Culiciomyia) nested within Cx. (Culex), while BI analyses of PCG123 and PCG123RNA datasets consistently recover Cx. (Culiciomyia) being the sister group to Cx. (Culex) with both Lutzia and Cx. (Oculeomyia) nested within Cx. (Culex) clade.

The divergence time is estimated based on the phylogeny of dataset PCGAA and the chronogram is shown in Figure 8. Mean age values and 95% high posterior density (HPD) intervals for each node are presented in Table S27.

The family Culicidae diverged during Late Triassic (~219.34 Mya; 95% HPD = 129.57–305.38 Mya). Within this family, Culicinae diverged from Anophelinae in Late Jurassic (~147.45 Mya; 95% HPD = 87.03–216.38 Mya). The clade Sabethini + Mansoniini split from other lineages during mid-Cretaceous (~119.90 Mya; 95% HPD = 69.43–177.94 Mya), followed by the divergence of Mansoniini from Sabethini later in mid-Cretaceous (~98.67 Mya; 95% HPD = 49.99–153.91 Mya). The lineage comprising Aedeomyiini + Toxorhynchitini emerged during mid-Cretaceous (~106.60 Mya; 95% HPD = 61.96–161.42 Mya), while Culisetini + Culicini diverged from Aedini later in Cretaceous (~86.42 Mya; 95% HPD = 47.41–131.23 Mya). Significant Neogene divergences include the separation of Toxorhynchitini from Aedeomyiini in the Early Neogene (~22.00 Mya; 95% HPD = 14.23–33.20 Mya). Additionally, the split between Culicini and Culisetini occurred during the Late Cretaceous (~67.12 Mya; 95% HPD = 35.79–102.36 Mya).

Within Culex, the clade Cx. (Culex) + Lutzia diverged from Cx. (Lophoceraomyia) at ~28.62 Mya (95% HPD = 15.07–44.73 Mya), with Cx. (Culex) differentiating from Lutzia at ~22.69 Mya (95% HPD = 13.13–35.15 Mya). The subgenera Cx. (Oculeomyia) and Cx. (Culiciomyia) diverged at ~11.09 Mya (95% HPD = 4.54–18.23 Mya) and ~11.65 Mya (95% HPD = 6.61–17.71 Mya), respectively.