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

Abstract

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.

About the Authors:

Matteo Montagna

* E-mail: [email protected]

Affiliation: Dipartimento di Scienze Agrarie e Ambientali—Università degli Studi di Milano, Via Celoria 2, I-20133, Milano, Italy

Valeria Mereghetti

Affiliation: Dipartimento di Scienze Agrarie e Ambientali—Università degli Studi di Milano, Via Celoria 2, I-20133, Milano, Italy

Valeria Lencioni

Affiliation: MUSE—Museo delle Scienze, Corso del Lavoro e della Scienza 3, I-38122, Trento, Italy

Bruno Rossaro

Affiliation: Dipartimento di Scienze per gli Alimenti, la Nutrizione e l’Ambiente—Università degli Studi di Milano, Via Celoria 2, I-20133, Milano, Italy

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–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–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–24] and others in which the adopted approach failed to identify the species [25, 26].

Materials and Methods

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 phenol-xylene 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).

[Figure omitted. See PDF.]

Fig 1. Geographical location of collecting sites in Northern Italy.

Localities in which samples were collected are denoted by red dots while black squares indicate the cities of Turin (to the west) and of Milan (towards the center of the map). The inset shows a magnification of the collecting localities within the Rhaetian Alps.

https://doi.org/10.1371/journal.pone.0149673.g001

[Figure omitted. See PDF.]

Table 1. Analyzed species of chironomids.

https://doi.org/10.1371/journal.pone.0149673.t001

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–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⁸ 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 cross-validation (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⁵ 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 Matlab® vers. R2015a.

Results and Discussion

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. The 112 specimens analyzed belong to six genera of midges: Macropelopia Thienemann 1916, Diamesa Meigen 1835, Pseudodiamesa Goetghebuer 1939, Pseudokiefferriella Zavrel 1941; Sympotthastia Pagast 1947 and Orthocladius van der Wulp 1874. The species belonging to genera other than Diamesa were included in the analyses as outgroups. Among the 93 specimens ascribed to the genus Diamesa, 89 are morphologically attributed to the following eleven species: D. bertrami, D. cinerella, D. dampfi (Kieffer 1924), D. goetghebueri Pagast 1947, D. incallida, D. latitarsis (Goetghebuer 1921), D. modesta Serra-Tosio 1968, D. steinboecki (Goetghebuer 1921), D. tonsa, D. vaillanti Serra-Tosio 1972 and D. zernyi. In the case of D. tonsa, D. cinerella, D. zernyi, D. goetghebueri, D. bertrami, D. steinboecki and D. vaillanti adult males and pharate pupae, of unequivocal attribution, are present. Larvae belonging to the latitarsis group were identified from adult specimens collected from the same localities. Regarding the remaining four specimens: MR-36 is a larva with a yellow head capsule, which could be identified as D. cinerella; MR-98 showed a black head and could be identified as D. zernyi; while MR-13 and MR-113 possess intriguing features. MR-13 is a larva collected from the River Vermigliana (Baita Velon, Trento; 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).

[Figure omitted. See PDF.]

Fig 2. Micrographs of contrasting morphological characters harbored by the three discussed specimens.

The upper micrographs report the morphological characters of MR-13: labrum with bifid S_III setae (top left) and procercus with six anal setae (top right). Below are micrographs reporting details of the hypopygium, respectively of MR-115 (bottom left) and of MR-113 (bottom right) specimens.

https://doi.org/10.1371/journal.pone.0149673.g002

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).

[Figure omitted. See PDF.]

Fig 3. Species delimitation analysis based on cox1 gene sequences.

A Bayesian ultrametric tree inferred from the cox1 gene sequence dataset and used as input for GMYC and bGMYC models. Specimen identifiers are reported on tips (MR as an acronym of the collection identifier plus the id number); ^§: possible 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: black-edged 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.

