Morphological and molecular diversity patterns of the genus Tropodiaptomus Kiefer, 1932 (Copepoda, Calanoida, Diaptomidae) in Thailand

Tropodiaptomus is one of the most specious genera in the family Diaptomidae, but it is often rare in terms of distribution and abundance. Moreover, Tropodiaptomus species show a noteworthy variability in some of the morphological characters considered of prime importance in diaptomid taxonomy, and the presence of cryptic or pseudocryptic species is likely. Thus, through a geographically-wide sampling in Thailand, we aimed to investigate the local diversity of the genus and to compare the morphological and molecular diversity pattern based on mitochondrial and nuclear genes sequences. DNA taxonomy was also implemented in order to check whether the Tropodiaptomus lineages were independent species according to the “evolutionary genetic species concept”. Six Tropodiaptomus morphospecies were found, three of which are putative species new to Science pending a formal description. The finding of such a high incidence of undescribed species stresses the existence of a significant “Linnean shortfall” affecting Thai diaptomids. The molecular results showed that most of the studied species could be identified consistently with their morphology-based taxonomy. However, Tropodiaptomus vicinus and T. cf. lanaonus showed a high level of genetic diversity, suggesting that traditional morphological techniques might be inadequate for correctly assessing their taxonomical status.

Sequence alignment and phylogenetic analyses. 12S and 28S chromatograms were checked and manually edited using the software Chromas v.2.6.2 (Technelysium, Pty. Ltd. 1998, Queensland, Australia). For ITS2, the forward and reverse sequences were assembled and the contig was examined through SEQUENCHER v.4.1.4 as suggested by Fontaneto et al. 29 . For ITS2, homozygotes consensus sequences were obtained directly from the contig. Heterozygotes sequences were carefully inspected and manually constructed when the chromatograms exhibited double peaks at a single position. Conversely, the two haplotypes were reconstructed using CHAMPURU v.1.0 30 (http:// jfflot. mnhn. fr/ champ uru/) for those sequences that showed multiple double peaks at several positions. All sequences were aligned using the Clustal W method 31 as implemented in the software MEGAX 32 .
In order to test whether the mitochondrial and nuclear fragments could be combined for joint analyses, the incongruence length difference test (ILD) 33 as implemented in PAUP* v. 4.0b10 34 was used. According to Cunningham 35 , if p > 0.01, pooling the data improves the phylogenetic accuracy and thus it is admissible to merge the tested datasets into a single matrix. This condition was fulfilled both for the concatenated 12S and ITS2 sequences (p = 0.63), and for the concatenation of all the genetic markers analyzed in the frame of this study (i.e., 12S, ITS2 and 28S; p = 1). Accordingly, two combined mito-nuclear datasets (i.e., the "12S-ITS2 dataset" and the "12S-ITS2-28S dataset") were analyzed in the frame of this work. The first dataset includes all the available novel 12S and ITS2 sequences; those individuals which were heterozygous at the ITS2 locus are here represented by two concatenated sequences indicated by different letters, bearing the alternative ITS2 haplotype. The second dataset includes a subset of the 12S and ITS2 sequences and the 28S sequences. For both combined datasets the software packages MrBayes v. 3.2.6 36 and PhyML v. 3 37 were used for inferring phylogenetic relationships through Bayesian inference of phylogeny (BI) and maximum likelihood (ML) analysis. As support measures for the nodes, bootstrap values (BS) 38 were calculated with 1000 replicates in the ML trees, whereas in the BI tree the posterior probability values (PP) reported. The best evolution model for each dataset was selected in the software Parti-tionFinder ver. 1.0.1 39 under the "Akaike Information Criterion" (AIC) 40 . The General Time-Reversible model of sequence evolution with a proportion of invariable sites was used for the 28S dataset (GTR + I). Instead, a GTR with a proportion of invariable sites and gamma-distributed rate variation among sites (GTR + I + Γ; nst = 6) was selected as the best evolutionary model for both the 12S and ITS2 datasets. In the BI analyses, two independent Markov chain Monte Carlo analyses were carried out for 1,000,000 generations (temp.: 0.2; default priors) with www.nature.com/scientificreports/ Species delimitation methods. In the frame of this paper, we followed the 'unified species concept' described by De Queiroz 42 . Accordingly, morphological similarity or identity was not considered per se sufficient evidence for conspecificity, and single-locus DNA taxonomy approaches were implemented to explore the possible presence of groups of putative species rank within each studied morphospecies 42,43 . We applied two independent methods of species delimitation based on the 12S and ITS2 datasets. The "assemble species by automatic partitioning" (ASAP) 44 method was implemented using the online ASAP server (https:// bioin fo. mnhn. fr/ abi/ public/ asap) with the following settings: fixed seed value = − 1 (i.e., no fixed seed value was used), and simple distance. The "multiple rate Poisson Tree Processes" (mPTP) 45 model was run through the online mPTP server (https:// mptp.h-its. org/). ASAP and mPTP analyses were performed on Tropodiaptomus spp. sequences, with the exclusion of the outgroup. Genetic distances were calculated within and between the Tropodiaptomus clades of putative species-level using the Kimura two-parameter model with pairwise deletion in MEGAX 46 . Ethics statement. The present study was approved by the ethics committee of Kasetsart University (approval no. ACKU61-SCI-004) for collecting the Tropodiaptomus specimens.

Results
Morphological analysis. In the frame of this study, six Tropodiaptomus morphospecies were identified in the 23 sampled locations in Thailand. These are Tropodiaptomus oryzanus, T. vicinus, T. cf. lanaonus, Tropodiaptomus sp.1, Tropodiaptomus sp.2, and Tropodiaptomus sp.3 (Table 1). Studied specimens mostly showed variation in i) the ornamentation of the basis and second exopod segment of adult male right P5, ii) the inner margin of the exopod in adult male left P5, and iii) the length of the process on the antepenultimate segment of adult male right antennule (see Table 2). Moreover, Lai et al. 16 also reported that in the inner margin of adult male left P5 exopod in T. vicinus two or three lobes might be present; however, Thai specimens consistently showed the presence of two lobes. Conversely, our T. vicinus specimens showed variability in the ornamentation on second exopod of adult male right P5 and length of spinous process on the antepenultimate segment of adult male right antennule. www.nature.com/scientificreports/ The morphology of the specimens here ascribed to Tropodiaptomus cf. lanaonus do not agree with the original description of the species 22 in two characters: (i) the length of the spinous process in the antepenultimate segment of adult male right antennule is longer than segment 21 in the original description, whereas in our specimens (from sites KPK1, NT, NER and NMK4, see Table 1) it is 1/2 to 3/4 of segment 21, and (ii) the ornamentation on the basis of adult male right P5 has one apophysis and one hyaline lamella in the original description, which is in accordance with Thai specimens from KPK1 and NT, whereas studied specimens from NER and NMK4 have only the hyaline lamella and no apophysis. In addition, Thai T. cf. lanaonus also showed variability in the morphology of the basis and second exopod of adult male right P5 (Table 2). Moreover, specimens from NER have a group of spinules near inner margin lobe on the exopod of adult male left P5, this character does not show in specimens from KPK, NMK4 and NT (see Fig. 2).    www.nature.com/scientificreports/ We refrained from ascribing to any known species the Thai Tropodiaptomus populations from DKT (here referred to as Tropodiaptomus sp.1) and NP2 and NP3 (here referred to as Tropodiaptomus sp.2) because the morphology of adult male P5 and right antennule in Tropodiaptomus sp.1 and Tropodiaptomus sp.2 are different from all the Tropodiaptomus species known to date. Moreover, Tropodiaptomus sp.1 and Tropodiaptomus sp.2 differ one from the other for the following characteristics: the basis of the right P5 of Tropodiaptomus sp.1 has two processes and one hyaline lamella but it has one process and one hyaline lamella in Tropodiaptomus sp.2; the length of the spinous process in antepenultimate segment of adult male right antennule in Tropodiaptomus sp.1 is 3/4 of or equal to segment 21, while it is longer than segment 21 in Tropodiaptomus sp.2; and the shape of second exopod of the right P5 is trapezoidal in Tropodiaptomus sp.1 and rectangular in Tropodiaptomus sp.2. Based on all these characters, we suggest that they are two new putative species pending a formal description.
In addition, a single female whose morphology was not ascribable to any known Tropodiaptomus species was collected from the site KSM1 and here reported as Tropodiaptomus sp.3. Unfortunately, no Tropodiaptomus males were collected from this site.
The most common morphospecies collected in the frame of present study was Tropodiaptomus vicinus, which was found in 13 localities; it was found co-occurring with T. cf. lanaonus (site KPK1, 13 October 2017) and T. oryzanus (sites KDP1 and SNG4, 3 June 2019) ( Table 1). All the other Tropodiaptomus taxa and populations were found with no co-occurring congeneric species in each locality in this study.

