Integrated Taxonomy and DNA Barcoding of Alpine Midges (Diptera: Chironomidae)

Rapid and efficient DNA-based tools are recommended for the evaluation of the insect biodiversity of high-altitude streams. In the present study, focused principally on larvae of the genus Diamesa Meigen 1835 (Diptera: Chironomidae), the congruence between morphological/molecular delimitation of species as well as performances in taxonomic assignments were evaluated. A fragment of the mitochondrial cox1 gene was obtained from 112 larvae, pupae and adults (Diamesinae, Orthocladiinae and Tanypodinae) that were collected in different mountain regions of the Alps and Apennines. On the basis of morphological characters 102 specimens were attributed to 16 species, and the remaining ten specimens were identified to the genus level. Molecular species delimitation was performed using: i) distance-based Automatic Barcode Gap Discovery (ABGD), with no a priori assumptions on species identification; and ii) coalescent tree-based approaches as the Generalized Mixed Yule Coalescent model, its Bayesian implementation and Bayesian Poisson Tree Processes. The ABGD analysis, estimating an optimal intra/interspecific nucleotide distance threshold of 0.7%-1.4%, identified 23 putative species; the tree-based approaches, identified between 25–26 entities, provided nearly identical results. All species belonging to zernyi, steinboecki, latitarsis, bertrami, dampfi and incallida groups, as well as outgroup species, are recovered as separate entities, perfectly matching the identified morphospecies. In contrast, within the cinerella group, cases of discrepancy arose: i) the two morphologically separate species D. cinerella and D. tonsa are neither monophyletic nor diagnosable exhibiting low values of between-taxa nucleotide mean divergence (0.94%); ii) few cases of larvae morphological misidentification were observed. Head capsule color is confirmed to be a valid character able to discriminate larvae of D. zernyi, D. tonsa and D. cinerella, but it is here better defined as a color gradient between the setae submenti and genal setae. DNA barcodes performances were high: average accuracy was ~89% and precision of ~99%. On the basis of the present data, we can thus conclude that molecular identification represents a promising tool that could be effectively adopted in evaluating biodiversity of high-altitude streams.


