# Phylogenetic Information Content of Copepoda Ribosomal DNA Repeat Units: ITS1 and ITS2 Impact

Authors: Maxim V. Zagoskin; Valentina I. Lazareva; Andrey K. Grishanin; Dmitry V. Mukha
Journal: BioMed Research International
 (2014)

---

## Abstract

The utility of various regions of the ribosomal repeat unit for phylogenetic analysis was examined in 16 species representing four families, nine genera, and two orders of the subclass Copepoda (Crustacea). Fragments approximately 2000 bp in length containing the ribosomal DNA (rDNA) 18S and 28S gene fragments, the 5.8S gene, and the internal transcribed spacer regions I and II (ITS1 and ITS2) were amplified and analyzed. The DAMBE (Data Analysis in Molecular Biology and Evolution) software was used to analyze the saturation of nucleotide substitutions; this test revealed the suitability of both the 28S gene fragment and the ITS1/ITS2 rDNA regions for the reconstruction of phylogenetic trees. Distance (minimum evolution) and probabilistic (maximum likelihood, Bayesian) analyses of the data revealed that the 28S rDNA and the ITS1 and ITS2 regions are informative markers for inferring phylogenetic relationships among families of copepods and within the Cyclopidae family and associated genera. Split-graph analysis of concatenated ITS1/ITS2 rDNA regions of cyclopoid copepods suggested that theMesocyclops, Thermocyclops, and Macrocyclops genera share complex evolutionary relationships. This study revealed that the ITS1 and ITS2 regions potentially represent different phylogenetic signals.


---

## Body

## 1. Introduction

Copepods are important components of zooplankton and the food chain in marine and freshwater ecosystems. The subclass Copepoda is believed to contain approximately 13,000 morphospecies; however, the actual number of species in this subclass might be much greater [1]. The majority of the freshwater copepod species belong to the order Cyclopoida, which includes the free-living species (approximately 800) in the family Cyclopidae. The other two free-living families (Oithonidae and Cyclopinidae) contain mainly marine species except for a few species in Oithonidae [2].Systematic analyses of cyclopoid copepods (order Cyclopoida) have primarily focused on morphological characteristics [2–8], and the majority of molecular studies have targeted marine copepods [9–41]. The phylogenetic history of freshwater cyclopoid copepods is not well understood. A few studies on Cyclopidae have used molecular and morphological analyses on theMesocyclops genus (Crustacea: Cyclopidae) [42],E. serrulatus group [43],Acanthocyclops vernalis-robustus species complex [44, 45], and 11 populations ofMacrocyclops albidus [46]; other phylogenetic analyses have focused only on molecular markers inDiacyclops spp., which are found in western Australia [47] and Lake Baikal [48].Molecular markers such as genomic DNA fragments are used for phylogenetic analyses to elucidate the evolutionary history of living organisms, and the region of genomic DNA analyzed is critical. Mitochondrial DNA fragments (genes encoding the cytochrome c oxidase subunit I (COI), 16S rRNA, and cytochrome b) [9–21, 39, 40, 49–52] and/or nuclear rDNA regions have been used for the phylogenetic analysis of cyclopoid copepods [22–37, 41, 53–58]. Mitochondrial DNA fragments might be less useful for the analysis of copepod phylogeny compared to the phylogeny of other taxa; furthermore, amplification of COI is difficult in some copepods [59–61]. However, these DNA fragments might be informative for analyses of population differentiation or cryptic speciation [9, 11–15, 38–40, 50, 52]. Therefore, comparison of nuclear rDNA regions might be informative for the phylogenetic analysis of copepods.In most eukaryotes, rRNA genes are located in a multigene family of genomic clusters of repeated sequences. Within these clusters, the 18S, 5.8S, and 28S rRNA genes are separated by internal transcribed spacers (ITS1 and ITS2) and an intergenic spacer [62] (Figure 1). Ribosomal DNA (rDNA) is a reliable and informative phylogenetic marker [63] that contains sequences with different rates of evolutionary variability. In most eukaryotes, the most evolutionarily conserved genes are the rRNA genes; comparison of their sequences allows estimation of the evolutionary distances at intergeneric and higher taxonomic levels. Comparison of the more evolutionarily variable spacer sequences enables the study of phylogenetic relationships at the species and population levels [63–65]. Therefore, comparison of different regions of rDNA enables the phylogenetic analysis of organisms over extended evolutionary distances.Figure 1
Organization of eukaryotic tandemly repeated rDNA clusters. 18S, 5.8S, and 28S ribosomal RNA genes; ITS1 and ITS2 internal transcribed spacers; IGS intergenic spacer. Arrows indicate the locations of the DAMS18 and DAMS28 primers.Phylogenetic relationships among cyclopoid copepods at higher systematic levels (ordinal, familial, and generic) have been resolved using the 18S and 28S nuclear rRNA genes [28, 29, 55–57], and the relationships at the lower taxonomic levels (species and populations) have been resolved using the ITS2 of the nuclear rDNA gene cluster [30, 40, 42, 52, 58].Notably, analysis of the evolutionary history of living organisms based on only one molecular marker can uncover bifurcating phylogenetic trees, revealing branched evolution. However, it recently becomes evident that evolution is not always tree-like. Comparisons of gene trees based on different genetic loci often reveal conflicting tree topologies. These discrepancies are not always due to the problems with the sampling and the gene tree reconstruction methods. Reticulation events such as horizontal gene transfer (HGT) and hybridization may be responsible for contradictions in lineages. During an HGT event, a DNA segment is transferred from one organism to another which is not its offspring, whereas hybridization describes the origin of a new species through an interspecies mating. Both processes yield genomes that are mixtures of DNA regions derived from different species. Consequently, evolutionary relationships between species whose past includes reticulation can often be better represented by using phylogenetic networks rather than trees [65–67].In view of the above, comparing phylogenetic trees based on different molecular markers may be used for the analysis of evolutionary events caused by reticulate evolution. Phylogenetic signals from various molecular markers are potentially divergent during reticulate evolution, resulting in phylogenetic trees with alternative positions for the individual branches [68–71]. Comparative analysis of the rDNA ITS1 and ITS2 sequences is suitable for studying phylogenetic relationships in terms of branching and reticulate evolution [63, 68, 70, 72–82].Reticulate evolution is primarily driven by hybrid speciation, which is common among plants [83] but also occurs among animals, particularly including fish [84], amphibians, and several invertebrates [85–87]. In both mammals [88] and arthropods [89, 90], a single instance of hybrid speciation has been well described. Interspecies hybridization typically results in complicated relationships within species complexes, characterized by indistinct species borders. Reticulate evolution among crustaceans has been observed only within species complexes of daphnids [91–94].In this study, we analyzed the phylogenetic relationships within a small group of cyclopoid copepods representing several genera of freshwater (Cyclops,Thermocyclops,Diacyclops,Megacyclops,Macrocyclops, andMesocyclops) and marine (Oithona andParacyclopina) organisms. Specific freshwater species were selected for analysis because these species are important for the maintenance of food chains in Russian freshwater ecosystems. The aim of this study was to analyze the sequence characteristics of the rDNA 28S gene, ITS1, and ITS2 regions as phylogenetic markers for the selected group of organisms.


## 2. Materials and Methods

### 2.1. Samples Collection