Molecular analyses.
Overall, a total of 62 Tropodiaptomus individuals were molecularly analyzed. The mitochondrial 12S and the nuclear ITS2 markers were successfully amplified in all the 62 individuals; conversely, 28S sequences were produced for a subset of 18 individuals only. See Table 1 for a synopsis of the fragments amplified for each studied specimen and their GenBank Accession Numbers. Unfortunately, no 28S amplicons could be obtained for representatives of clades VI, VIII, and XII (see Fig. 3), which are thus not represented in the analyses which include this fragment. In addition, 12S, ITS2 and 28S sequences were obtained from the Italian Eudiaptomus intermedius specimen used as outgroup (ANs: OL584216 for 12S, OL584119 for 28S, OL630143 for ITS2).
The phylogenetic trees (BI and ML analysis) based on the aligned 899 bp-long fragment of the concatenated "12S-ITS2 dataset" (367 bp for 12S, 533 bp for ITS2; Fig. 3) show a congruent topology, with the occurrence of 13 major clades with strong to moderate node support, hereafter indicated with roman numerals (i.e., clades I-XIII, see Fig. 3). The Thai populations ascribed to the morphospecies Tropodiaptomus vicinus resulted to be paraphyletic, since they were split into two major monophyletic groups with no sister-clades relationship and further subdivided in five clades. "Clade I" includes 15 specimens that belong to five different sites located in the northeastern part of Thailand (KDP1, SNG1, SNG2, SNG3 and TP; PP = 1, BS = 98%). Clade II includes two specimens from a single site located in eastern part of Thailand (RBG13; PP = 1; BS = 100%). Clade III includes four specimens from two locations in eastern and southern part of Thailand (SMR1 and TMP1; PP = 0.99, BS = 80%). Clade IV includes eight specimens from five locations in central, northeastern, and western parts (NH1, NN04, TMG1 and KPK1; PP = 1, BS = 100%). Clade V includes three specimens from one location in the northeastern part of Thailand (KDB; PP = 1, BS = 100%). All the other morphospecies proved to constitute monophyletic groups, although four well-characterized subclades are present within Tropodiaptomus cf. lanaonus (clades VI, VII, VIII, and IX, see Fig. 3).
The BI and ML phylogenetic trees based on the aligned 1709 bp-long 12S-ITS2-28S dataset (367 bp for 12S, 533 bp for ITS2, 809 bp for 28S; Fig. 4A) group the concatenated sequences in 10 major clades with moderate to strong statistical support (Fig. 4A). These ten clades are consistent with the thirteen clades singled out based on the 12S-ITS2 dataset (Fig. 3), with the exception of clades VI, VIII and XII, which are missing from this analysis due to the failure of 28S amplification (see above).
The Tropodiaptomus spp. haplotypes network based on the 28S nuclear marker shows the occurrence of seven haplotypes (Fig. 4B). The nine available sequences of T. vicinus belong to three different haplotypes. Tropodiaptomus vicinus specimens belonging to clades II, II and V shared the same haplotype with T. cf. lanaonus (clades VII, IX). One single haplotype was observed for T. oryzanus and Tropodiaptomus sp.1, respectively, whereas two different haplotypes were found within Tropodiaptomus sp.2.