Introduction
Recent climatic warming has had a strong impact on habitats and species occurring at high elevation [1]. Amongst other effects, extensive environmental change occurs during glacial retreat, which affects hydrological and thermal regimes of glacier-fed streams [2]. Stream biodiversity is expected to dramatically change in relation to retreating glaciers, favoring an upstream shift of lowland species associated with local extinction of kryal species [3][4][5][6][7]. Unfortunately, a monitoring method and dedicated biotic indices, able to evaluate the biotic components of high-altitude habitats, have not yet been developed. There are several reasons for this, including and most importantly difficulties in identifying larvae belonging to the genus Diamesa Meigen 1835 (Diptera; Chironomidae), which dominate species richness and abundance in glacial streams and cold spring habitats [7][8][9]. Identification of midges based on morphology can be accurately achieved for adult males [10][11] or to a lesser extent for pupal exuviae [10], as clear discriminating features are visible for these stages. In addition, at present, despite updated keys identifying Diamesa larvae being available [12], the implementation of innovative tools, able to accurately identify larvae of the genus Diamesa by integrating different sources of information (e.g., molecular and morphological diagnostic characters), are required. Such approaches will promote the exploitation of ecological information provided by the presence/ absence of these species in the habitats under study [3].
The West Palaearctic species belonging to the genus Diamesa have been separated into nine different groups [13] according to adult male and pupal morphology. A combination of qualitative and quantitative characters observable in larvae (head capsule color, mouth-parts and posterior body appendages) only partially confirm this separation. This is because D. aberrata Lundbeck, 1898 and D. incallida (Walker 1856), included by Serra-Tosio in the aberrata group, have very different larvae, suggesting the separation into two distinct groups, while D. bertrami, Edwards, 1935, included by Serra-Tosio in the zernyi-insignipes-cinerella groups, has a larva very similar to the larvae of the latitarsis group [12]. Within each group, determination to species level is generally hampered by the lack of diagnostic characters or by the degradation of valid taxonomic characters, such as mental and mandibular teeth, in field-collected samples [12]. Moreover, quantitative characters should be used with caution due to the intraspecific variability (even within the fourth larval instar) present in different populations adapted to different environmental conditions [14].
Larvae belonging to zernyi-insignipes-cinerella groups [13] share the presence of a very reduced procercus bearing four anal setae of moderate length (~200-300 μm) and short posterior pseudopods [12,15]. At present, species belonging to these groups are separated from each other only according to head capsule color, from yellow (D. insignipes Kieffer in Kieffer and Thienemann 1908 and D. cinerella Meigen in Gisti 1835) to dark brown (D. zernyi Edwards 1933 and D. vaillanti Serra-Tosio 1972). D. tonsa (Haliday in Walker 1856) represents an intermediate case, possessing a yellow head capsule with variably extended brown areas [12,15].
The larvae belonging to steinboecki, latitarsis, bertrami and aberrata groups differ to those of the zernyi-insignipes-cinerella group as they possess very elongated posterior pseudopods, reduced anal setae (<120 μm), and are characterized by the absence of procerci and a dark head capsule. In addition, identification of species belonging to these groups is extremely difficult, despite recent proposal of valid diagnostic traits such as the number, length and diameter of anal setae and the shape of the mentum [12]. Only the dampfi and incallida groups are easily separated from other species at the larval stage: the former due to the presence of well developed procerci bearing six anal setae, the latter due to the characteristic annulation of the third antennal segment.
The present study is designed to test the congruence between species identifications on the basis of morphological diagnostic characters (i.e., morphospecies), such as head capsule color and anal setae length [12,15,16], and putative molecular species (operational taxonomic units, evolutionary species and phylospecies) delimitated using different approaches on the basis of a single-gene marker (i.e., the mitochondrial cytochrome oxidase subunit I-cox1). In addition, the effectiveness of DNA barcoding for species-level identification is tested. In recent years, molecular-based approaches have been successfully adopted to delimit midge species [17,18] and facilitate species identification (e.g., DNA barcoding studies). Contrasting results have been achieved by DNA barcoding, including cases in which its utility has been demonstrated [17][18][19][20][21][22][23][24] and others in which the adopted approach failed to identify the species [25,26].

Ethics Statement
No species of Diptera Chironomidae are listed in any national or regional law as protected or endangered. All the specimens were collected in state-owned properties. The collection of these invertebrates is not subjected to restriction by Italian law and does not require permission; permission to collect biological specimens in protected areas was provided by the competent authorities (prot. N. 2342/V/8/2-2014; prot. N. 2598/10.10-2015).

Sampling, Specimen Manipulation and Morphological Identification
Chironomid samples (larvae, pupae and adults) were collected by using drift net, Surber net and malaise traps during several field trips between 2013 and 2015 in nine localities within the Alps and Apennines (Table 1; Fig 1). The collected specimens were directly placed in absolute ethanol and sorted to the genus level in the laboratory by stereomicroscopy (Leica DM LS B2 and Leica MS 5). DNA was extracted from the body of full-grown larvae (4 th instar) after the removal of the head capsule and the caudal part, while DNA was extracted from pupae and adults preserving the whole morphology. The removed larval parts were mounted on a microscope slide in Canada Balsam, after dehydration with acetic acid and clarification with phenolxylene 3:1 [27,28], then identified to the species level [12] on the basis of morphological features including all semaphoronts (adults, pupae and larvae; morphological species concept), whenever possible. Pupae and adult males were mounted as usual, and identified using available identification keys [10,29]. Species ecology and distribution, as well as association of larvae with pupae and adults collected in the same locality were also considered as additional information to identify the specimens, e.g., adult males of D. insignipes have never been collected within the Alps [3], and so larvae with a yellow head were not assigned to this species. Measures were acquired by using optical microscopy at different magnifications (× 40 -× 1000), including photography using a digital camera (Leica DFC320).