Nine freshwater species of the Cyclopidae family were collected near the Borok settlement in the Yaroslavskaya region of Russia:Mesocyclops leuckarti (Claus, 1857),Cyclops strenuus (Fischer, 1851), andCyclops insignis (Claus, 1857) (population no. 1) from the Barskiy Pond (58°3′59.35′′N; 38°15′10.16′′E);Thermocyclops oithonoides (Sars, 1863) from the Sunoga pond (58°2′34.66′′N; 38°14′41.29′′E);Thermocyclops crassus (Fischer, 1853),Macrocyclops distinctus (Richard, 1887),Macrocyclops albidus (Jurine, 1820),Diacyclops bicuspidatus (Claus, 1857), andMegacyclops viridis (Jurine, 1820) (population no. 1) from the Ikhteologichesky Canal (58°3′55.62′′N; 38°15′21.05′′E); andMegacyclops viridis (Jurine, 1820) (population no. 2) from a pond in the flood zone of the Rybinsk Reservoir (58°4′4.70′′N; 38°15′39.88′′E).Cyclops kolensis (Lilljeborg, 1901) andCyclops insignis (Claus, 1857) (population no. 2) were collected from the Andreevsky small pond in Vorob’evy Gory, Moscow, Russia (55°42′35.40′′N; 37°34′6.61′′E). Two marine species,Oncaea sp. (Claus) andOithona similis(Claus, 1866), were collected from the Norwegian Sea (68°52′36.67′′N; 3°8′21.91′′E).Individuals of each species were collected for further analysis at the specified locations over a 0.5- to 1-hour period.No specific permission was required to collect samples at these locations. None of the studied species is endangered or protected.A fragment of the 28S gene from each of the four marine speciesParacyclopina nana (Smirnov, 1935) (GenBank accession number FJ214952),Oithona nana (Giesbrecht, 1893) (GenBank accession number FM991727),Oithona simplex (Farran, 1913) (GenBank accession number AF385458), andOithona helgolandica (Claus, 1863) (GenBank accession number FM991724.1) was also used for the analysis. The DNA was extracted from either samples preserved in 70% ethanol or raw materials (Moscow populations).


### 2.2. DNA Extraction, PCR Amplification, and Sequencing

The genomic DNA was isolated from 10–20 individuals of each collected species using the DNeasy Blood & Tissue kit (QIAGEN, Hilden, Germany) according to the manufacturer’s instructions and was frozen at −20°C. The rDNA region (approximately 2000 bp) was amplified from the genomic DNA by polymerase chain reaction (PCR) using the universal eukaryotic rDNA primers DAMS18 and DAMS28 [95–97] (Figure 1). The amplified rDNA regions contained the ITS1 (261–388 bp) and ITS2 (188–262 bp) regions, the 5.8S gene (157 bp), and approximately 200 and 1000 bp of the 18S and 28s genes, respectively. The amplification was performed in 50 μL reactions using a PCR Master Mix (2X) (Fermentas, Vilnius, Lithuania) according to the manufacturer’s instructions; the reactions were performed in a Primus 25 advanced Thermocycler (PEQLAB, Erlangen, Germany) using previously published rDNA-specific parameters [98]. The PCR products were resolved on 1.0% agarose gels, and DNA was extracted from the observed unique bands using the QIAquick Gel Extraction Kit (QIAGEN, Hilden, Germany). The extracted products were cloned into the pGEM-T Easy vector (Promega, USA), and the resulting plasmids were used to transformEscherichia coli JM109 competent cells (Promega, USA) according to the manufacturer’s instructions. For each species, the amplified product and five clones were sequenced. Automated sequences were generated on an ABI PRISM 310 Genetic Analyzer according to Sanger et al. [99] with a BigDye Termination kit (Applied Biosystems, USA). The sequences generated in this study were deposited in GenBank under the accession numbers KF153689–KF153701.


### 2.3. Phylogenetic Analyses

