Copepods are omnipresent in aquatic ecosystems. They are extremely abundant in freshwater and constitute a major component of most planktonic, benthic and groundwater communities, including semi-terrestrial situations such as damp mosses and leaf litter on forest floors (Boxshall and Jaume, 2000; Boxshall and Defaye, 2008). Of an estimated 13,000 morphospecies known, about 2,800 species inhabit freshwater. Freshwater copepods are comprised of five orders (Calanoida, Harpacticoida, Cyclopoida, Gelyelloida, and Siphonostomatoida; Boxshall and Halsey, 2004). Among them, calanoids show interesting ecological characteristics (Andronov, 1974; Park, 1986; Blanco-Bercial et al., 2011). Diaptomidae (Sars, 1903) is the dominant family in inland waters of Europe, Asia, North America, Africa, and northern low-altitude South America.

Most freshwater Diaptomidae are planktonic, inhabiting the pelagic (Boxshall and Defaye, 2008). They are filter-feeding and an important regulator of bottom-up (available resources) and top-down (consumers) effects in freshwater zooplankton communities (Gliwicz, 2002; Lin et al., 2018). They respond quickly to environmental changes and most are endemic to well-defined regions (Matsumura-Tundisi and Tundisi, 2003; Boxshall and Defaye, 2008). As copepods have a high dispersal potential, it is interesting to know how much their populations differentiate within species’ ranges (Marrone et al., 2013). Molecular phylogenetic analysis with mitochondrial DNA has been used in a variety of studies including evolutionary genomics, systematics, and molecular evolution has become a robust and reliable tool due to such unique features of mitogenomes as high mutation rates, an absence of introns, and maternal inheritance (Ki et al., 2009). In calanoids, several partial mitochondrial sequences, such as Cytb, COI, and 16S rRNA, have been utilized for species diagnoses (Marrone et al., 2013; Perbiche-Neves et al., 2015; Previšić et al., 2016), phylogenetic analysis (Blanco-Bercial et al., 2011) and population genetics (Makino and Tanabe, 2009; Chen and Hares, 2011). For a phylogenetic relationship multiple gene analysis is more powerful than analysis using single markers (Blanco-Bercial et al., 2011). Complete mitochondrial sequences can provide improved resolution and sensitivity in such investigations. The complete mitogenomes of copepods are typically circular molecules, approximately 16kb in length (Jooste et al., 2019). Overall, they consist of 37 genes, including 13 protein coding genes (PCGs), two ribosomal RNA genes (rRNAs), 22 transfer RNA genes (tRNAs), and a control region.

To date, only two complete mitochondrial genomes have been isolated from Diaptomidae, a tiny fraction of the number of known species (Jooste et al., 2019; Zhang et al., 2020b). Here, we discuss another eight geographical mitogenomes of Phyllodiaptomus tunguidus, a freshwater species, widespread in a variety of waters including ponds, lakes and reservoirs in southern China (Chen, 1990; Dumont and Reddy, 1993; Xue et al., 2006). Shen and Tai (1964) first described it from Guangdong Province, southern China. Dumont and Reddy (1993) re-described the species from the River Li in Guangxi and Lake Erh in Yunnan Province, also southern China, and they characterized it as a Chinese endemic (Dumont and Reddy, 1993). Literature and our own field investigation show that most of P. tunguidus is currently distributed in the Pearl River Basin, but also in Hunan, Sichuan and Fujian provinces (Chen, 1990; Dumont and Reddy, 1993). From lowland (Guangdong) to highland (Yunnan), this includes a complex landscape with an elevation gradient of 2000m, mainly pertaining to the Pearl River Basin. This large area not only has a strong environmental gradient, but also a complex geological history, affected by the uplift of the Tibetan Plateau (Hall, 2002; Clark et al., 2004; Zheng, 2015). Therefore, a deep genetic differentiation is expected between at least the most extreme of its geographic populations. To document this, we here sequenced eight mitogenomes of far-apart geographical origin plus the mitogenome we previously sequenced (Zhang et al., 2020b). We searched for differences in genomic structure, base composition, codon usage, gene order, and other genetic information. We also conducted an extensive phylogenetic analysis that included 11 species from other copepod groups. Such a comparative analysis of whole mitogenomes was expected to provide new insights into the evolutionary relationships in Copepoda and to reaffirm the taxonomic position of P. tunguidus.