DNA Extraction, PCR Amplification and Sequencing
DNA was extracted using DNeasy Blood and Tissue Kit (Qiagen, Heidelberg) following the manufacturer's instructions. A fragment of 658 bp of the mitochondrial cox1 gene was amplified by PCR using universal primers for metazoa LCO1490/HCO2198 [30]. The concentration of reagents used for cox1 amplification and thermal profile followed [31]. Successful amplification was determined by gel electrophoresis and PCR products were bidirectionally sequenced by ABI technology (Applied Biosystems, Foster City, CA, USA). The obtained electropherograms were manually edited and assembled into a consensus sequence using Geneious Pro 5.3 (Biomatters Ltd., Auckland, New Zealand); consensus sequences were deposited in ENA archive (accession numbers: LN897576-LN897687).

Bioinformatic and Species Delimitation Analyses
The obtained cox1 gene sequences were checked and aligned at the amino acid level using MUSCLE [32] and then back translated to the nucleotide sequence. The previously obtained alignment was used as input for species delimitation analyses. Independent methods requiring no a priori information on the existing morphospecies were adopted: i) the automatic barcode gap discovery tool (ABGD; [33]), which attempts to delimit species (here equivalent to operational taxonomic units) by estimating the optimal distance threshold (OT) for the given set of data; ii) coalescent tree-based methods as the generalized mixed Yule-coalescent model (GMYC; [34,35]) associated with its Bayesian implementation (bGMYC) [36] and the Poisson tree process model (PTP; [37]) in order to identify phylospecies and evolutionary species. Molecular approaches delimiting evolutionary units have been successfully adopted in several case studies in insects [38][39][40]. ABGD analyses were performed using the web-based interface (http://wwwabi.snv.jussieu.fr/public/abgd) with the Kimura-2-parameter model (K2P; [41]) as the model of nucleotide evolution. Prior maximum divergence of intraspecific diversity was settled from the value corresponding to a single nucleotide difference (i.e., 0.00153) to 0.1, relative gap width of 0.5. The remaining parameters were left with default settings. Despite the extensive use of K2P nucleotide distance in the scientific literature, this measure be inadequate to properly delimit species [42][43]. In order to avoid such problems a further ABGD analysis was performed, adopting uncorrected nucleotide distance and the results were compared with those of previous analyses. The single threshold GMYC method was implemented in the R package "splits" (SPecies LImits by Threshold Statistics) while the bGMYC method was performed in the R package "bGMYC". Bayesian inference analysis was performed by MrBayes 3.2 [44] in order to obtain the phylogram used, after conversion in ultrametric (see below for the adopted procedure), as input for GMYC and bGMYC analyses. Nucleotide substitution models were estimated using jModelTest 2 [45] and the model best-fitting the sequence was selected as General Time Reversible (GTR; [46]) with gamma distribution and proportion of invariable sites according to Bayesian Information Criterion. Two independent runs were performed using the following parameters: length of the Markov chain settled to 1 Ã 10 8 generations; trees and parameters sampled every 1000 generations; models of nucleotide evolution as obtained by model selection. The convergence of the two runs was checked using Tracer [47] and the burn-in fraction estimated accordingly. The Bayesian majority-rule consensus tree was converted to ultrametric in r8s 1.7 [48] using penalized likelihood with a smoothing parameter of 0.1, selected after crossvalidation (as described in [38,49]). The coalescent tree-based PTP method was performed using the web interface available at http://species.h-its.org/ptp/ with the following parameters: the Bayesian majority-rule consensus non-ultrametric tree as input, 5 Ã 10 5 MCMC generations, thinning every 100 generations, burning fraction = 0.20.
Maximum likelihood tree was inferred, adopting previous model parameters and approximate likelihood ratio test as node support (aLRT; [50]), by using PhyML [51].

DNA Barcoding and Nucleotide Distance Matrix
In its original meaning, DNA barcoding is designed to identify organisms on the basis of a DNA sequence adopting a fixed threshold of nucleotide distance [52]. In order to increase the success of specimen identification, the optimal intra-interspecific nucleotide distance threshold (OT; [53,54]) estimated from the analyzed dataset of sequences was adopted instead of a fixed, a priori defined, threshold. OT corresponds to the values of nucleotide distances at which the sum of false positive (FP; type I errors) and false negative (FN; type II errors) identifications reached minimum values.
All DNA barcoding analyses, including OT optimization, were performed on different cox1 sequence datasets (hereafter reported as ds plus a number from one to six on the basis of their features) using functions implemented in the R package "spider" [55]. For each dataset, a K2P [41] distance matrix was calculated. The accuracy and precision of DNA barcoding was calculated on the basis of the obtained data as defined by [54].
Pairwise nucleotide mean distance, box plots and heat map were estimated using the R package vegan [56], K2P [41] was adopted as the model of nucleotide substitutions.

New Diagnostic Character and Image Analysis
A new diagnostic character, represented by the brown-yellow color gradient in the area joining setae submenti and genal setae [57], has been considered as operative in identifying larvae of D. zernyi, D. tonsa and D. cinerella. The ventral and dorsal part of the head capsule were separated with fine tungsten needles and mounted so that the area between setae submenti and genal setae was easily visible. The RGB color profile of the segment joining setae submenti and genal setae was analyzed using the functions imread.m, imshow.m and improfile.m from the Image Analysis toolbox of Matlab1 vers. R2015a.

Morphological Identification of Analyzed Specimens
The DNA was extracted from a total of 112 specimens (subfamily Diamesinae, and few Orthocladiinae and Tanypodinae as outgroups) collected in 16 localities in the Alps and Apennines (Fig 1). Morphological identification, geographical coordinates and altitude of collecting localities, developmental stages and the overall head capsule color (the last feature only for larvae belonging to zernyi and cinerella groups) are reported in Table 1 Table 1) exhibiting a yellow head capsule but harboring contrasting characters that hampered its identification. MR-13 exhibits six setae on each procercus and bifid S III setae on the labrum (Fig 2), the former feature suggests that this specimen should belong to the dampfi group while the latter suggests its ascription to the zernyi-cinerella groups. MR-113 is an adult male collected at the Amola glacier that, on the basis of morphological characters, resembles D. nowickiana Kownacki & Kownacka 1975 (Fig 2).