https://doi.org/10.1371/journal.pone.0149673.g003

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 currently-used 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, while at the higher value (26 entities) the specimens MR-13 and MR-36 are recovered as separate units. Regarding the specimens belonging to zernyi group, GMYC and bGMYC achieved the same results obtained by ABGD, clearly discriminating specimens of D. zernyi (belonging to the same evolutionary unit) by the unit composed of D. tonsa, D. cinerella and D. vaillanti specimens. In addition, the same cases of discordance were confirmed: MR-14 and MR-40 identified initially as D. tonsa, and MR-115 a male pupa identified as D. vaillanti. The new morphological character (head capsule color gradient) developed in the present work (see paragraph below) revealed that MR-14 and MR-40 should be considered as D. zernyi. Morphological identification of specimen MR-115 is confirmed; nevertheless, as this organism exhibits contrasting characters it has been reported as a possible hybrid between D. vaillanti and D. tonsa (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).

[Figure omitted. See PDF.]

Fig 4. Bayesian consensus tree inferred from an alignment of 112 cox1 gene sequences.

On the nodes of main the lineages the support values are expressed as bpp (above) and aLRT (below); * denotes support values ≤ 0.65 bpp and ≤ 65% aLRT. Vertical dashed lines indicate species groups. The scale bar at the bottom indicates the distance in substitutions per site.

https://doi.org/10.1371/journal.pone.0149673.g004

[Figure omitted. See PDF.]

Fig 5. Pairwise Kimura two-parameter nucleotide distance.

a. Heat map of the K2P pairwise distance matrix; values of nucleotide distance are proportional to color intensity, with low and high values of pairwise nucleotide distance indicated respectively by dark and light colors; morphological species are reported on the axis. Box-plot representing intra- and inter-specific K2P (b) and observed (c) pairwise nucleotide distance.

https://doi.org/10.1371/journal.pone.0149673.g005

[Figure omitted. See PDF.]

Table 2. Within and between Kimura 2 parameter nucleotides mean distances and mean values of observed nucleotide differences^a.

https://doi.org/10.1371/journal.pone.0149673.t002

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. 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–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).

[Figure omitted. See PDF.]

Fig 6. A novel morphological diagnostic character: color gradient between submenti and genal setae.

Photo of Diamesa head capsule, the area of interest is highlighted by a rectangle. For each species on the left side is reported a micrograph of the head capsule; on the right side a graph reporting the RGB color profile of the analyzed region, embedded in the graph a picture reporting the color gradient from the analyzed specimens. A. Diamesa zernyi. B. Diamesa tonsa. C. Diamesa cinerella. SSm: setae submenti; S_9-10: genal setae.

https://doi.org/10.1371/journal.pone.0149673.g006

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²), 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.

[Figure omitted. See PDF.]

Table 3. DNA Barcoding statistics and performances.

https://doi.org/10.1371/journal.pone.0149673.t003

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

[Figure omitted. See PDF.]

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).

https://doi.org/10.1371/journal.pone.0149673.s001

(DOCX)

Acknowledgments

The authors thanks everyone involved in specimen collection, in particular Alessandra Franceschini and Luca Toldo; Simon Pierce and David Neale for English revision.

Author Contributions

Conceived and designed the experiments: MM BR. Performed the experiments: MM BR VM VL. Analyzed the data: MM BR VM VL. Contributed reagents/materials/analysis tools: BR VL. Wrote the paper: MM. Discussed the results, revised and commented on the manuscript: MM BR VM VL.

Citation: Montagna M, Mereghetti V, Lencioni V, Rossaro B (2016) Integrated Taxonomy and DNA Barcoding of Alpine Midges (Diptera: Chironomidae). PLoS ONE 11(3): e0149673. https://doi.org/10.1371/journal.pone.0149673

References

1. Dirnböck T, Essl F, Rabitsch W. Disproportional risk for habitat loss of high-altitude endemic species under climate change. Glob Chang Biol. 2011;17: 990–996.

2. Beniston M, Stoffel M. Assessing the impacts of climatic change on mountain water resources. Sci Total Environ. 2014;493: 1129–1137. pmid:24360916

3. Rossaro B, Lencioni V, Boggero A, Marziali L. Chironomids from southern Alpine running waters: ecology, biogeography. Hydrobiologia. 2006;562: 231–246.

4. Brown LE, Hannah DM, Milner AM. Vulnerability of alpine stream biodiversity to shrinking glaciers and snowpacks. Glob Chang Biol. 2007;13: 958–966.

5. Jacobsen D, Milner AM, Brown LE, Dangles O. Biodiversity under threat in glacier-fed river systems. Nat Clim Chang. 2012;2: 361–364.