The mitogenomes of P. tunguidus were found to be similar to those of other copepods in terms of gene quantity and organizational structure (Ki et al., 2009; Jooste et al., 2019). Their characteristics conformed to the typical invertebrate mitogenome without gene insertion and deletion but with rearrangements. A common feature in metazoan mitochondrial genomes is bias toward a higher representation of nucleotides As and Ts leading to a subsequent bias in corresponding encoded amino acids (Wang et al., 2018; Ding et al., 2020; Liu et al., 2020; Sun et al., 2020). At this point, DNA with a high AT-content makes the double helix less stable because the AT pair is bound by two hydrogen bonds while GC pairs are bound by three such bonds (Yakovchuk et al., 2006); therefore, the high AT-content of mitogenome could have a higher mutation rate such that may accelerate the evolutionary process for Diaptomidae (Triant and DeWoody, 2006; Ding et al., 2020). Interestingly, the A+T content of PCGs at the 3rd codon position is higher than that of 1st and 2nd codon positions, and the G content at the 3rd codon position was lower than that 1st and 2nd codon positions. In addition, bias against T is prominent in the whole genome, all PCGs, the first, the second and the third codon of PCGs, especially at the second codon. A possible explanation is that the AT bias can selectively avoid the formation of stop codons and the loss of amino acids, and results in T-rich codons (Charneski et al., 2011; Ding et al., 2020). We determine preference for specific synonymous codons by evaluating the relative synonymous codon usage (RSCU; Granthan et al., 1980; Meganathan et al., 2012). The codon usage patterns were similar among the four P. tunguidus, with both 2- and 4-fold degenerate codons exhibiting over-usage of A and T at the third codon positions. The over-usage of A and T at the third codon is consistent with the observation in other animal species (Mindell et al., 1998; Ding et al., 2020; Sun et al., 2020). This might be related to genome bias, optimal selection of tRNA usage, or DNA repair efficiency (Crozier and Crozier, 1993).

The Ka/Ks ratio is an indicator of selective pressure acting on PCG (Yang and Bielawski, 2000; Hurst, 2002). It changed in the 13 PCGs, reflecting different functional constraints between genes (Muse, 2000). Specifically, COI had the lowest Ka/Ks ratio, indicating that this gene was under the strongest selective pressure. The ND3 gene had the highest Ka/Ks ratio, i.e., it was under the least selective pressure (Yang and Nielsen, 2000; Sun et al., 2011). Molecular evolution (Ka/Ks ratio) of single PCGs seems to be different in all nine P. tunguidus. However, it was consistent for the two P. tunguidus from NB and TB in Guangxi Province (near Vietnam). We failed to detect significant adaptive mutations in the nine mitogenomes. Therefore, purifying selection may be the predominant force governing mitogenome evolution in the species. Previous studies showed that energetic constraints are the main factors contributing to the evolution of mitochondrial-encoded proteins (Sun et al., 2011; Wagner et al., 2018). Combining the ratio of Ka/Ks and the genetic distance, we considered the P. tunguidus from NB and TB a unique lineage compared to the seven from the other sites, and the PCGs gene plays a crucial role in the evolutionary adaption under selection.

The ATP8 gene had the smallest genetic distance among the 13 single PCGs, whereas the ND4 gene had the largest distance, revealing different mutation pressures among genes (Sun et al., 2011). Additionally, the pair distance based on amino acid sequences was found to be higher than those based on nucleotides sequences. In other words, there are fewer synonymous substitutions than non-synonymous substitution in PCGs of the P. tunguidus mitogenomes, suggesting that some PCGs have undergone positive selection.

There is no evidence of resting eggs from P. tunguidus until now, thus it may be passively dispersed mainly following water flow. Isolation by distance (IBD) is expected in the dispersal limitation. However, the result showed no significantly correlation between the geographical distance and genetic differentiation (mitogenome, amino acid, PCGs, ribosomal RNA). This spatial pattern without IBD indicated that the current geographic distribution of P. tunguidus could have a complex historical origin. The geographical topography is indeed extremely complex in southern China (Clark et al., 2004; Xiang et al., 2007; Zheng, 2015), where Yunnan-Guizhou Plateau, Yunkai mountain and Shiwan mountain structure the landscape for P. tunguidus. Genetic differentiation between adjacent populations such as GZR and GZP, LXH and ZWY indicated that there are strong environmental filters or selections.