Species Delimitation Analyses
A fragment of 658 bp of the mitochondrial cox1 gene was obtained from 112 specimens (16 identified morphospecies and 10 specimens identified at genus level); no indels were observed.
The aligned cox1 gene sequences were subjected to ABGD analysis designed to delimit species estimating the OT from the data. The frequency distributions of pairwise K2P distance highlighted the existence of a clear gap in pairwise comparisons (Fig A, B in S1 File). The perfect match between the initial and the recursive partitions occurred at nucleotide divergence values ranging from 0.7% to 1.4% and twenty-three groups (or putative molecular species) were identified (Fig C in S1 File). ABGD analysis implementing the observed nucleotide distance lead to comparable results and 26 groups were identified at the match between initial and recursive partitions. These results showed a high level of congruence between groups identified on the basis of cox1 gene sequences and the identified morphospecies. Specifically, all species belonging to D. steinboecki, latitarsis, bertrami, dampfi and incallida as well as all outgroup species are recovered as members of separate entities by ABGD analysis performed with K2P distance, perfectly matching morphospecies. Whereas, within the zernyi-cinerella groups, all specimens morphologically identified as D. zernyi, D. tonsa and D. cinerella where grouped into two clusters: i) a group composed by specimens identified as D. tonsa and D. cinerella (including adult male of both species), five unidentified larvae at the 3 rd developmental stage, a larva identified as D. zernyi according to head capsule color (MR-9) and a male pupa identified as D. vaillanti (MR-115; Fig 2); ii) a group composed of specimens of D. zernyi (with an adult male), all larvae initially identified as D. zernyi on the basis of the overall color of the head capsule. Thus, on the basis of the adopted distance-based approach, D. tonsa-D. cinerella-D. vaillanti (only one) specimens of certain morphological identification (adult males and larvae of clear attribution) belong to the same unit. Interestingly, the specimen MR-13 showed contrasting characters (MR-13; Fig 2) and, being collected from the River Vermigliana clustered with MR-36 from the same locality, in a single, separate group. ABGD analysis, performed using observed nucleotide distance, identified the specimens D. zernyi MR-40, D. steinboecki MR-32 and D. tonsa MR-108 as entities separated from groups harboring conspecific specimens. Putative molecular species recovered by ABGD analysis adopting the K2P model of evolution are more congruent with morphology with respect to those achieved by the same approach using observed nucleotide distance.
Species delimitation analyses performed by implementing the coalescent tree-based approach (i.e., GMYC, bGMYC and bPTP) led to almost identical results but some differences were apparent relative to ABGD (Fig 3). The GMYC model exhibited a significantly better likelihood than the null model (p-value < 0.001; logL GMYC = 612.6, logL NULL = 575.5), indicating that a boundary between and within species was identified. Twenty-five maximum likelihood entities (95% CI [24,26]) were identified at the threshold between Yule and Coalescent models (Fig 3). Similar results were obtained by bGMYC, which identified 25-26 evolutionary units, and by the bPTP method with 26 maximum likelihood partitions (estimated number of species between 23 and 34, average 26.2) (Fig 3).
Discrepancy relative to the distance-based ABGD was recovered for specimens of D. goetghebueri, for which three separate evolutionary units were identified. No complete congruence between collecting localities and clustering pattern was recovered. Indeed, at the Amola glacier organisms belonging to all three identified evolutionary units of D. goetghebueri coexist in sympatry, whereas at Trobio and Careser only specimens belonging to one unit were found. The result that at Amola at least three independent lineages of D. goetghebueri coexist, which do not possess a recent common ancestor, could be interpreted as the result of the antiquity of this population or as the result of repeated events of colonization by organisms from different populations. Due to the small sample size and to the use of a single-gene marker, we cannot reach any reliable conclusions on the basis of the present data. A possible alternative scenario could be that larval stages of species phylogenetically close to D. goetghebueri, such as D. lindrothi Goetghebuer 1931 and D. laticauda Serra-Tosio 1964, are not distinguishable by currentlyused morphological characters but segregate at the molecular level. At the lower value of the GMYC confidence interval (24 entities) two evolutionary units of D. goetghebueri collapse,  (Fig 2; see paragraph below). Interestingly, on the basis of cox1 sequences, it is not possible to discriminate between specimens of D. tonsa and D. cinerella: the two taxa were determined to be paraphyletic on the basis of the cox1 gene (Figs 3 and 4) and possess values of pairwise K2P nucleotide distance (average 0.94%, SD = 0.22%) comparable with the average of intraspecific nucleotide distance (K2P-intra avg = 0.88%, SD = 0.64%; K2P-inter avg = 11.79%, SD = 3.58%; Fig 5, Table 2). The graphical representation of pairwise K2P nucleotide distance matrix through the heat map allows the identification of a group of specimens on the basis of their pairwise K2P nucleotide distance values (Fig 5A). Comparisons between morphologically conspecific specimens are denoted by darker boxes (low values of pairwise nucleotide distance) while non-conspecific comparisons, with some already discussed exceptions, are characterized by light boxes (high values of pairwise nucleotide distance; Fig 5). The non-overlapping distribution of intra-and inter-specific pairwise K2P/observed distances confirmed the existence of a clear gap in pairwise comparisons (box plots in Fig 5A-5C).