6. Jacobsen D, Cauvy-Fraunie S, Andino P, Espinosa R, Cueva D, Dangles O. Runoff and the longitudinal distribution of macroinvertebrates in a glacier-fed stream: implications for the effects of global warming. Freshwater Biol. 2014;59: 2038–2050.

7. Rossaro B, Montagna M, Lencioni V. Environmental traits affect chironomid communities in glacial areas of the Southern Alps: evidence from a long-lasting case study. Insect Conserv Divers. 2016.

8. Lods-Crozet B, Lencioni V, Olafsson J, Snook D, Velle G, Brittain J, et al. Chironomid (Diptera: Chironomidae) communities in six European glacier-fed streams. Freshwater Biol. 2001;46: 1791–1809.

9. Lencioni V, Rossaro B. Microdistribution of chironomids (Diptera: Chironomidae) in Alpine streams: an autoecological perspective. Hydrobiologia. 2005;533: 61–76.

10. Pagast F. Systematik und Verbreitung der um die Gattung Diamesa gruppierten Chironomiden. Archiv für Hydrobiologie. 1947;41: 435–596.

11. Hansen DC, Cook EF. The systematics and morphology of the Nearctic species of Diamesa Meigen, 1835 (Diptera: Chironomidae) Mem Am Entomol Soc. 1976;30: 1–203.

12. Rossaro B, Lencioni V. A key to larvae of Diamesa Meigen, 1835 (Diptera: Chironomidae), well known as adult males and pupae from Alps (Europe). J Entomol Acarol Res. 2015;47: 123–138.

13. Serra-Tosio B. Ecologìe et biogéographie des Diamesini d’Europe (Diptera, Chironomidae). Travaux du Laboratorie d’ Hydrobiologie et de pisciculture de l’Université de Grenoble. 1973;63: 5–175.

14. Willassen E. A review of Diamesa davisi Edwards and the davisi group (Diptera, Chironomidae). In: Fittkau EJ (ed) Beiträge zur Systematik der Chironomidae, Diptera. Spixiana Supp. 1985;11: 109–137.

15. Thienemann A. Bestimmungstabelle für die Larven der mit Diamesa nächst verwandten Chironomiden. Beitr Entomol. 1952;2: 244–256.

16. Rossaro B. Description of some unknown larvae of Diamesa genus and corrections to previous descriptions (Diptera, Chironomidae). Arch Hydrobiol. 1980;90: 298–308.

17. Carew ME, Hoffmann AA. Delineating closely related species with DNA barcodes for routine biological monitoring. Freshwater Biol. 2015;60: 1545–1560.

18. Lin X, Stur E, Ekrem T. Exploring Genetic Divergence in a Species-Rich Insect Genus Using 2790 DNA Barcodes PLoS One. 2015;10: e0138993. pmid:26406595

19. Sharley DJ, Pettigrove V, Parson YM. Molecular Identification of Chironomus spp. (Diptera) for biomonitoring of aquatic ecosystems. Aust J Entomol. 2004;43: 359–365.

20. Carew ME, Pettigrove V, Hoffmann AA. The utility of DNA markers in classical taxonomy: using cytochrome oxidase I markers to differentiate Australian Cladopelma (Diptera: Chironomidae) midges. Ann Entomol Soc Am. 2005;98: 587–594.

21. Carew ME, Pettigrove V, Cox RL, Hoffmann AA. DNA identification of urban Tanytarsini chironomids (Diptera: Chironomidae). J North Am Benthol Soc. 2007;26: 587–600.

22. Pfenninger M, Nowak C, Kley C, Steinke D, Streit B. Utility of DNA taxonomy and barcoding for the inference of larval community structure in morphologically cryptic Chironomus (Diptera) species. Mol Ecol. 2007;16: 1957–1968. pmid:17444904

23. Al-Shami SA, Siti MN Md., Rawi CS, Ahmad AH. Preliminary study of phylogenetic relationship of rice field Chironomidae (Diptera) inferred from DNA sequences of mitochondrial cytochrome oxidase subunit I. Am J Appl Sci. 2009;6: 1004–1009.