The exact position of 22 tRNA genes of the eight P. tunguidus was identified as being the same genomic position observed for our published P. tunguidus mitogenome (Zhang et al., 2020b). The tRNA gene in all the nine P. tunguidus had anti-codons that match the invertebrate mitochondrial code (Weydmann et al., 2017). However, tRNA-Arg of the two P. tunguidus from NB and TB in Guangxi Province had the complete clover-leaf structure. The structure of tRNA-Ser1 was not the typical clover-leaf structure that occurs in other invertebrates and vertebrates (Wolstenholme, 1992; Hong et al., 2008; Yang et al., 2019). The secondary structure of lacking a dihydrouridine arm is likely related its structural compensation mechanism among tRNA arms (Steinberg and Cedergren, 1994). Mismatched and wobble pairs are commonly found in invertebrate tRNAs and can be corrected by post-transcriptional RNA editing processes (Lavrov et al., 2000). Two rRNA genes were identified for P. tunguidus through alignment with other Diaptomidae mitogenomes (Jooste et al., 2019; Zhang et al., 2020b). The rRNA genes exhibited a heavy AT nucleotide bias in the P. tunguidus mitogenomes. The RNA-Seq analysis in many species demonstrates that a TAA stop codon is created by post-transcriptional polyadenylation that changes T or TA residues to a complete TAA stop codon. Therefore, the rRNA with a high AT-content species has a less stable double helix, which helps to implement the above mechanisms (Yakovchuk et al., 2006).

Phylogenetic analysis indicates that Copepoda species fall into three orders. Calanoida and Cyclpoida for a monophyletic clade, but Harpacticoida do not. Calanoida is identified as the most recently diverging order of the three, and Diaptomidae is the most recently diverging family of the Calanoida. Morphological taxonomic studies are generally consistent the current molecular phylogenies except for the A. atopus in Harpacticoida.

In this study, we clarified the phylogeny and hidden diversity of P. tunguidus in southern China. Focusing on the mitochondrial genomes of nine P. tunguidus, the phylogeny showed that the two P. tunguidus from NB and TB (Clade 1) are separated from the other seven. A phylogenetic tree based on mitogenomes and Cytb showed that Clade 1 is restricted to the southern Guangxi Province. In addition, two P. tunguidus populations (TB and NB) near Vietnam can be considered older, if we look at the fact that most species in Phyllodiaptomus are limited to Southeast Asia (Dumont and Reddy, 1993, 1994; Dumont et al., 1996; Sanoamuang and Yindee, 2001; Sanoamuang and Teeramaethee, 2006; Sanoamuang and Watiroyram, 2020). In contrast, Clade 2 is widely distributed in many relatively isolated basins, including the Lancnagjiang River (Yunnan Dali), Upper Pearl River (Yunnan, Kunming), Lower Pearl River (Gaozhou, Guangzhou and Shaoguan), the Minjiang River (Fujian, Fuzhou) and the Yangtze River (Hunan, Changsha; Supplementary Table S2 and Supplementary Figure S1). The spatial proximity among these basins allows P. tunguidus populations in the Pearl River basin to expand into nearby basins easily, especially in flooding seasons. At the same time, the hypothesis of active river capture in southwestern China may explain the close relationship of populations between the Yangtze and Upper Pearl River (Clark et al., 2004; Zheng, 2015). In addition to the two strongly supported mitogenome Clades, a high level of mtDNA genetic differentiation, observed between the two Clades (5.2–13.4%) supports a deep divergence among P. tunguidus populations. This finding is shared with several small-sized freshwater fish species in China, e.g., Hemiculter leucisculus (Chen et al., 2017), Opsariichthys bidens (Perdices et al., 2005), Rhynchocypris oxycephalus (Yu et al., 2014).

Similar to the phylogenetic relationship, the gene arrangement of mitochondrial genomes demonstrates that Calanoida have a more complete classification system, while the classification of Harpacticoida is currently unresolved.

The nine mitogenomes of P. tunguidus have similar size and gene arrangement, but their base composition, genetic distance and tRNA structure indicate a large genetic differentiation between mitogenomes. Deep genetic differentiation occurs between two clades of widely different geographic origins. This differentiation is also shown in the difference in the structures of tRNA and rRNA. Two clades have reached the level of subspecies or cryptic species. The mitogenome-based phylogeny of Copepoda further shows that Diaptomidae is the most recently diverging family of Calanoida and that P. tunguidus is at the evolutionary apex of the family.

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

B-PH: conceptualization. X-LZ, PL, and S-LX: methodology and software. X-LZ: formal analysis and writing – original draft preparation. X-LZ, QZ, and PL: data curation. X-LZ, B-PH, ER, and HD: writing – review and editing. All authors have read and agreed to the published version of the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.