The rDNA sequences were aligned using ClustalW 2.1 [100, 101] with some manual adjustments. The boundaries of the ITS1 and ITS2 regions and the 28S gene were identified by comparing the primer-delimited sequences against sequences in the GenBank database using BLAST analysis. The boundaries of the conserved sequences were considered to represent the 5.8S, 18S, and 28S gene flanking regions if they were 100% similar to the boundaries of rDNA sequences in the GenBank database. The initial sequence alignment flanked by DAMS18/DAMS28 primers was divided into ITS1 and ITS2 alignments. The rDNA genetic distances were estimated using the MEGA V5.2 software [102]. DAMBE (Data Analysis in Molecular Biology and Evolution) software was used to analyze substitution saturation [103–105]. This method computes the entropy-based index of substitution saturation and its critical value. If the index of substitution saturation (Iss) approaches 1 or if the Iss is not smaller than the critical Iss value (Iss.c), then sequences are considered to contain substantial saturation. As is known, the substitution saturation decreases phylogenetic information contained in sequences and has plagued the phylogenetic analysis involving deep branches. In the extreme case when sequences have experienced full substitution saturation, the similarity between the sequences will depend entirely on the similarity in nucleotide frequencies, which often does not reflect phylogenetic relationships [106].The rDNA-based phylogenetic trees were estimated using probabilistic (maximum likelihood (ML), Bayesian) and distance (minimum evolution (ME)) methods [107–110]. ML and ME analyses of ITS1, ITS2, and 28S data were performed using the program MEGA V5.2. Branch support was assessed using the bootstrap method [111] (1,000 replicates) with the close-neighbor-interchange (CNI) algorithm at a search level of 1 for ME analysis and heuristic search for ML analysis. The Bayesian information criterion (BIC), as implemented in MEGA V5.2, was used to identify the best-fit model of sequence evolution for the trees estimated using ML. The evolutionary history was inferred using the ML method based on the general time reversible with the gamma distribution shape parameter (GTR+G) model for 28S and the Hasegawa-Kishino-Yano with gamma distribution shape parameter (HKY+G) model [112] for the ITS1, ITS2, and concatenated ITS1/ITS2 alignments. In addition to these methods, ITS1 and ITS2 alignments were constructed using the MAFFT version 7 (http://mafft.cbrc.jp/alignment/server/) [113] and Gblocks version 0.91b (http://www.phylogeny.fr/version2_cgi/one_task.cgi?task_type=gblocks) software programs [114–116] to eliminate poorly aligned and highly divergent regions. Default parameters were used for both of these methods. The Tamura 3-parameter model (T92) [117] and HKY with evolutionary invariable (HKY+I) for Gblocks-treated MAFFT ITS1 and ITS2 data, respectively, were used to infer evolutionary history inference using the ML method.The Bayesian analysis was performed using MrBayes version 3.1.2 software [118, 119]. Two replicate analyses of 1 million generations each were performed for each dataset, with sampling every 10 generations. The hierarchical likelihood ratio test (hLRT) implemented in MrModeltest version 2.3 software [120] was used to identify the model of best fit (Hasegawa-Kishino-Yano with invariant sites and gamma distribution shape parameter (HKY+I+G) [112] for ITS1 and the HKY+G model for ITS2). Trees from the first 53,000 and 118,000 generations were discarded as burn-in for ITS1 and ITS2, respectively. The Bayesian tree was estimated from the majority-rule consensus of the post-burn-in trees.A reticulogram [121] was constructed using the T-REX version 4.01a software [122] with the distance matrix computed using the Kimura 2-parameters model (ignoring missing bases); the weighted least-squares method was used for tree reconstruction [123], and addition of reticulation branches stopped when K=1 branches were added.Network reconstruction was performed using Splits Tree 4 version 4.11.3 software [65]. The neighbor-net network method and uncorrectedp-distances were used to analyze and visualize reticulate relationships. All gaps were excluded for analysis. Network robustness was tested using 1,000 bootstrap replicates.


## 2.1. Samples Collection

Nine freshwater species of the Cyclopidae family were collected near the Borok settlement in the Yaroslavskaya region of Russia:Mesocyclops leuckarti (Claus, 1857),Cyclops strenuus (Fischer, 1851), andCyclops insignis (Claus, 1857) (population no. 1) from the Barskiy Pond (58°3′59.35′′N; 38°15′10.16′′E);Thermocyclops oithonoides (Sars, 1863) from the Sunoga pond (58°2′34.66′′N; 38°14′41.29′′E);Thermocyclops crassus (Fischer, 1853),Macrocyclops distinctus (Richard, 1887),Macrocyclops albidus (Jurine, 1820),Diacyclops bicuspidatus (Claus, 1857), andMegacyclops viridis (Jurine, 1820) (population no. 1) from the Ikhteologichesky Canal (58°3′55.62′′N; 38°15′21.05′′E); andMegacyclops viridis (Jurine, 1820) (population no. 2) from a pond in the flood zone of the Rybinsk Reservoir (58°4′4.70′′N; 38°15′39.88′′E).Cyclops kolensis (Lilljeborg, 1901) andCyclops insignis (Claus, 1857) (population no. 2) were collected from the Andreevsky small pond in Vorob’evy Gory, Moscow, Russia (55°42′35.40′′N; 37°34′6.61′′E). Two marine species,Oncaea sp. (Claus) andOithona similis(Claus, 1866), were collected from the Norwegian Sea (68°52′36.67′′N; 3°8′21.91′′E).Individuals of each species were collected for further analysis at the specified locations over a 0.5- to 1-hour period.No specific permission was required to collect samples at these locations. None of the studied species is endangered or protected.A fragment of the 28S gene from each of the four marine speciesParacyclopina nana (Smirnov, 1935) (GenBank accession number FJ214952),Oithona nana (Giesbrecht, 1893) (GenBank accession number FM991727),Oithona simplex (Farran, 1913) (GenBank accession number AF385458), andOithona helgolandica (Claus, 1863) (GenBank accession number FM991724.1) was also used for the analysis. The DNA was extracted from either samples preserved in 70% ethanol or raw materials (Moscow populations).


## 2.2. DNA Extraction, PCR Amplification, and Sequencing

The genomic DNA was isolated from 10–20 individuals of each collected species using the DNeasy Blood & Tissue kit (QIAGEN, Hilden, Germany) according to the manufacturer’s instructions and was frozen at −20°C. The rDNA region (approximately 2000 bp) was amplified from the genomic DNA by polymerase chain reaction (PCR) using the universal eukaryotic rDNA primers DAMS18 and DAMS28 [95–97] (Figure 1). The amplified rDNA regions contained the ITS1 (261–388 bp) and ITS2 (188–262 bp) regions, the 5.8S gene (157 bp), and approximately 200 and 1000 bp of the 18S and 28s genes, respectively. The amplification was performed in 50 μL reactions using a PCR Master Mix (2X) (Fermentas, Vilnius, Lithuania) according to the manufacturer’s instructions; the reactions were performed in a Primus 25 advanced Thermocycler (PEQLAB, Erlangen, Germany) using previously published rDNA-specific parameters [98]. The PCR products were resolved on 1.0% agarose gels, and DNA was extracted from the observed unique bands using the QIAquick Gel Extraction Kit (QIAGEN, Hilden, Germany). The extracted products were cloned into the pGEM-T Easy vector (Promega, USA), and the resulting plasmids were used to transformEscherichia coli JM109 competent cells (Promega, USA) according to the manufacturer’s instructions. For each species, the amplified product and five clones were sequenced. Automated sequences were generated on an ABI PRISM 310 Genetic Analyzer according to Sanger et al. [99] with a BigDye Termination kit (Applied Biosystems, USA). The sequences generated in this study were deposited in GenBank under the accession numbers KF153689–KF153701.


## 2.3. Phylogenetic Analyses

The rDNA sequences were aligned using ClustalW 2.1 [100, 101] with some manual adjustments. The boundaries of the ITS1 and ITS2 regions and the 28S gene were identified by comparing the primer-delimited sequences against sequences in the GenBank database using BLAST analysis. The boundaries of the conserved sequences were considered to represent the 5.8S, 18S, and 28S gene flanking regions if they were 100% similar to the boundaries of rDNA sequences in the GenBank database. The initial sequence alignment flanked by DAMS18/DAMS28 primers was divided into ITS1 and ITS2 alignments. The rDNA genetic distances were estimated using the MEGA V5.2 software [102]. DAMBE (Data Analysis in Molecular Biology and Evolution) software was used to analyze substitution saturation [103–105]. This method computes the entropy-based index of substitution saturation and its critical value. If the index of substitution saturation (Iss) approaches 1 or if the Iss is not smaller than the critical Iss value (Iss.c), then sequences are considered to contain substantial saturation. As is known, the substitution saturation decreases phylogenetic information contained in sequences and has plagued the phylogenetic analysis involving deep branches. In the extreme case when sequences have experienced full substitution saturation, the similarity between the sequences will depend entirely on the similarity in nucleotide frequencies, which often does not reflect phylogenetic relationships [106].The rDNA-based phylogenetic trees were estimated using probabilistic (maximum likelihood (ML), Bayesian) and distance (minimum evolution (ME)) methods [107–110]. ML and ME analyses of ITS1, ITS2, and 28S data were performed using the program MEGA V5.2. Branch support was assessed using the bootstrap method [111] (1,000 replicates) with the close-neighbor-interchange (CNI) algorithm at a search level of 1 for ME analysis and heuristic search for ML analysis. The Bayesian information criterion (BIC), as implemented in MEGA V5.2, was used to identify the best-fit model of sequence evolution for the trees estimated using ML. The evolutionary history was inferred using the ML method based on the general time reversible with the gamma distribution shape parameter (GTR+G) model for 28S and the Hasegawa-Kishino-Yano with gamma distribution shape parameter (HKY+G) model [112] for the ITS1, ITS2, and concatenated ITS1/ITS2 alignments. In addition to these methods, ITS1 and ITS2 alignments were constructed using the MAFFT version 7 (http://mafft.cbrc.jp/alignment/server/) [113] and Gblocks version 0.91b (http://www.phylogeny.fr/version2_cgi/one_task.cgi?task_type=gblocks) software programs [114–116] to eliminate poorly aligned and highly divergent regions. Default parameters were used for both of these methods. The Tamura 3-parameter model (T92) [117] and HKY with evolutionary invariable (HKY+I) for Gblocks-treated MAFFT ITS1 and ITS2 data, respectively, were used to infer evolutionary history inference using the ML method.The Bayesian analysis was performed using MrBayes version 3.1.2 software [118, 119]. Two replicate analyses of 1 million generations each were performed for each dataset, with sampling every 10 generations. The hierarchical likelihood ratio test (hLRT) implemented in MrModeltest version 2.3 software [120] was used to identify the model of best fit (Hasegawa-Kishino-Yano with invariant sites and gamma distribution shape parameter (HKY+I+G) [112] for ITS1 and the HKY+G model for ITS2). Trees from the first 53,000 and 118,000 generations were discarded as burn-in for ITS1 and ITS2, respectively. The Bayesian tree was estimated from the majority-rule consensus of the post-burn-in trees.A reticulogram [121] was constructed using the T-REX version 4.01a software [122] with the distance matrix computed using the Kimura 2-parameters model (ignoring missing bases); the weighted least-squares method was used for tree reconstruction [123], and addition of reticulation branches stopped when K=1 branches were added.Network reconstruction was performed using Splits Tree 4 version 4.11.3 software [65]. The neighbor-net network method and uncorrectedp-distances were used to analyze and visualize reticulate relationships. All gaps were excluded for analysis. Network robustness was tested using 1,000 bootstrap replicates.


## 3. Results and Discussion

### 3.1. Characteristics and Analysis of an rDNA Sequence Dataset

In each species, the nucleotide sequences of the amplified rDNA region and five clones were not significantly different from each other. The frequency of variable nucleotides did not exceed the average rate of nucleotide substitutions caused by DNA polymerase errors, which is approximately one substitution per 1,000 nucleotides. The compared sequences contained both relatively evolutionarily conserved (fragments of 18S and 28S rDNA and the complete 5.8S rDNA) and evolutionary variable genomic regions (ITS1 and ITS2). For different taxa, the ITS1 and ITS2 sequences vary significantly among individuals at the inter- and intrapopulation levels; furthermore, these sequences can exhibit intragenomic variability [25, 41, 53, 54]. Recently, a high level of intrapopulation polymorphism of the 28S rDNA sequences was observed withinOithona spp. [22]. However, there are instances of strong evolutionary conservation of the 28S and ITS sequences [15, 23, 30, 34]. Notably, theM. leuckarti ITS2 sequence obtained in this study did not exhibit any nucleotide substitutions compared to theM. leuckarti ITS2 sequence described previously (GenBank accession number GQ848499) [42]. Therefore, the strong evolutionary conservation of ITS1 and ITS2 sequences is a characteristic feature of the copepod species analyzed in this study.In this study, the applicability of different segments of rDNA containing the ITS1 and ITS2 regions, the 5.8S RNA gene, and fragments of the 18S and 28S rRNA genes was examined for reconstruction of the phylogenetic relationships among freshwater cyclopoid copepods. The 5.8S gene and the analyzed fragment of the 18S gene were not considered for phylogenetic reconstruction due to their short length and strong evolutionary conservation: only a few nucleotide substitutions were detected by comparing these sequences with evolutionary distant species (data not shown).For 15 specimens of Cyclopoida species (including the two marine species), the average length of the 28S gene fragment sequenced was 1051 bp. We trimmed these sequences to 703 bp and compared them with the 28S gene sequences of marine species available in GenBank. These 703 bp of 28S rDNA sequences were aligned, and 342 variable sites were observed.Oncaea sp. (Oncaeidae family) was used as the out group.The ITS1 sequence lengths varied from 267 to 388 bp among the 13 Cyclopidae specimens. The ITS1 sequence alignments possessed 442 characters, and among them, 283 were variable. The ITS2 sequence lengths varied from 188 to 262 bp among the 13 Cyclopidae specimens. ITS2 sequence alignment possessed 302 characters, and 190 were variable.All alignment sets were examined for homogeneity of base frequencies and substitution saturation. The average base frequencies of the 28S gene fragment (A=20.52, C=25.33, G=32.75, and T=21.40%) differed from the ITS1 (A=14.27, C=30.68, G=27.55, and T=27.50%) and ITS2 (A=13.08, C=29.64, G=30.08, and T=27.19%) regions. Gaps were excluded while estimating the average base frequencies of the ITS sequences. Using the chi-squared test, no significant differences were observed in the base compositions of the 28S (χ2 = 39.77, df=54, and P=0.93), ITS1 (χ2=22.04, df=36, and P=0.97), and ITS2 (χ2=22.65, df=36, and P=0.96) sequences among different taxa.To analyze whether the divergence of 28S, ITS1, and ITS2 rDNA fragments among species was saturated, we performed a substitution saturation test and generated saturation plots. Using DAMBE, the substitution saturation test revealed an Iss value that was significantly(P<0.001) lower than the Iss.c in all cases (Table 1). This result indicated the suitability of the data for phylogenetic analysis. The total numbers of transition and transversion substitutions were plotted individually against model-corrected maximum-likelihood pairwise distances for the 28S, ITS1, and ITS2 sequences (see Supplementary Figure 1 in Supplementary Material available online at http://dx.doi.org/10.1155/2014/926342). Using linear regression analysis on the 28S, ITS1, and ITS2 saturation graphs, the coefficients of determination (R2) were calculated for both classes of substitutions: for transitions, the R2 values were 0.79, 0.68, and 0.74 for the 28S, ITS1, and ITS2 sequences, respectively; for transversions, the R2 values were 0.95, 0.93, and 0.91 for the 28S, ITS1, and ITS2 sequences, respectively. The R2 values indicated that no less than 70% of the total variation in pairwise transitions and transversions could be explained by the linear relationship between pairwise distances and the total number of transitions and transversions. All saturation plots showed significant linear correlations (Supplementary Figure 1). Therefore, both transitions and transversions steadily accumulated as the corrected pairwise divergence increased, indicating that saturation was not reached.Table 1
False test of substitution saturation.
Alignment
Iss
Iss.c
Std. error
28S
0.200
0.739
0.019
ITS1
0.427
0.679
0.044
ITS2
0.376
0.665
0.043
Testing whether the observed Iss is significantly(P<0.001, two-tailed t-test) lower than the Iss.c for a symmetrical tree.


### 3.2. Distance Analyses and Phylogenetic Tree Reconstruction

Phylogenetic analysis of the cyclopoid copepods species based on rDNA showed that the 28S rDNA sequences are informative for the phylogeny of both higher-level and closely related Copepoda species, whereas the ITS1 and ITS2 sequences are highly informative for reconstruction of the evolutionary history of closely related species. The ITS1 and ITS2 sequences are known to evolve more rapidly than the ribosomal RNA genes. Consistent with this observation, in this study, the pairwise ITS1/ITS2p-distances were significantly higher thanthe 28Sp-distances (compare Tables 2 and 3). These data are consistent with other studies showing considerable variation in ITS1 and ITS2 divergence levels among different groups of copepods [40, 42, 52, 124, 125]. In this study, fragments of 28S rDNA sequences were used for the analysis of marine and freshwater cyclopoid copepods species, whereas ITS1 and ITS2 sequences were used exclusively for the analysis of freshwater cyclopoid copepods species.Table 2
Pairwise uncorrected genetic distance among 18 sequences of Cyclopoida based on comparison of 28S rRNA gene.
#
Species name
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
1
Oncaeasp.
2
Paracyclopina nana
0.229
3
Oithona similis
0.230
0.148
4
Oithona helgolandica
0.227
0.143
0.011
5
Oithona nana
0.233
0.181
0.155
0.146
6
Oithona simplex
0.244
0.152
0.175
0.166
0.177
7
Cyclops kolensis
0.209
0.223
0.197
0.195
0.244
0.212
8
Cyclops strenuus
0.207
0.220
0.197
0.195
0.242
0.207
0.005
9
Cyclops insignisMoscow
0.216
0.223
0.204
0.203
0.245
0.201
0.037
0.032
10
C. insignis Borok
0.209
0.216
0.198
0.197
0.241
0.197
0.027
0.023
0.014
11
Diacyclops bicuspidatus
0.186
0.200
0.183
0.180
0.224
0.194
0.090
0.085
0.098
0.087
12
Thermocyclops oithonoides
0.220
0.203
0.191
0.186
0.216
0.184
0.128
0.123
0.120
0.116
0.098
13
Thermocyclops crassus
0.218
0.203
0.181
0.178
0.220
0.204
0.127
0.125
0.127
0.125
0.093
0.084
14
Mesocyclops leuckarti
0.198
0.183
0.174
0.171
0.215
0.181
0.133
0.128
0.123
0.122
0.090
0.087
0.079
15
Macrocyclops albidus
0.221
0.198
0.197
0.191
0.215
0.178
0.145
0.143
0.137
0.134
0.099
0.131
0.114
0.105
16
Macrocyclops distinctus
0.230
0.216
0.192
0.192
0.218
0.191
0.152
0.148
0.143
0.142
0.113
0.127
0.143
0.131
0.123
17
Megacyclops viridis Borok1
0.213
0.207
0.195
0.192
0.230
0.191
0.137
0.134
0.134
0.123
0.105
0.116
0.117
0.114
0.130
0.139
18
M. viridis Borok2
0.207
0.210
0.206
0.203
0.236
0.198
0.123
0.122
0.127
0.119
0.085
0.114
0.110
0.105
0.125
0.130
0.050
All positions containing gaps and missing data were eliminated (complete deletion option). Final dataset, 657 positions.Table 3
Pairwise uncorrected genetic distance among 12 Cyclopidae sequences based on comparison of concatenated ITS1/ITS2 rRNA sequences.
#
Species name
1
2
3
4
5
6
7
8
9
10
11
1
Cyclops kolensis
2
Cyclops strenuus
0.006
3
Cyclops insignis Moscow
0.054
0.048
4
C. insignis Borok
0.048
0.042
0.006
5
Diacyclops bicuspidatus
0.250
0.244
0.232
0.232
6
Thermocyclops oithonoides
0.333
0.330
0.330
0.330
0.313
7
Thermocyclops crassus
0.301
0.298
0.295
0.295
0.262
0.188
8
Mesocyclops leuckarti
0.324
0.321
0.310
0.310
0.286
0.280
0.241
9
Macrocyclops distinctus
*
0.414
0.411
0.408
0.405
0.393
0.429
0.408
0.369*
10
Macrocyclops albidus
0.393
0.390
0.381
0.378
0.357
0.369
0.366
0.336
0.321
11
Megacyclops viridis Borok1*
0.381
0.384
0.396
0.390
0.387
0.387
0.387
0.390*
0.455
0.438
12
M. viridis Borok2*
0.375
0.378
0.384
0.378
0.363
0.366
0.357
0.381*
0.438
0.429
0.149
All positions containing gaps and missing data were eliminated (complete deletion option). Final datasets, 336 positions for concatenated ITS1/ITS2. The names of the three taxa most evolutionarily distant fromMesocyclops and the values of the corresponding genetic distances are shown with asterisk.The cladogram based on comparison of the 28S rDNA sequences reflected the evolutionary history of the analyzed species (Figure2).Oncaea sp. was used as the out group. Similar topologies and levels of support at most nodes were obtained for all 28S phylogenetic trees constructed using the ML and ME methods. The specimens belonging to the order Cyclopoida with high bootstrap support (ML/ME 99) formed two major clades on the tree (Figure 2). One clade combined the marine cyclopoid copepods species, whereas the freshwater species specimens formed the second clade. Thep-distance between these two clades varied in the range of 0.171–0.245 (Table 2). The 28S phylogenetic tree revealed detailed relationships among theOithona spp. with high bootstrap support (ML 79, 100 and ME 76, 100). However,P. nana (Cyclopettidae family) andOithona spp. (Oithonidae family) were poorly resolved. Notably, this study is the second on the molecular phylogenetics of theOithona spp.; the previous study described the phylogenetic relationships between threeOithona spp.:O. similis,O. atlantica, andO. nana [22].Figure 2
Phylogenetic relationships of Cyclopoida based on ~700 bp of the 28S rRNA gene. The consensus cladogram inferred from the 28S ribosomal DNA fragment sequence data of 16 Podoplea superorder species using maximum likelihood (ML) analysis under the HKY+G model and minimum evolution (ME) analysis. The numbers above branches indicate bootstrap percentages. The values are listed for ML/ME.The cladogram based on the comparison of the concatenated ITS1/ITS2 sequences is shown in Figure3. Notably, the 28S and ITS1/ITS2 cladograms had several common features, reflecting the evolutionary history of the analyzed freshwater cyclopoid copepods species. Both cladograms revealed thatD. bicuspidatus and specimens of theCyclops genus with high bootstrap values (>80) are separated from other studied freshwater copepods in a distinct clade. Thep-distance betweenD. bicuspidatus andCyclops spp. calculated based on ITS1/ITS2 analysis varied in the range of 0.232–0.250, whereas thep-distance betweenD. bicuspidatus andThermocyclops spp. varied in the range of 0.298–0.333, and thep-distance betweenD. bicuspidatus and other analyzed freshwater species varied in the range of 0.310–0.405. This result is consistent with a previous phylogenetic study based on 18S rDNA sequence analysis [48]. Notably, the systematic position of this species, based solely on the analysis of morphological characteristics, remained unclear.Diacyclops bicuspidatus is considered to be evolutionarily closer toThermocyclops spp. [8].Figure 3
Phylogenetic relationships of Cyclopoida based on ~500 bp of concatenated ITS1/ITS2 rDNA sequences. The consensus cladogram inferred from the ITS1-ITS2 ribosomal DNA fragment sequence data of 10 species of the Cyclopidae family using maximum likelihood analysis under the HKY+G model and minimum evolution (ME) analysis. The numbers above branches indicate bootstrap percentages. The values are listed for ML/ME.Another important conclusion from the analysis of 28S and ITS1/ITS2 cladograms relates to the systematic position ofC. strenuus. TheCyclops genera subclade was divided into theC. kolensis:C. strenuus subsubclade and theC. insignis subsubclade (Figures 2 and 3). Thep-distance betweenC. strenuus andC. kolensis calculated based on ITS1/ITS2 analysis is 0.006, whereas thep-distance betweenC. strenuus andC. insignisvaried in the range of 0.042–0.054. Therefore,C. strenuus is more closely related toC. kolensis than toC. insignis. Notably, the phylogenetic relationships between the studiedCyclops species could not be elucidated solely on the basis of morphological characteristics.The only difference between the 28S and ITS1/ITS2 cladograms within freshwater copepods was the position ofM. leuckarti (Figures 2 and 3). The cladogram based on comparison of the 28S rDNA sequences showed that theM. leuckarti andThermocyclops cluster together to form a separate subclade (Figure 2). This result is consistent with the previous observation that theMesocyclops andThermocyclops genera are phylogenetically closely related, which was confirmed by the similarity of morphological characteristics and using molecular data [42]. ITS1/ITS2 analysis revealed thatM. leuckarti is located separately fromThermocyclops and other clades (Figure 3). Using phylogenetic networks, we analyzed whether theM. leuckarti position in the ITS1/ITS2 cladogram was caused by different contributions of ITS1 and ITS2 sequences to the phylogenetic signal.


### 3.3. Phylogenetic Networks

A reticulogram-based phylogenetic network inference approach was used to verify the reticulate evolution of the studied copepods. Concatenated ITS1/ITS2 sequences of 10 species from the Cyclopidae family were used for reticulogram reconstruction. The reticulogram revealed a network withMesocyclops andThermocyclops clustered together and a reticulation (lateral branch) connectingM. leuckarti to theMacrocyclops clade node (Figure 4). Therefore, the reticulogram indicated the reticulation inMesocyclops evolution.Figure 4
Reticulogram for concatenated ITS1/ITS2 sequences of 10 species of the Cyclopidae family. The red dashed line indicates the reticulation event connectingM. leuckarti to theMacrocyclops clade node. The number of internal vertices begins with n+1, where n is the number of leaves. The order of internal vertices distribution corresponds to the increasing lengths of the 22 reticulogram edges.A split network represents incompatible edges of trees as a band of parallel edges. Parallel edges split a network into two sets of nodes. Split-graph analysis of concatenated ITS1/ITS2 sequences of 10 species from the Cyclopidae family revealed a reticulate relationship betweenMesocyclops,Thermocyclops, andMacrocyclops with high reliability (Figure 5). All principal splits were well supported. Two splits were observed in the ITS1/ITS2 split network. The first split (parallel edges highlighted with bold red) separatedM. leuckarti,T. oithonoides, andT. crassus with 80.2% bootstrap support. The second split (parallel edges highlighted with bold blue) separatedM. leuckarti,M. distinctus, andM. albidus with 65.0% bootstrap support.Figure 5
Split networks for concatenated ITS1/ITS2 sequences of 10 species of the Cyclopidae family. Split network based on concatenated ITS1/ITS2 sequences; the split separatingM. leuckarti,T. oithonoides, andT. crassus is indicated in bold red; the split separatingM. leuckarti,M. distinctus, andM. albidus is indicated in bold blue; purple indicates theM. leuckartireticulate relationship withThermocyclops andMacrocyclops. The values on the branches indicate bootstrap percentages.In addition to the network data, we performed phylogenetic reconstruction based on independent ITS1 and ITS2 analyses using probabilistic and distance methods. Irrespective of the method used, the main difference between the topologies of the ITS1 and ITS2 phylogenetic trees was as follows: based on the ITS1 analysis,M. leuckarti is clustered withThermocyclops, whereas the ITS2 analysis revealed thatM. leuckarti clustered withMacrocyclops (Figures 6(a)–6(d)).Phylogenetic relationships of Cyclopidae based on ITS1 and ITS2 sequences. The ITS1 consensus Clustal tree of ten Cyclopidae species constructed using (a) maximum likelihood and minimum evolution and (b) Bayesian inference. The ITS2 consensus Clustal tree of ten Cyclopidae species constructed using (c) maximum likelihood and minimum evolution and (d) Bayesian inference. The consensus Gblocks-treated MAFFT trees of eight Cyclopidae species constructed using maximum likelihood and minimum evolution: (e) ITS1, (f) ITS2. The numbers above the branches indicate bootstrap percentages and Bayesian posterior probabilities. The values are listed for ML/ME. Missing or weakly supported nodes (<50% or 0.5) are denoted by a ‘‘—’’.
(a)
(b)
(c)
(d)
(e)
(f)The impact of the chosen DNA sequence on the clustering ofM. leuckarti might reflect the different evolutionary histories of ITS1 and ITS2, which indicates the potential hybrid origin ofM. leuckarti. However, the values of bootstrap support for the clustering ofMesocyclops andThermocyclops and ofMesocyclops andMacrocyclops depended on the method used for phylogenetic tree reconstruction and varied over a wide range (Figures 6(a)–6(d)).Phylogenetic trees can be inconsistent due to the so-called long-branch attraction (LBA) phenomenon, which occurs when two nonadjacent taxa share many homoplastic character states along long branches and/or from uncorrected sequence alignments. Interpretation of the observed similarity depends on the method used for phylogenetic analysis, and this similarity can often be interpreted as homology. Model-based methods are most resistant to LBA, but these methods can exhibit LBA if their assumptions are seriously violated or if there are insufficient taxa in the analysis to accurately estimate the parameters of the evolutionary model [126]. Taxon sampling is a crucial factor for avoiding LBA in phylogenetic analysis [127]. The inclusion of additional taxa in phylogenetic analysis increases the accuracy of the inferred topology by dispersing homoplasty across the tree and reducing the effect of LBA. The LBA effect might also be revealed by exclusion of the long-branched taxon from the analysis [127].To reduce the possible effects of LBA and correctness of the sequence alignment, we used the following approaches: (1) three taxa the most evolutionarily distant fromM. leuckarti (M. distinctus,M. viridis Borok1, andM. viridis Borok2) were removed from the list of species used for ITS1 and ITS2 phylogenetic tree reconstruction (the taxa selection was based on the data presented in Table 3) and (2) ITS1 and ITS2 sequences were aligned using new multiple sequence alignment programs to eliminate poorly aligned and highly divergent regions (see Section 2). The final ITS1 and ITS2 alignments are shown in Supplementary Figures 2(a) and 2(b), and the resulting phylogenetic trees are shown in Figures 6(e) and 6(f). Based on the ITS1 analysis,M. leuckarti clustered withThermocyclops, whereas the ITS2 analysis revealed thatM. leuckarti clustered withMacrocyclops (Figures 6(e) and 6(f)). Notably, the topology of the new phylogenetic trees had high bootstrap support: 84 (ML)/77 (ME) for ITS1 and 77 (ML)/72 (ME) for ITS2 (Figures 6(e) and 6(f)).We think that one of the most intriguing explanations for the observed differences in the clustering ofM. leuckarti is the interspecific hybridization between extinct taxa (presumably closely related) that were ancestral to bothMesocyclops,Macrocyclops, andThermocyclops. However, a rigorous proof of this hypothesis requires further analysis of a larger number of species. This will be the subject of our further research.


## 3.1. Characteristics and Analysis of an rDNA Sequence Dataset

In each species, the nucleotide sequences of the amplified rDNA region and five clones were not significantly different from each other. The frequency of variable nucleotides did not exceed the average rate of nucleotide substitutions caused by DNA polymerase errors, which is approximately one substitution per 1,000 nucleotides. The compared sequences contained both relatively evolutionarily conserved (fragments of 18S and 28S rDNA and the complete 5.8S rDNA) and evolutionary variable genomic regions (ITS1 and ITS2). For different taxa, the ITS1 and ITS2 sequences vary significantly among individuals at the inter- and intrapopulation levels; furthermore, these sequences can exhibit intragenomic variability [25, 41, 53, 54]. Recently, a high level of intrapopulation polymorphism of the 28S rDNA sequences was observed withinOithona spp. [22]. However, there are instances of strong evolutionary conservation of the 28S and ITS sequences [15, 23, 30, 34]. Notably, theM. leuckarti ITS2 sequence obtained in this study did not exhibit any nucleotide substitutions compared to theM. leuckarti ITS2 sequence described previously (GenBank accession number GQ848499) [42]. Therefore, the strong evolutionary conservation of ITS1 and ITS2 sequences is a characteristic feature of the copepod species analyzed in this study.In this study, the applicability of different segments of rDNA containing the ITS1 and ITS2 regions, the 5.8S RNA gene, and fragments of the 18S and 28S rRNA genes was examined for reconstruction of the phylogenetic relationships among freshwater cyclopoid copepods. The 5.8S gene and the analyzed fragment of the 18S gene were not considered for phylogenetic reconstruction due to their short length and strong evolutionary conservation: only a few nucleotide substitutions were detected by comparing these sequences with evolutionary distant species (data not shown).For 15 specimens of Cyclopoida species (including the two marine species), the average length of the 28S gene fragment sequenced was 1051 bp. We trimmed these sequences to 703 bp and compared them with the 28S gene sequences of marine species available in GenBank. These 703 bp of 28S rDNA sequences were aligned, and 342 variable sites were observed.Oncaea sp. (Oncaeidae family) was used as the out group.The ITS1 sequence lengths varied from 267 to 388 bp among the 13 Cyclopidae specimens. The ITS1 sequence alignments possessed 442 characters, and among them, 283 were variable. The ITS2 sequence lengths varied from 188 to 262 bp among the 13 Cyclopidae specimens. ITS2 sequence alignment possessed 302 characters, and 190 were variable.All alignment sets were examined for homogeneity of base frequencies and substitution saturation. The average base frequencies of the 28S gene fragment (A=20.52, C=25.33, G=32.75, and T=21.40%) differed from the ITS1 (A=14.27, C=30.68, G=27.55, and T=27.50%) and ITS2 (A=13.08, C=29.64, G=30.08, and T=27.19%) regions. Gaps were excluded while estimating the average base frequencies of the ITS sequences. Using the chi-squared test, no significant differences were observed in the base compositions of the 28S (χ2 = 39.77, df=54, and P=0.93), ITS1 (χ2=22.04, df=36, and P=0.97), and ITS2 (χ2=22.65, df=36, and P=0.96) sequences among different taxa.To analyze whether the divergence of 28S, ITS1, and ITS2 rDNA fragments among species was saturated, we performed a substitution saturation test and generated saturation plots. Using DAMBE, the substitution saturation test revealed an Iss value that was significantly(P<0.001) lower than the Iss.c in all cases (Table 1). This result indicated the suitability of the data for phylogenetic analysis. The total numbers of transition and transversion substitutions were plotted individually against model-corrected maximum-likelihood pairwise distances for the 28S, ITS1, and ITS2 sequences (see Supplementary Figure 1 in Supplementary Material available online at http://dx.doi.org/10.1155/2014/926342). Using linear regression analysis on the 28S, ITS1, and ITS2 saturation graphs, the coefficients of determination (R2) were calculated for both classes of substitutions: for transitions, the R2 values were 0.79, 0.68, and 0.74 for the 28S, ITS1, and ITS2 sequences, respectively; for transversions, the R2 values were 0.95, 0.93, and 0.91 for the 28S, ITS1, and ITS2 sequences, respectively. The R2 values indicated that no less than 70% of the total variation in pairwise transitions and transversions could be explained by the linear relationship between pairwise distances and the total number of transitions and transversions. All saturation plots showed significant linear correlations (Supplementary Figure 1). Therefore, both transitions and transversions steadily accumulated as the corrected pairwise divergence increased, indicating that saturation was not reached.Table 1
False test of substitution saturation.
Alignment
Iss
Iss.c
Std. error
28S
0.200
0.739
0.019
ITS1
0.427
0.679
0.044
ITS2
0.376
0.665
0.043
Testing whether the observed Iss is significantly(P<0.001, two-tailed t-test) lower than the Iss.c for a symmetrical tree.


## 3.2. Distance Analyses and Phylogenetic Tree Reconstruction

Phylogenetic analysis of the cyclopoid copepods species based on rDNA showed that the 28S rDNA sequences are informative for the phylogeny of both higher-level and closely related Copepoda species, whereas the ITS1 and ITS2 sequences are highly informative for reconstruction of the evolutionary history of closely related species. The ITS1 and ITS2 sequences are known to evolve more rapidly than the ribosomal RNA genes. Consistent with this observation, in this study, the pairwise ITS1/ITS2p-distances were significantly higher thanthe 28Sp-distances (compare Tables 2 and 3). These data are consistent with other studies showing considerable variation in ITS1 and ITS2 divergence levels among different groups of copepods [40, 42, 52, 124, 125]. In this study, fragments of 28S rDNA sequences were used for the analysis of marine and freshwater cyclopoid copepods species, whereas ITS1 and ITS2 sequences were used exclusively for the analysis of freshwater cyclopoid copepods species.Table 2
Pairwise uncorrected genetic distance among 18 sequences of Cyclopoida based on comparison of 28S rRNA gene.
#
Species name
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
1
Oncaeasp.
2
Paracyclopina nana
0.229
3
Oithona similis
0.230
0.148
4
Oithona helgolandica
0.227
0.143
0.011
5
Oithona nana
0.233
0.181
0.155
0.146
6
Oithona simplex
0.244
0.152
0.175
0.166
0.177
7
Cyclops kolensis
0.209
0.223
0.197
0.195
0.244
0.212
8
Cyclops strenuus
0.207
0.220
0.197
0.195
0.242
0.207
0.005
9
Cyclops insignisMoscow
0.216
0.223
0.204
0.203
0.245
0.201
0.037
0.032
10
C. insignis Borok
0.209
0.216
0.198
0.197
0.241
0.197
0.027
0.023
0.014
11
Diacyclops bicuspidatus
0.186
0.200
0.183
0.180
0.224
0.194
0.090
0.085
0.098
0.087
12
Thermocyclops oithonoides
0.220
0.203
0.191
0.186
0.216
0.184
0.128
0.123
0.120
0.116
0.098
13
Thermocyclops crassus
0.218
0.203
0.181
0.178
0.220
0.204
0.127
0.125
0.127
0.125
0.093
0.084
14
Mesocyclops leuckarti
0.198
0.183
0.174
0.171
0.215
0.181
0.133
0.128
0.123
0.122
0.090
0.087
0.079
15
Macrocyclops albidus
0.221
0.198
0.197
0.191
0.215
0.178
0.145
0.143
0.137
0.134
0.099
0.131
0.114
0.105
16
Macrocyclops distinctus
0.230
0.216
0.192
0.192
0.218
0.191
0.152
0.148
0.143
0.142
0.113
0.127
0.143
0.131
0.123
17
Megacyclops viridis Borok1
0.213
0.207
0.195
0.192
0.230
0.191
0.137
0.134
0.134
0.123
0.105
0.116
0.117
0.114
0.130
0.139
18
M. viridis Borok2
0.207
0.210
0.206
0.203
0.236
0.198
0.123
0.122
0.127
0.119
0.085
0.114
0.110
0.105
0.125
0.130
0.050
All positions containing gaps and missing data were eliminated (complete deletion option). Final dataset, 657 positions.Table 3
Pairwise uncorrected genetic distance among 12 Cyclopidae sequences based on comparison of concatenated ITS1/ITS2 rRNA sequences.
#
Species name
1
2
3
4
5
6
7
8
9
10
11
1
Cyclops kolensis
2
Cyclops strenuus
0.006
3
Cyclops insignis Moscow
0.054
0.048
4
C. insignis Borok
0.048
0.042
0.006
5
Diacyclops bicuspidatus
0.250
0.244
0.232
0.232
6
Thermocyclops oithonoides
0.333
0.330
0.330
0.330
0.313
7
Thermocyclops crassus
0.301
0.298
0.295
0.295
0.262
0.188
8
Mesocyclops leuckarti
0.324
0.321
0.310
0.310
0.286
0.280
0.241
9
Macrocyclops distinctus
*
0.414
0.411
0.408
0.405
0.393
0.429
0.408
0.369*
10
Macrocyclops albidus
0.393
0.390
0.381
0.378
0.357
0.369
0.366
0.336
0.321
11
Megacyclops viridis Borok1*
0.381
0.384
0.396
0.390
0.387
0.387
0.387
0.390*
0.455
0.438
12
M. viridis Borok2*
0.375
0.378
0.384
0.378
0.363
0.366
0.357
0.381*
0.438
0.429
0.149
All positions containing gaps and missing data were eliminated (complete deletion option). Final datasets, 336 positions for concatenated ITS1/ITS2. The names of the three taxa most evolutionarily distant fromMesocyclops and the values of the corresponding genetic distances are shown with asterisk.The cladogram based on comparison of the 28S rDNA sequences reflected the evolutionary history of the analyzed species (Figure2).Oncaea sp. was used as the out group. Similar topologies and levels of support at most nodes were obtained for all 28S phylogenetic trees constructed using the ML and ME methods. The specimens belonging to the order Cyclopoida with high bootstrap support (ML/ME 99) formed two major clades on the tree (Figure 2). One clade combined the marine cyclopoid copepods species, whereas the freshwater species specimens formed the second clade. Thep-distance between these two clades varied in the range of 0.171–0.245 (Table 2). The 28S phylogenetic tree revealed detailed relationships among theOithona spp. with high bootstrap support (ML 79, 100 and ME 76, 100). However,P. nana (Cyclopettidae family) andOithona spp. (Oithonidae family) were poorly resolved. Notably, this study is the second on the molecular phylogenetics of theOithona spp.; the previous study described the phylogenetic relationships between threeOithona spp.:O. similis,O. atlantica, andO. nana [22].Figure 2
Phylogenetic relationships of Cyclopoida based on ~700 bp of the 28S rRNA gene. The consensus cladogram inferred from the 28S ribosomal DNA fragment sequence data of 16 Podoplea superorder species using maximum likelihood (ML) analysis under the HKY+G model and minimum evolution (ME) analysis. The numbers above branches indicate bootstrap percentages. The values are listed for ML/ME.The cladogram based on the comparison of the concatenated ITS1/ITS2 sequences is shown in Figure3. Notably, the 28S and ITS1/ITS2 cladograms had several common features, reflecting the evolutionary history of the analyzed freshwater cyclopoid copepods species. Both cladograms revealed thatD. bicuspidatus and specimens of theCyclops genus with high bootstrap values (>80) are separated from other studied freshwater copepods in a distinct clade. Thep-distance betweenD. bicuspidatus andCyclops spp. calculated based on ITS1/ITS2 analysis varied in the range of 0.232–0.250, whereas thep-distance betweenD. bicuspidatus andThermocyclops spp. varied in the range of 0.298–0.333, and thep-distance betweenD. bicuspidatus and other analyzed freshwater species varied in the range of 0.310–0.405. This result is consistent with a previous phylogenetic study based on 18S rDNA sequence analysis [48]. Notably, the systematic position of this species, based solely on the analysis of morphological characteristics, remained unclear.Diacyclops bicuspidatus is considered to be evolutionarily closer toThermocyclops spp. [8].Figure 3
Phylogenetic relationships of Cyclopoida based on ~500 bp of concatenated ITS1/ITS2 rDNA sequences. The consensus cladogram inferred from the ITS1-ITS2 ribosomal DNA fragment sequence data of 10 species of the Cyclopidae family using maximum likelihood analysis under the HKY+G model and minimum evolution (ME) analysis. The numbers above branches indicate bootstrap percentages. The values are listed for ML/ME.Another important conclusion from the analysis of 28S and ITS1/ITS2 cladograms relates to the systematic position ofC. strenuus. TheCyclops genera subclade was divided into theC. kolensis:C. strenuus subsubclade and theC. insignis subsubclade (Figures 2 and 3). Thep-distance betweenC. strenuus andC. kolensis calculated based on ITS1/ITS2 analysis is 0.006, whereas thep-distance betweenC. strenuus andC. insignisvaried in the range of 0.042–0.054. Therefore,C. strenuus is more closely related toC. kolensis than toC. insignis. Notably, the phylogenetic relationships between the studiedCyclops species could not be elucidated solely on the basis of morphological characteristics.The only difference between the 28S and ITS1/ITS2 cladograms within freshwater copepods was the position ofM. leuckarti (Figures 2 and 3). The cladogram based on comparison of the 28S rDNA sequences showed that theM. leuckarti andThermocyclops cluster together to form a separate subclade (Figure 2). This result is consistent with the previous observation that theMesocyclops andThermocyclops genera are phylogenetically closely related, which was confirmed by the similarity of morphological characteristics and using molecular data [42]. ITS1/ITS2 analysis revealed thatM. leuckarti is located separately fromThermocyclops and other clades (Figure 3). Using phylogenetic networks, we analyzed whether theM. leuckarti position in the ITS1/ITS2 cladogram was caused by different contributions of ITS1 and ITS2 sequences to the phylogenetic signal.


## 3.3. Phylogenetic Networks

A reticulogram-based phylogenetic network inference approach was used to verify the reticulate evolution of the studied copepods. Concatenated ITS1/ITS2 sequences of 10 species from the Cyclopidae family were used for reticulogram reconstruction. The reticulogram revealed a network withMesocyclops andThermocyclops clustered together and a reticulation (lateral branch) connectingM. leuckarti to theMacrocyclops clade node (Figure 4). Therefore, the reticulogram indicated the reticulation inMesocyclops evolution.Figure 4
Reticulogram for concatenated ITS1/ITS2 sequences of 10 species of the Cyclopidae family. The red dashed line indicates the reticulation event connectingM. leuckarti to theMacrocyclops clade node. The number of internal vertices begins with n+1, where n is the number of leaves. The order of internal vertices distribution corresponds to the increasing lengths of the 22 reticulogram edges.A split network represents incompatible edges of trees as a band of parallel edges. Parallel edges split a network into two sets of nodes. Split-graph analysis of concatenated ITS1/ITS2 sequences of 10 species from the Cyclopidae family revealed a reticulate relationship betweenMesocyclops,Thermocyclops, andMacrocyclops with high reliability (Figure 5). All principal splits were well supported. Two splits were observed in the ITS1/ITS2 split network. The first split (parallel edges highlighted with bold red) separatedM. leuckarti,T. oithonoides, andT. crassus with 80.2% bootstrap support. The second split (parallel edges highlighted with bold blue) separatedM. leuckarti,M. distinctus, andM. albidus with 65.0% bootstrap support.Figure 5
Split networks for concatenated ITS1/ITS2 sequences of 10 species of the Cyclopidae family. Split network based on concatenated ITS1/ITS2 sequences; the split separatingM. leuckarti,T. oithonoides, andT. crassus is indicated in bold red; the split separatingM. leuckarti,M. distinctus, andM. albidus is indicated in bold blue; purple indicates theM. leuckartireticulate relationship withThermocyclops andMacrocyclops. The values on the branches indicate bootstrap percentages.In addition to the network data, we performed phylogenetic reconstruction based on independent ITS1 and ITS2 analyses using probabilistic and distance methods. Irrespective of the method used, the main difference between the topologies of the ITS1 and ITS2 phylogenetic trees was as follows: based on the ITS1 analysis,M. leuckarti is clustered withThermocyclops, whereas the ITS2 analysis revealed thatM. leuckarti clustered withMacrocyclops (Figures 6(a)–6(d)).Phylogenetic relationships of Cyclopidae based on ITS1 and ITS2 sequences. The ITS1 consensus Clustal tree of ten Cyclopidae species constructed using (a) maximum likelihood and minimum evolution and (b) Bayesian inference. The ITS2 consensus Clustal tree of ten Cyclopidae species constructed using (c) maximum likelihood and minimum evolution and (d) Bayesian inference. The consensus Gblocks-treated MAFFT trees of eight Cyclopidae species constructed using maximum likelihood and minimum evolution: (e) ITS1, (f) ITS2. The numbers above the branches indicate bootstrap percentages and Bayesian posterior probabilities. The values are listed for ML/ME. Missing or weakly supported nodes (<50% or 0.5) are denoted by a ‘‘—’’.
(a)
(b)
(c)
(d)
(e)
(f)The impact of the chosen DNA sequence on the clustering ofM. leuckarti might reflect the different evolutionary histories of ITS1 and ITS2, which indicates the potential hybrid origin ofM. leuckarti. However, the values of bootstrap support for the clustering ofMesocyclops andThermocyclops and ofMesocyclops andMacrocyclops depended on the method used for phylogenetic tree reconstruction and varied over a wide range (Figures 6(a)–6(d)).Phylogenetic trees can be inconsistent due to the so-called long-branch attraction (LBA) phenomenon, which occurs when two nonadjacent taxa share many homoplastic character states along long branches and/or from uncorrected sequence alignments. Interpretation of the observed similarity depends on the method used for phylogenetic analysis, and this similarity can often be interpreted as homology. Model-based methods are most resistant to LBA, but these methods can exhibit LBA if their assumptions are seriously violated or if there are insufficient taxa in the analysis to accurately estimate the parameters of the evolutionary model [126]. Taxon sampling is a crucial factor for avoiding LBA in phylogenetic analysis [127]. The inclusion of additional taxa in phylogenetic analysis increases the accuracy of the inferred topology by dispersing homoplasty across the tree and reducing the effect of LBA. The LBA effect might also be revealed by exclusion of the long-branched taxon from the analysis [127].To reduce the possible effects of LBA and correctness of the sequence alignment, we used the following approaches: (1) three taxa the most evolutionarily distant fromM. leuckarti (M. distinctus,M. viridis Borok1, andM. viridis Borok2) were removed from the list of species used for ITS1 and ITS2 phylogenetic tree reconstruction (the taxa selection was based on the data presented in Table 3) and (2) ITS1 and ITS2 sequences were aligned using new multiple sequence alignment programs to eliminate poorly aligned and highly divergent regions (see Section 2). The final ITS1 and ITS2 alignments are shown in Supplementary Figures 2(a) and 2(b), and the resulting phylogenetic trees are shown in Figures 6(e) and 6(f). Based on the ITS1 analysis,M. leuckarti clustered withThermocyclops, whereas the ITS2 analysis revealed thatM. leuckarti clustered withMacrocyclops (Figures 6(e) and 6(f)). Notably, the topology of the new phylogenetic trees had high bootstrap support: 84 (ML)/77 (ME) for ITS1 and 77 (ML)/72 (ME) for ITS2 (Figures 6(e) and 6(f)).We think that one of the most intriguing explanations for the observed differences in the clustering ofM. leuckarti is the interspecific hybridization between extinct taxa (presumably closely related) that were ancestral to bothMesocyclops,Macrocyclops, andThermocyclops. However, a rigorous proof of this hypothesis requires further analysis of a larger number of species. This will be the subject of our further research.


## 4. Conclusion

We evaluated the utility of a ~2000 bp fragment of rDNA (easily amplified by universal primers) for the phylogenetic reconstruction of the relationships of Copepoda species. Our data showed that the 28S rDNA and the ITS1 and ITS2 regions are highly informative for the phylogeny of both higher-level and closely related Copepoda species. Comparative analysis of the ITS1 and ITS2 nucleotide sequences among closely related Copepoda species revealed an unusual evolutionary history of these spacer sequences; therefore, the ITS1 and ITS2 regions might contain different phylogenetic signals.



---

*Source: 926342-2014-08-18.xml*