24. Ekrem T, Stur E, Hebert PDN. Female do count: documenting Chironomidae (Diptera) species diversity using DNA barcoding. Org Divers Evol. 2010;10: 397–408.

25. Martin J, Guryev V, Blinov A. Population variability in Chironomus (Camptochironomus) species (Diptera, Nematocera) with a Holarctic distribution: evidence of mitochondrial gene flow. Insect Mol Biol. 2002;11: 387–398. pmid:12230538

26. Proulx I, Martin J, Carew M, Hare L. Using various lines of evidence to identify Chironomus species (Diptera: Chironomidae) in eastern Canadian lakes. Zootaxa. 2013;3741: 401–458. pmid:25113001

27. Wirth WW, Marston N. A method for mounting small insects on microscope slides in Canada Balsam. Ann Entomol Soc Am. 1968;61: 783–784.

28. Lencioni V, Marziali L, Rossaro B. I Ditteri Chironomidi morfologia, tassonomia, ecologia, fisiologia e zoogeografia. Quaderni del Museo Tridentino Scienze Naturali; 2007 pp. 1–172.

29. Andersen T, Cranston PS, Epler JH. Chironomidae of the Holarctic Region. Keys and diagnoses–Larvae. Insect Syst Evol Supp. 2013;66: 573.

30. Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol. 1994;3: 294–299. pmid:7881515

31. Porretta D, Mastrantonio V, Mona S, Epis S, Montagna M, Sassera D, et al. The integration of multiple independent data reveals an unusual response to Pleistocene climatic changes in the hard tick Ixodes ricinus. Mol Ecol. 2013;22: 1666–1682. pmid:23398505

32. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32: 1792–1797. pmid:15034147

33. Puillandre N, Lambert A, Brouillet S, Achaz G. ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Mol Ecol. 2012;21: 1864–1877. pmid:21883587

34. Pons J, Barraclough TG, Gomez-Zurita J, Cardoso A, Duran DP, Hazell S, et al. Sequence-based species delimitation for the DNA taxonomy of undescribed insects. Syst Biol. 2006;55: 595–609. pmid:16967577

35. Fujisawa T, Barraclough TG. Delimiting species using single-locus data and the Generalized Mixed Yule Coalescent (GMYC) approach: a revised method and evaluation on simulated datasets. Syst Biol. 2013;62: 707–724. pmid:23681854

36. Reid NM, Carstens BC. Phylogenetic estimation error can decrease the accuracy species delimitation: A Bayesian implementation of the General Mixed Yule Coalescent model. BMC Evol Biol. 2012;12: 196. pmid:23031350

37. Zhang J, Kapli P, Pavlidis P, Stamatakis A. A General Species Delimitation Method with Applications to Pylogenetic Placements. Bioinformatics. 2013;22: 2869–2876.

38. Montagna M, Sassi D, Giorgi A. Pachybrachis holerorum (Coleoptera: Chrysomelidae: Cryptocephalinae), a new species from the Apennines, Italy, identified by integration of morphological and molecular data. Zootaxa. 2013;3741: 243–253. pmid:25112986

39. Cranston PS, Krosh MN. DNA sequences and austral taxa indicate generic synonymy of Paratrichocladius Santos-Abreu with Cricotopus Wulp (Diptera: Chironomidae). Syst Entomol. 2015;40: 719–732.

40. Lecocq T, Dellicour S, Michez D, Dehon M, Dewulf A, De Meulemeester T, et al. Methods for species delimitation in bumblebees (Hymenoptera, Apidae, Bombus): towards an integrative approach. Zool Scr. 2015;44: 281–297.

41. Kimura M. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. Journal of Molecular Evolution. 1980;16: 111120.

42. Collins RA, Boykin LM, Cruickshank RH, Armstrong KF. Barcoding's next top model: an evaluation of nucleotide substitution models for specimen identification. Methods Ecol Evol. 2012;3: 457–465.

43. Srivathsan A, Meier R. On the inappropriate use of Kimura‐2‐parameter (K2P) divergences in the DNA‐barcoding literature. Cladistics. 2012;28: 190–194.

44. Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61: 539–542. pmid:22357727

45. Darriba D, Taboada GL, Doallo R, Posada D. J ModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9: 772.