Species delimitation.
Based on a 365 bp-long fragment of the mitochondrial 12S gene, ASAP analysis suggested the existence of 12 groups of putative species rank within the ingroup, splitting both Tropodiaptomus vicinus and T. cf. lanaonus in four groups of putative species rank each; conversely mPTP analysis suggested the existence of 13 groups, splitting Tropodiaptomus vicinus and T. cf. lanaonus in five and four groups of putative species rank, respectively (Supplementary Figs. S1, S2, S3). The known geographical distribution of the clades found within T. vicinus and T. cf. lanaonus are reported in Fig. 5. Results of DNA-based classification based on the 12S gene matched morphology rather well, except for the case where putative cryptic species were found, i.e., within T. vicinus (clades I-V) and T. cf. lanaonus (clades VI-IX).
Based on a 533 bp-long fragment of ITS2, analyses suggested the existence of five (ASAP) or four (mPTP) groups of putative species rank, lumping in the same group specimens and populations characterized by different morphologies (Supplementary Figs. S4, S5, S6 47 , Indonesia (6 species) 48 , Philippines (4 species) 49 , Malaysia (3 species) 15 , and Vietnam (3 species) 17 . Most species found in present and previous studies seem to prefer habitats characterized by a temporary hydroperiod. However, inhabited habitats might differ. For example, T. oryzanus was recorded from rice fields, swamps and manmade ponds in the present study but it was found also in temporary canals and ponds 50 , and roadside canals 51,52 . The co-occurrence of different Tropodiaptomus species has been seldom observed in Thailand, and it always implied the co-occurrence of T. vicinus with other species, i.e., Tropodiaptomus cf. lanaonus, T. megahyaline and T. oryzanus (Table 1; see also: 13,14,51 Kiefer, 1937;Neodiaptomus Kiefer, 1932 andPhyllodiaptomus Kiefer, 1936, with multi-species diaptomid coexistences ranging from two to eight species per site 13,[52][53][54] .
Out of the six morphospecies identified in the frame of present survey, three taxa are putative species new to Science pending a formal description. The finding of such a high incidence of undescribed species suggests the existence of a significant "Linnean shortfall" 55 affecting Thai diaptomid fauna, as also confirmed by the recent description of a new species of the genus 13 and the report of a further new species pending a formal description by Sanoamuang and Dabseepai 14 . Current knowledge about Thai and Oriental Tropodiaptomus species is thus likely largely incomplete. Such a shortfall prevents from getting an exhaustive picture of the morphological and genetic diversity patterns of the genus, and of their taxonomical value.
Comparison between morphological and molecular diversity patterns. The existence of a noteworthy morphological variability within the Philippine and Indian species belonging to the genus Tropodiaptomus was reported by Lai et al. 16 and Ambedkar 18 , respectively. They concluded that both intra-and interspecific variations of adult male P5 can be observed. In the frame of the present study, a certain degree of morphological variability was observed within Tropodiaptomus vicinus and Tropodiaptomus cf. lanaonus only, whereas a negligible variability was observed in the other species ( Table 2). The molecular analyses carried out on the studied species showed that four of the six morphospecies found in the frame of present study, i.e., Tropodiaptomus oryzanus, Tropodiaptomus sp.1, Tropodiaptomus sp.2, and Tropodiaptomus sp.3, could be consistently identified based on their morphology and DNA sequences. Conversely, the morphospecies T. vicinus and T. cf. lanaonus showed a high level of genetic diversity, which suggests that the resolution of traditional morphological techniques may be insufficient for correctly assess their statuses.
The branching patterns of the phylogenetic trees built upon both the "12S-ITS2" and "12S-ITS2-28S" datasets (Figs. 3, 4a) show a noteworthy structuring within these last two morphospecies, which, along with their morphological variability (Table 2) and the paraphyly of T. vicinus (Fig. 3), is suggestive of the possible existence of multiple taxa lumped under these binomia. The two implemented DNA taxonomy analyses based on 12S sequences suggest the occurrence of four (ASAP) and five (mPTP) groups of putative species rank within Tropodiaptomus vicinus, and of four groups of the same rank within Tropodiaptomus cf. lanaonus for both ASAP and mPTP (Fig. 3, Supplementary Figs. S2, S3). Interestingly, such a deep genetic structure corresponds to the high variability in morphological characters observed in these two morphospecies. Conversely, species delimitation analyses based on ITS2 sequences produced results, which are counter-intuitive and not in agreement Table 3. Intra-clade and inter-clades genetic diversity, assessed by Kimura two-parameter distance based on the 12S rRNA-ITS2 dataset. Minimum and maximum values of genetic distance between clades are reported in bold. The presence of n/c in the results denotes cases in which it was not possible to estimate evolutionary distances. Taxa  Clade  n  Within Clade   Between Clade   I  II  III  IV  V  VI  VII  VIII  IX  X  XI  XII (Table 2); however, the observed characters do not allow to consistently distinguish among the clades of putative species rank suggested by DNA taxonomy analyses (Fig. 3, Table 2), and these clades are currently impossible to be told apart based on morphology. However, no clear geographic pattern of the observed molecular diversity pattern is evident (Fig. 5a), so that the existence of a geographical cline of genetic diversity seems not to be supported. Accordingly, the DNA-based clustering of the studied populations in five clades not forming a monophyletic group is neither in agreement with their morphology nor geographical distribution.  www.nature.com/scientificreports/ Tropodiaptomus lanaonus was described as an endemic species from Lake Lanao in Philippines by Kiefer 22 . Since then, the species has not been recorded until the twenty-first century, when Sanoamuang 23 recorded the occurrence of this species in central and northeastern Thailand, and from Laos. Lopez et al. 48 reported that this species has not been found from samples collected in the Philippines during 2008-2015 and hypothesize that the species could be locally extinct. Therefore, until now, this species has been reported only from its type locality (Lake Lanao, Philippines), from Laos, and from the eastern and northeastern parts of Thailand 14 . In the frame of present work, a Tropodiaptomus species close to T. lanaonus was collected also in southern Thailand, but due to its morphological peculiarities, which are suggestive of a possible species-level differentiation between the studied Thai populations and T. lanaonus s.s., we here conservatively ascribed these populations to Tropodiaptomus cf. lanaonus. A morphological re-analysis of the populations reported as Tropodiaptomus lanaonus in recent Thai literature is advisable in order to compare them with the Tropodiaptomus cf. lanaonus sampled in the frame of this survey. Based on molecular data, Tropodiaptomus cf. lanaonus consists of four clades of putative species rank (Fig. 3) with an allopatric distribution ( Fig. 5b) but, as also observed for Tropodiaptomus vicinus (see above), morphological characters do not allow to consistently distinguish among clades.
The spreading of molecular techniques and DNA taxonomy revealed that cryptic or overlooked species are an evolutionary constant also among diaptomids (e.g. [3][4][5][6][7]57 ). The possible occurrence of cryptic, or simply overlooked, species within the widespread, but locally rare, Tropodiaptomus vicinus and T. cf. lanaonus is thus verisimilar and somehow expected. DNA taxonomy methods are powerful tools for the recognition of cryptic species and, according to the International Commission on Zoological Nomenclature (ICZN), "new species can be described on the basis of DNA sequences" 58 (ICZN, 1999; accessed 27 July 2020). However, as stressed by Fontaneto et al. 29 and Dellicour and Flot 43 , the clades of putative species rank singled out by DNA taxonomy methods must be considered just as "primary species hypotheses", and the implementation of an integrative approach including the search for a consensus among the evidence provided by different data sources, e.g., morphology, ecology, biogeography, and ethology, must be adopted before formally proposing their actual species status.
The small inter-clades distances observed for some of the clades of putative species-level found in the frame of present survey (Table 3) suggest caution in order to avoid the risk of an unsupported oversplitting of the www.nature.com/scientificreports/ taxa. In particular, the possible incomplete sampling of the actual genetic diversity of the Thai populations of Tropodiaptomus vicinus and T. cf. lanaonus might lead to the finding of spurious inter-clades gaps, which would bias subsequent DNA taxonomy analyses (e.g. 59,60 ). This is particularly evident for Tropodiaptomus cf. lanaonus, whose four clades of putative species rank correspond in fact to four clusters of geographically isolated populations (Fig. 5b), so that the occurrence of intermediate haplotypes in geographically intermediate areas, which would reveal a scenario of a clinal distribution of their genetic diversity, cannot be excluded. Moreover, even in the absence of intermediate haplotypes, the genetic differentiation of allopatric populations might be ascribable to local adaptation, genetic drift, or simple isolation by distance phenomena instead of to actual speciation processes. Conversely, it is not likely that a wider sampling effort dedicated to T. vicinus would change the paraphyletic status suggested for the Thai populations of this morphospecies (Fig. 3), which seems to actually harbour a relevant cryptic diversity, as also suggested by the absence of a clear geographical pattern of molecular diversity (Fig. 5a). In fact, when genetically distinct clades observed within a morphospecies are found in sympatry or constitute a paraphyletic group, this provides strong indirect evidence that these entities are actually independent evolutionary units (cf. 9 ).
Pending other studies and evidence, we thus here refrain from considering the highlighted clades within Tropodiaptomus vicinus and T. cf. lanaonus as taxa of species rank, at the same time stressing the urgency to an in-depth analysis of the diversity pattern observed within these two morphospecies with dedicated surveys covering the whole potential distribution area of these taxa in the Oriental region.
The apparent decoupling between morphological-and molecular-based assessments of Thai Tropodiaptomus diversity is here ascribed to the morphological conservatism of the studied taxa, coupled with an inadequate understanding of the taxonomical value of the characters traditionally used when describing or identifying species. Accordingly, the taxonomic diagnostic morphological characters proposed in recent keys to the genus (e.g. 13,18,22,47 ) should be used with caution and coupled, when possible, with a molecular characterization of the studied populations.
Concluding remarks. Currently knowledge of the distribution of Thai diaptomid species is largely incomplete and more dedicated studies should be realised to get a sound picture of the actual distribution of the studied taxa. This problem affecting biodiversity studies, known as "Wallacean shortfall" 55 , severely hinders our understanding of the diversity, ecology, and natural history of diaptomids. Present results thus stress that the diversity of Thai diaptomid copepods is currently largely underestimated due both to the Linnean and Wallacean shortfalls. Stabilizing selection and the focus of traditional morphological studies on a few "classical" characters along with the practical difficulties linked with sampling throughout vast and sometimes difficult-to-reach areas, might in fact have been leading to a gross underestimate of the diversity of Diaptomidae fauna of Thailand. The realization of further systematic faunal surveys in currently undersampled areas, coupled with the integrative morphological and molecular study of the collected species, is desirable before any taxonomical act involving cryptic species is done.
We hope that this work might pave the way to inspire future work aimed at a better knowledge of Thai diaptomids.