Topology Inferred from Cox1 Gene Sequences
Although a single DNA marker can fail to produce a reliable phylogeny among organisms, we believe that the results achieved in the present study on the basis of cox1 gene sequences have merit. A Bayesian consensus tree was inferred as an input for the tree-based species delimitation methods (GMYC, bGMYC and bPTP; Fig 4). Interestingly, the inferred topology clearly determines the close relationships of taxa belonging to the same species group, as defined by morphological synapomorphies exhibited by males and pupal exuviae [13]. All the determined species groups were well supported (BI!0.97, aLRT!91). In contrast, relationships among the species groups are not resolved by cox1 gene sequences. These results could be explained by the limitation of cox1 in recovering the cladogenesis among the species groups under analysis or, alternatively, that almost simultaneous cladogenetic events led to the formation of the main species groups of Diamesa. The approach adopted in the present study is not adequate to test the latter hypothesis; further investigations is currently in progress.
On the basis of cox1 gene sequences over a total of 16 analyzed morphospecies, 11 were monophyletic while three are represented by a single specimen; the two sympatric species D.
hybrid specimens between D. vaillanti and D. tonsa; + : larvae at third instar. The vertical green line identifies the between/within species GMYC threshold. M: vertical black lines indicating the identified morphospecies. bG: putative species identified by bGMYC are represented by vertical solid colored boxes, colors indicate support values of Bayesian posterior probability (bpp) as follow: 0.05-0.5 in red, 0.5-0.9 in orange and 0.95-1 in yellow. G: vertical solid light-grey boxes represent putative species identified by GMYC. bPTP: blackedged boxes indicate the putative species (corresponding to the maximum likelihood partition) identified by the bPTP approach; values of bpp supporting putative species are reported, * = bpp of 1. Solid dark grey and light grey texture boxes indicate putative species identified by the ABGD approach, respectively implementing K2P and observed pairwise distance.    tonsa-D. cinerella were paraphyletic and clustered in a single group. The branching pattern obtained for D. tonsa-D. cinerella is congruent with a scenario of recent origin for these two taxa, in which a complete lineage sorting has not yet been achieved; episodes of gene flow between these species represent a further possible explanation. The clear differences in the morphology of male genitalia, associated with a small between-taxa nucleotide mean distance (0.94%, SD = 0.22%; Fig 5A, Table 2), is further evidence that morphological features could evolve more rapidly than neutral/semi neutral genetic markers [57][58][59][60].