46. Lanave C, Preparata G, Saccone C, Serio G. A new method for calculating evolutionary substitution rates. J Mol Evol. 1984;20:86–93. pmid:6429346

47. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29: 1969–1973. pmid:22367748

48. Sanderson MJ. r8s: inferring absolute rates of molecular evolution, divergence times in the absence of a molecular clock. Bioinformatics. 2003;19: 301–302. pmid:12538260

49. Montagna M, Sassera D, Epis S, Bazzocchi C, Vannini C, Lo N, et al. "Candidatus Midichloriaceae" fam. nov. (Rickettsiales), an ecologically widespread clade of intracellular alphaproteobacteria. Appl Environ Microbiol. 2013;79: 3241–3248. pmid:23503305

50. Anisimova M, Gascuel O. Approximate likelihood-ratio test for branches: a fast, accurate, and powerful alternative. Syst Biol. 2006;55: 539–552. pmid:16785212

51. Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New Algorithms and Methods to Estimate Maximum-Likelihood Phylogenies: Assessing the Performance of PhyML 3.0. Syst Biol. 2010;59: 307–321. pmid:20525638

52. Hebert PDN, Cywinska A, Ball SL, De Waard JR. Biological identifications through DNA barcodes. Proc R Soc Lond B Biol Sci. 2003;270: 313–321.

53. Meyer CP, Paulay G. DNA barcoding: error rates based on comprehensive sampling. PLoS Biol. 2005;3: 2229–2238.

54. Sonet G, Jordaens K, Nagy ZT, Breman FC, De Meyer M, Backeljau T, et al. Adhoc: an R package to calculate ad hoc distance thresholds for DNA barcoding identification. ZooKeys. 2013;365: 329–336. pmid:24453565

55. Brown SDJ, Collins RA, Boyer S, Lefort MC, Malumbres-Olarte J, Vink CJ, et al. SPIDER: an R package for the analysis of species identity and evolution, with particular reference to DNA barcoding. Mol Ecol Resour. 2012;12: 562–565. pmid:22243808

56. Oksanen J, Kindt R, Legendre P, O’Hara B. Vegan: community ecology package. R package version 1.8–5. 2007. Available: http://cran.r-project.org/.

57. Sæther OA. A glossary of chironomid morphology terminology (Diptera: Chironomidae). Entomol Scand Supp. 1980;14: 1–51.

58. Paterson HEH. Evolution and the Recognition Concept of Species. Johns Hopkins University Press, Baltimore, Maryland; 1993.pp. 234.

59. Arnqvist G. Comparative evidence for the evolution of genitalia by sexual selection. Nat. 1998;393: 784–786.

60. McPeek MA, Shen L, Torrey JZ, Farid H. The tempo and mode of three-dimensional morphological evolution in male reproductive structures. Am Nat. 2008;171: 158–178.

61. Serra-Tosio B. Description et écologie de Diamesa vaillanti n. sp. (Diptera, Chironomidae). Trav scient Parc natn Vanoise 1972; 2:9–25.

62. Mavárez J, Salazar CA, Bermingham E, Salcedo C, Jiggins CD, Linares M. Speciation by hybridization in Heliconius butterflies. Nat. 2006;441: 868–871.

63. Gómez-Zurita J, Sassi D, Cardoso A, Balke M. Evolution of Cryptocephalus leaf beetles related to C. sericeus (Coleoptera: Chrysomelidae) and the role of hybridization in generating species mtDNA paraphyly. Zool Scr. 2012;41: 47–67.

64. Michailova P. Cytogenetic analysis of a hybrid, Glyptotendipes pallens Mg. × Glyptotendipes glaucus Mg. (Diptera, Chironomidae): evolutionary considerations. J Zoolog Syst Evol Res. 1998;36: 1439–1469.

65. Nowak C, Hankeln T, Schmidt ER, Schwenk K. Development and localization of microsatellite markers for the sibling species Chironomus riparius and Chironomus piger (Diptera: Chironomidae). Mol Ecol Notes. 2006;6: 915–917.

66. Meier R, Shiyang K, Vaidya G, Ng PKL. DNA barcoding and taxonomy in Diptera: a tale of high intraspecific variability and low identification success. Syst Biol. 2006;55: 715–728. pmid:17060194