New Diagnostic Character and Image Analysis
The cases of discrepancies observed within zernyi and cinerella groups with regard to morphological and molecular identification methods prompted us to re-analyse slides of all specimens to accurately explore the color of the head capsule. Detailed analyses lead to the discovery of a new and, in our view, more accurate diagnostic character represented by the color gradient from the genal setae to the setae submenti (Fig 6). Reared larvae of D. zernyi exhibit a color profile from a darker color in the genal region (from genal setae) to a lighter color in the submental region (from submenti setae) (Fig 6A), whereas those of D. tonsa exhibit the opposite trend with a lighter color in the genal region and a darker color in the submental region ( Fig  6B). The color gradient was not observed in reared larvae of D. cinerella: indeed they possess a consistent pale color in both genal and submental regions ( Fig 6C).
All specimens belonging to the zernyi group were then re-analyzed and identified according to the newly developed character. Interestingly, all specimens previously identified as not clustering with conspecific specimens were misidentified according with the newly discovered character. Only the identification of mature male pupa MR-115, assigned on the basis of hypopygium to D. vaillanti, but clustering within the D. tonsa-D. cinerella clade on the basis of the cox1 gene sequence, produced a conflict. On the basis of the achieved results and considering the sister relationship between D. zernyi and D. vaillanti [61] we can hypothesize that the specimen MR-115 represents a hybrid between D. tonsa/D. cinerella and D. vaillanti. This result does not affect the status of the species since the capability of closely related taxa to hybridize is regarded as a plesiomorphic state that occurs among insects (e.g. [62,63]) and it has been demonstrated in Chironomidae (e.g., [64,65]). Crucially, the possible event of hybridization occurred at a very small glacier (area < 1 km 2 ), the Trobio, in the Orobian Alps, where the limited living and breeding habitats improve the probability of contact amongst organisms. Analyses that include information provided by nuclear genes are required to rigorously test the possible hybridization event.

DNA Barcoding
A total of six sequence datasets were analyzed. Features and DNA barcoding performances of each dataset are reported in Table 3. The analyses for the estimation of intra-/inter-specific nucleotide OT achieved contrasting results depending on the dataset analyzed ( Table 3). The estimated OT ranges from a minimum value of 0.7% in the case of ds2, where sequences of larvae at the 3 rd developmental stage and singletons were excluded, to a maximum value of 1.4% for ds1, where only sequences of larvae at the 3 rd developmental stage were excluded. Values of OT estimated from Alpine non-biting midges included in this study were much lower than those obtained for the delineation of species belonging to the genus Tanytarsus (4-5%; [18]). For the estimated OTs the cumulative error of misidentification ranges from 0, in the case of ds4 and ds6, to 26 in the case of ds1. Twenty-five out of the 26 misidentifications are due to specimens morphologically identified as D. tonsa and D. cinerella. Previous results can be explained by: i) the value of pairwise K2P nucleotide mean distance between the two species (0.94%; Table 2, Fig 5) being lower than the estimated OT (1.4% in the case of ds1), and by ii) the paraphyletic status of the two species on the base of cox1 gene sequences (Figs 3 and 4). The remaining case of misidentification is represented by the apparent hybrid between D. vaillanti and D. tonsa. Nevertheless, the near neighbor analysis [66] highlighted that the majority of tested sequences (from~80% to 100%) showed a conspecific individual as closest. Overall, the DNA barcoding approach on the six analyzed datasets achieved on average an accuracy of 89% [74%, 100%] and a precision of 99% [92%, 100%]. Results obtained by the DNA barcoding approach on the midge species under study here are very promising and confirm its enormous utility in supporting rapid and large-scale for the evaluation of insect biodiversity of high-altitude stream habitats.

Conclusion
The present study, mainly focused on testing the congruence between species identified using "traditional" morphological characters and putative molecular species (identified by a~650 bp fragment of the mitochondrial cytochrome oxidase I), has demonstrated an almost complete congruence in the results achieved by both approaches. Cases of discrepancies between the two methods were recovered within zernyi and cinerella groups, where some larval specimens were found to be misidentified on the basis of the traditionally used morphological characters (i.e., the overall color of the head capsule) and a possible hybrid between D. vaillanti and D. tonsa was collected at Trobio glacier. In addition, identification methods based on cox1 gene sequences failed to distinguish between specimens belonging to the clade D. tonsa-D. cinerella as these two taxa were found to be paraphyletic on the base of this marker. Further analyses, based on a multi-gene approach or on more innovative methods such as RAD sequencing, are required to disentangle the intricate relationships between these two sister species.
The above-mentioned critical cases determined within zernyi and cinerella groups prompted us to analyze more carefully the color of the head capsule. This larval character, despite being influenced by several factors such as the developmental stage and specimens conservation and preparation, has been extensively used to identify larvae of Diamesa [12,15,16]. The procedure led to the discovery of a more accurate trait related to head capsule color. Larvae of D. zernyi, D. tonsa and D. cinerella, in the 4 th developmental stage, are more accurately identifiable on the base of the color gradient between the setae submenti and genal setae. This character, focusing on the color gradient, in influenced to a lesser extent, with respect to overall head color, by storage conditions and specimens preparation. We believe that it is important to remark that the discovery of the new diagnostic character has been possible only because species delimitation analyses performed on molecular data highlighted cases in which information provided by molecules were in contrast with those supported by morphology. This result represents further evidence that "traditional" taxonomy benefits from the molecular tools and that conclusive results can only be achieved by adopting integrated approaches. On average, performances of molecular identifications through DNA barcoding were found to be elevated, with a recovered accuracy in specimen identification of~89% and a precision of 99%. These values reached 100% after the removal of specimens identified as D. tonsa or as D. cinerella and of the possible hybrid specimen recovered at the Trobio glacier. The results allow us to conclude that the cox1 gene sequence is an useful aid in species identification and paves the way for the use of molecular taxonomy, through DNA barcoding or DNA metabarcoding protocols, in support of biological studies aiming to monitor and evaluate the biodiversity of midges, and more generally invertebrates, inhabiting high-altitude streams and cold spring habitats. Nevertheless, it is fair to remember that in some cases (i.e., D. tonsa and D. cinerella) molecular tools fail in specimen identification, and thus the support of specialist entomologists is still required.
Supporting Information S1 File. Automatic Barcode Gap Discovery analysis. Histogram of pairwise nucleotide distances (Fig A). Plot of ranked pairwise nucleotide distances (Fig B). Automatic partition of the analyzed data set (Fig C). (DOCX)