Since 1981, the biocide Bacillus thuringiensis var. israelensis (Bti) is widely applied for mosquito (Culicidae, Diptera) control in temporary wetlands of the German Upper Rhine valley to minimize nuisance of local residents (Becker 1998). Bti is considered as the most environmentally friendly alternative to chemical pesticides for mosquito control due to a supposedly high specificity to mosquito larvae and negligible non-target effects even on closely related dipterans (Boisvert and Boisvert 2000). This is important as large areas of both aquatic and terrestrial habitats of the Upper Rhine valley are protected (bird sanctuaries, nature reserves and Natura 2000 sites) and comprise of biodiversity hotspots (Biggs et al. 2005, Lukács et al. 2013).

However, several studies have shown that Bti non-target effects are possible (reviewed in Boisvert and Boisvert (2000)). Non-biting midges (Chironomidae, Diptera) are the most Bti-sensitive non-target family (Boisvert and Boisvert 2000). Controlled experiments revealed varying mortality rates on chironomid larvae with older larvae being typically less sensitive to Bti (Ali et al. 1981, Treverrow 1985, Ping et al. 2005). They also reported different sensitivities among species (Yiallouros et al. 1999) and subfamilies (Liber et al. 1998). A recent study found that first instar larvae of Chironomus riparius are highly susceptible to Bti treatment even at commonly used mosquito control application rates (Kästel et al. 2017). Consequently, Bti application might overproportionately affect chironomid species in early larval stages at the time of application. So far, field studies have yielded ambiguous data on possible side effects of Bti on chironomid abundances. These range from positive effects on chironomid larvae richness possibly due to reduced mosquito competition (Lundström et al. 2010), over no effect on chironomid abundance (Lagadic et al. 2016), to a 35-80% reduction of chironomids abundances (Rodcharoen et al. 1991, Hershey et al. 1995, Vaughan et al. 2008, Poulin et al. 2010, Jakob and Poulin 2016).

Chironomids are a taxonomically and ecologically highly diverse group and often dominate all kinds of lotic and lentic ecosystems in terms of species abundances and biomass (Ferrington 2008). With sometimes over 50% of the total macroinvertebrate fauna in aquatic ecosystems (Milošević et al. 2013, Puntí et al. 2009) chironomids are thus a key food resource in wetland ecosystems. They also constitute a central link between aquatic and terrestrial food webs as adult midges are prey for birds, bats, spiders and adult dragonflies (Niemi et al. 1999, Stav et al. 2005, Poulin et al. 2010, Pfitzner et al. 2015). Furthermore, temporal and spatial variability in chironomid community composition has been observed (Lindegaard and Brodersen 1995, Rossaro et al. 2006, Milošević et al. 2013) together with a high adaptability of the community to changing environmental conditions (Raunio et al. 2011). Given these particular ecological characteristics, chironomids have been proposed as important indicators for monitoring freshwater ecosystems (Moog 2002, Sᴂther 1979). However, their morphological determination is very challenging and the taxonomic expertise needed for species identification of chironomids is often lacking (Batzer and Boix 2016). This makes it difficult to study changes in chironomid composition and utilize this as a monitoring tool. DNA-based determination approaches such as DNA barcoding seem therefore promising to support and complement the taxonomic assessment of chironomid community composition.

During recent years, DNA metabarcoding of whole communities has become a new powerful tool for environmental monitoring of aquatic ecosystems (Hajibabaei et al. 2011, Carew et al. 2013, Elbrecht and Leese 2015, Gibson et al. 2015). The DNA-based assays to monitor species biodiversity proved to be a rapid and efficient tool that allows the recovery of a substantial amount of taxa (Sweeney et al. 2011, Taberlet et al. 2012, Yu et al. 2012, Carew et al. 2013, Elbrecht et al. 2017). For metabarcoding in animals, typically the mitochondrial cytochrome oxidase subunit I (COI) is used (Hebert et al. 2003). The COI barcoding region has usually a good taxonomic resolution and comparatively well-curated databases as reference for many taxa (Sweeney et al. 2011, Elbrecht et al. 2017). The species present in the sample are identified based on a comparison of retrieved COI sequences (summarized as operational taxonomic units; OTUs) with reference databases (e.g. NCBI or BOLD; Ratnasingham and Hebert 2007). Good species coverage in the database is necessary for taxonomic assignment from sequences, but the identification rate for the different taxa vary widely (Ekrem et al. 2007, Kwong et al. 2012). For Chironomidae only about 30% of the estimated 700 different species in Germany (Samietz 1996) have an entry in the public databases BOLD with formal barcodes (accessed on 3. October 2017, search terms: chironomid & Germany). However, common taxa might be well represented.

In this study, DNA metabarcoding was applied to assess the distribution and species richness of chironomids in Bti-treated vs. first year Bti-untreated temporary wetlands in the Upper Rhine Valley. The study sites were part of a mosquito control area that has received regular Bti treatments for approximately 20 years (http://www.kabsev.de/1/1_2/1_2_1/index.php, accessed on 11. August 2017). Our first aim was to study the effectiveness of metabarcoding for chironomid diversity assessment as important and often overlooked freshwater bioindicator. We expected to obtain more species-level identifications based on molecular methods as compared to traditional taxonomic determination and used these data to calculate the saprobic index for the respective sites. Our second aim was to test for possible changes in the chironomid community composition of the temporary wetlands in response to mosquito control actions. Based on the above-mentioned studies we expected:

• an overall reduction in chironomid abundance at Bti-treated sites,
• a reduction in species richness of chironomids at Bti-treated sites, and therefore:
• an overall effect of Bti-treatment on chironomid community composition.

Study sites

The study was conducted in Rhineland-Palatinate in southwest Germany close to Neustadt-Geinsheim (Fig. 1). The study sites are regularly flooded in spring and dry out in summer. Thus, the area can be classified as seasonal (= temporary) wetland, which is moreover partly protected as a key amphibian breeding area in the region (Williams 2006). The area has been subject to regular mosquito control management actions for over 20 years, with usually one to two helicopter-applications of Bti between March and June, depending on temperature and precipitation. The study sites Fig. 1 were: "Stiftungsfläche" (S): mainly flooded grassland with some small permanent water bodies, "Großwald" (G): alder carr with larger permanent water bodies, "Mitteltrumm" (M): alder /oak carr with some deeper trenches and ditches and flat sinks. Additionally, the site "Lachen-Speyerdorf" (CL; see Fig. 1) served as control site and was located approximately 7 km away from the sites S, G and M. The site CL was dominated by open alder and pine forest with an abandoned river course.

For the first time after 20 years of regular Bti treatment, parts of the study area were left Bti-untreated in spring 2013 allowing for a split field design. Accordingly, S, G and M were divided into Bti-treated (T; 20 years treated) and untreated (U; first season untreated) site pairs, and CL served as control site never been treated with Bti. The helicopter application took place on April 10, 2013 using IcyPearls (Vectobac WG^®, ValentBiosciences) at a concentration of 1.44 x 10⁹ ITU/ ha.

Emergence data

Insect imagines were collected weekly with emergence traps (N = 5 per site and treatment, in total 35 traps, 0.25 m² area each) over a period of four months (April – July 2013) for 13 weeks after application (WAA; WAA 1 – WAA 13) of Bti. The preserved emergence was determined to order level and the order Diptera to family level using a Leica M80 microscope and a 10x magnification and counted per trap and sampling week. All chironomid specimens were selected for further analyses. All specimens were conserved in 70% ethanol and stored at room temperature for up to two years until DNA extraction.

Chironomid samples of all emergence traps per WAA were pooled for Bti-treated and Bti-untreated sites. For specific emergence peaks (N = 18, see Fig. 2 and text in results) these pooled samples were selected for metabarcoding to determine whether abundance differences can be attributed to a shift in species community composition.

Laboratory methods

Chironomids of emergence peaks (N = 18) were selected based on taxonomy and dried overnight at 60 °C. The specimens of each sample were grinded using a Tissue Lyser II (Qiagen, Hilden, Germany) at 30 Hz for 3 x 1 min with two sterile metal beads (3 mm, Hobbyfix, Opitec, Giebelstadt) with a brief centrifugation in between. DNA was extracted following a high salt DNA extraction protocol after (Aljanabi and Martinez 1997). Extraction success was verified using a Nanodrop (ND-1000 Spectrophotometer, Wilmington, USA). 50 µL of DNA from each sample were treated with 1.1 µL RNase (10 mg/mL, Roth, Karlsruhe, Germany) at 37°C for 30 min, followed by purification using a MinElute Reaction Clean up Kit (Qiagen, Hilden, Germany). The DNA concentration after clean-up was again measured using the Nanodrop and DNA concentrations of all samples were adjusted to 25 ng/µL.

A 322-bp fragment of the mitochondrial COI gene was amplified using the BF2 + BR1 primer set (Elbrecht et al. 2017). The used primer set was developed and evaluated with mock and in silico methods and does incorporate the needed degeneracy to amplify macroinvertebrates (including chironomids) reliably (Elbrecht and Leese 2017). The used fusion primers included Illumina adapter tails for sequencing (P5 or P7) and inline barcodes of different lengths for sample multiplexing (Elbrecht and Leese 2015). For each of the 18 samples two PCRs were conducted using the same primer pair but switching P5 and P7 Illumina adapters (Elbrecht and Leese 2015, Suppl. material 1). Sample 04CL was run with PCR replicates to test the PCR variablity. PCR reactions consisted of 1× PCR buffer (including 2.5 mM Mg²⁺), 0.2 mM dNTPs, 0.5 μM of each primer, 0.025 U/μL of HotMaster Taq (5Prime, Gaithersburg, MD, USA), 25 ng DNA, and HPLC H₂O to a total volume of 50 μL. The PCR program included the following steps: 94 °C for 3 minutes, 35 cycles of 94 °C for 30 seconds, 50 °C for 30 seconds, 65 °C for 120 seconds and ended with 65 °C for 5 minutes. PCR success was checked on a 1% agarose gel. Since some samples exhibited low DNA quantity (Qubit 2.0, Life Technologies, Carlsbad, CA, USA; measured concentration below 1 ng/µL), PCR for the respective samples was repeated with cycle number increased to 40 (Suppl. material 1). Amplicons were purified and size selected (retaining fragments of >300 bp) with a left-sided size selection using magnetic beads (SpriSelect, Beckman Coulter, Bread, CA, USA, ratio: 0.76x). The DNA concentration was quantified using the Qubit and a high sensitivity (HS) Assay Kit. Purified PCR products were pooled proportionately according to the number of specimens used in the extraction into a library to ensure all specimens are sequenced with the same sequencing depth. After pooling, the library was sent to an external laboratory (Macrogen, Seoul, Korea) for 300 bp paired-end sequencing on a MiSeq Illumina system (v3) run.

Bioinformatic analysis

Following the bioinformatic pipeline as previously described in Elbrecht and Leese (2017), the sequence data were processed as follows. In brief, after demultiplexing using a custom R script paired-end reads were merged to one sequence (Usearch version 8.8.1756). Primer sequences were removed via cutadapt (version 1.9.1). Singletons in each sample were before clustering OTUs with the cluster_otus command at 97% identity Edgar 2013. All samples (including singletons) were matched against the OTUs (Usearch). To enhance data reliability, sequences matched to the respective OTU had to occur in both replicates and exceed the 0.003% threshold sequence abundance to being considered in downstream analysis (Elbrecht and Leese 2015). Finally, the obtained OTU sequences were matched against the BOLD database to retrieve taxon identification. All used metabarcoding pipeline and R scripts are available inSuppl. material 2.

Using metabarcoding data for chironomid diversity assessment

The retrieved chironomid species list was checked for biogeographic and ecologically plausibility, i.e., if the species names were listed for Germany and are representative for temporary wetlands. If more than one species name per OTU was retrieved from BOLD with over 98% identity, we carefully examined the resulting hit table. For most hits, we then selected the biogeographically plausible species for our study region, based on the known biogeographical distribution and chironomid expert knowledge, for further ecological interpretations. If no clear decision could be made together with expert taxonomists, we followed the conservative approach to select the species name already represented in our data by another OTU. The species names retrieved in that way were categorized based on larval morphology in the context of standard water quality assessments into morphological determination “possible” (i.e., determination under 80x magnification without preparation), “difficult” (i.e., some characteristics need to be prepared and checked under greater magnification) and “impossible” (i.e., for species where the larva is not described, or do not show morphological differences within one genus, or would demand highly elaborative preparation for species determination). Using this approach, we aimed to elucidate which proportion of the chironomid species pool is neglected in standard water quality assessments, where only the easily and quickly determinable chironomid larvae are considered. We then calculated the percentage of species retrieved via metabarcoding (here emergence data) in relation to species, which would have also been possible to determine morphologically in standard water quality assessments (usually larval data). Based on available saprobic indices of the chironomid species (Moog 1995, Moog 2002) retrieved via metabarcoding we also calculated the chironomid saprobic index (SI; Table 2) exemplarily for the four study sites across all WAA. When the same species name was retrieved from more than one OTU the respective abundances were summed up.

Bti effects on chironomid community composition

To test whether the abundance data of emerged chironomids differed among Bti-treated and Bti-untreated sites (including the control site) for pooled samples of WAA 1 - WAA 4 (first emergence peak) and WAA 1 - WAA 13 (whole sampling period) a generalized linear mixed effect model (GLMM) (packages "nlme" v. 3.1-117, Pinheiro et al. 2016, and "MASS" v. 7.3-31, Venables and Ripley 2002) was implemented. As error structure the quasi Poisson family was chosen, where “study site” was implemented as random factor.

To test for differences in chironomid species richness between Bti-treated and Bti-untreated sites (including the control site), a Welch test was applied to the total number of OTUs and based on retrieved species names.

To test the hypothesis that the chironomid community composition differed between treatments, an adonis analysis (nonmetrical permutational MANOVA equivalent, Anderson 2001) was performed. For this, the Bray-Curtis distances of the Hellinger-transformed (see below) OTU read abundance assemblages per sites were calculated between pairs of sites, and these pairwise distances between sites were combined to a distance matrix of all sites, using the command ‘vegdist’ in the package “vegan” v. 2.4-1 (Oksanen et al. 2016). Within this distance matrix the nonmetrical permutational MANOVA equivalent was calculated using “Bti-treatment” as the distinguishing factor, with the command ‘adonis’ from the vegan package. We then assumed that time (WAA) would have a dominant effect on chironomid communities, but that Bti-treatment would alter the community as well. Therefore, also the interaction time * Bti-treatment was tested. Due to low number of replicates, the samples from the sites M and CL of WAA 8 and 9 were combined to the same time period. For the sites S and G only WAA 8 and WAA 9, respectively, were available.

All analyses were conducted in R (https://www.R-project.org). For all multivariate analyses, the Hellinger transformation was chosen to give less weight to the few high abundant OTUs, since the abundance data were highly left-skewed with few taxa reaching abundances several orders of magnitude higher than those of the less abundant species (Legendre and Legendre 2012).

Emergence data

In total, 11,589 emerged insects were collected, comprising of 17 taxa groups (Chironomidae: 78%; Culicidae: 14%; Trichoptera: 4%; Chaoboridae: 2%; Brachycera: 1%, other: 1%). On the Bti-treated sites 27 mosquito individuals were collected in the emergence traps, while on the Bti-untreated sites 1,006 mosquitoes emerged. Based on morphological identification 9,033 adult chironomids were collected. The number of chironomid specimens per emergence trap across all sites varied from 1 (Bti-treated site) to 1,239 (first year Bti-untreated site).

Emergence of chironomids fluctuated over time with varying emergence peaks at the Bti-untreated sites (Fig. 2). In particular, we detected one spring emergence peak (WAA 4) and two summer emergence peaks (WAA 9 + 10) at site G. At site S one spring peak (WAA 4) and two summer peaks (WAA 8 + 10) were detected. At site M, two summer peaks (WAA 10 + 13) were identified. A spring peak (WAA 4) and a summer peak (WAA 10) were identified at the never Bti-treated site CL. Specifically pooled emergence peak samples (N = 18, Figure 2) were selected for metabarcoding to investigate, if the abundance difference between Bti treated and untreated site pairs can be attributed to a shift in the chironomid community composition. The amount of individuals per pooled emergence peak sample varied from 22 to 541 (Suppl. material 1).

Bioinformatic analyses

In total 20,805,626 raw reads were generated by the MiSeq run with good read quality (Q30 ≥ 76.7% of reads). Raw data are available at NCBI SRA archive (accession number SRR4244505; https://www.ncbi.nlm.nih.gov/sra/?term=SRR4244505). After demultiplexing, merging and trimming of PCR primers 8,869,048 sequences were used for further analysis. The number of sequences in each sample was significantly correlated (p < 0.001, adj. R²=0.942) with the abundance of specimens per sample (Suppl. material 3).

OTU clustering analysis resulted in 442 OTUs. After application of the previously defined quality standards (0.003% minimum abundance) 89 OTUs were retained and used for subsequent analyses. The BOLD database searches identified 54 of the 89 OTUs (60.7%) as chironomids (Suppl. material 4). Of those, 38 OTUs (68.5%) could be assigned a species identification with 98 - 100% similarity, leaving 17 OTUs (31.5%) without species level identification (Suppl. material 4).

Using metabarcoding data for chironomid diversity assessment

In total 30 chironomid species were detected in the metabarcoding data set, with six species being assigned to 2-3 OTUs respectively (Table 1). For 11 OTUs we retrieved more than one species name with a sequence similarity of 97.82–100% (Suppl. material 4). For those OTUs, we were able to select the only biogeographically plausible species for our study region, based on biogeography and chironomid expert knowledge, for further ecological interpretations (underlined species names in Suppl. material 4). Only for OTU_12 two species names were biogeographically plausible, namely Chironimus luridus and C. riparius. Here, we selected the latter species as C. riparius was also characterized by other OTUs in our data, whereas C. luridus was not represented otherwise.

Of the 30 retrieved species, seven species can be routinely determined (cost and time efficient) based on larval morphology in ecological water assessments, whereas the remaining species are difficult (N = 5) or impossible (N = 18) to determine based on larval morphology (Table 1. This resulted in a 73% increase of retrieved chironomid species names based on metabarcoding (using emergence data) in relation to morphological larvae determination.

The chironomid SI was calculated based on the available saprobic value for 14 chironomid species (45.2%) of our data set Table 2. The saprobic value per detected species ranged from 0.8 (Monopelopia tenuicalcar, Xenopelopia nigricans) up to 3.5 (C. riparius). For our study sites the chironomid SI ranged from 1.3 at the control site to 3.3 at the Bti treated sites (Table 2).

Bti effects on chironomid community composition

The abundance of emergent chironomids until WAA 4 at the Bti-treated sites was reduced by 64.99% compared to the abundance in the Bti-untreated sites (GLMM t = 11.29, p = 0.008, df = 2). After WAA 13, slightly more chironomids hatched at the Bti-treated vs. Bti-untreated sites (2,132 vs 1,800 individuals, respectively). However, this difference was not statistically significant (GLMM t = -0.239884, df = 2, p = 0.833).

Neither the number of OTU per sample (Welch two sample t-test, t = 1.33, p = 0.20) nor the number of species assigned from the OTU based on the data base (Welch two sample t-test, t = 1.45, p = 0.17) were significantly different between Bti-treated and Bti-untreated site pairs.

The adonis model of crossed Bti-treatment * site effects (with the variation from time implemented as groups (strata) within which permutations are constrained) explained 51% of the variation of the multivariate chironomid community composition, 34% of which were due to the differences within sites (Table 3). The effect from Bti treatment suggests a statistically significant, but only minor component explaining 12% of the variation (p = 0.02; Table 3).

Using metabarcoding data for chironomid diversity assessment

With our metabarcoding approach we detected 54 chironomid OTUs across all study sites, of which almost 70% could be identified to species level using the BOLD database. Even though we did not have a specific reference database for our study system (e.g. Carew et al. 2013), we have mainly extracted biogeographically and ecologically meaningful species names as many of these species are frequently found in periodically desiccative ponds as euryoecious ubiquists (e.g., Ablabesmyia monilis, Acricotopus lucens, Chironomus riparius, Corynoneura scutellata, Cricotopus sylvestris, Limnophyes pentaplastus, Paratendipes albimanus, Polypedilum uncinatum, Psectrotanypus varius and Tanytarsus pallidicornis). Despite selecting only chironomid specimens for metabarcoding, also non-chironomid OTUs (N = 35; 39.3%) were detected in low abundancies (Suppl. material 4). These records included other Dipterans (N = 14; 15.7%), Trichopterans (N = 2; 2.2%), Lepidopterans (N = 2; 2.2%), Arachnids (Pionidae: N = 6; 6.7%), Fungi (Sporidiobolale, N = 2, Eurotiales, N = 1, Tremellales, N = 1 and Microstromatares, N = 1; 5.6% in total) and Bacteria (Rickettsia: N =1; 1.1%). In addition, OTUs without hit in the BOLD-database (N = 6, 6.7%) were detected. However, the presence of other taxa in the data is most likely related to the universality of the primer pair used, which could have amplified traces of other taxa previously stored in the same collection tubes, thus causing the large number of low level OTUs detected. Previous experiences show that low abundance OTUs are often derived from PCR or sequencing errors as well as chimeras and non-target DNA of small organisms (Elbrecht and Leese 2015, Elbrecht and Leese 2017, Elbrecht et al. 2017). As DNA was extracted from tissue bulk samples, eDNA and small non-target specimens are not of interest for this study. Thus, the data was not further analyzed. It should also be noted that the lower the OTU abundance is the more there are stochastic effects between samples, making it difficult to analyze OTUs with very low abundance.

For 11 OTUs we retrieved more than one species name with a sequence similarity of 97.82–100% (Suppl. material 4). This could indicate 1) limited taxonomic resolution of the short fragment amplified by the used primer set, producing only 322 bp COI fragments opposed to the 658 bp COI fragments using the classical Folmer primers (Folmer 1994); 2) different taxonomic keys used thus having different synonyms included; or 3) potential taxonomic misidentification which are also discussed in (Elbrecht et al. 2017) who recommend better data curation in taxonomic databases. In contrast, six chironomid species names were retrieved from two or three different OTUs, respectively, and two OTUs (45, 29) were assigned to questionable species names due to their biogeography (Procladius cf. fuscus and C. curabilis, Table 1). This could suggest cryptic intraspecific diversity, as we used a species divergence rate of 3% which might be too low for some species (Ekrem et al. 2007, Carew et al. 2013). As little is known about the genetic lineage of Chironomidae, cryptic species or variable phenotypes could be possible and cumber the correct identification by taxonomist (Anderson et al. 2013, Carew et al. 2007, Stur and Ekrem 2011).

The information benefit of metabarcoding by obtaining species names strongly depends on the quality of the database. For 17 chironomid OTUs (31.5%) no species identification could be obtained (Suppl. material 4). BOLD holds 270,292 published records of Chironomidae forming 5,540 BINs (clusters) with specimens from 49 countries. Of these records (accessed on 14.07.2017) , 100,231 have species names, and represent 1,233 species for an estimated species diversity of 15,000 worldwide (Armitage et al. 1995). For Germany, with an estimated species richness of approx. 700 different Chironomidae (Samietz 1996), BOLD has 3,706 published records forming 217 BINs (clusters). Of these records (accessed on 14.07.2017), 3,683 have species names, representing only 208 species (around 30%). Metabarcoding can only be as good as the database on which it relies for OTU matching to species identifications. We thus encourage experienced chironomid taxonomists to increase the number chironomid species in the BOLD database to even improve the effectiveness of metabarcoding for chironomid diversity assessments.

By applying metabarcoding we obtained 70% more chironomid species identifications than would have been possible based on traditional taxonomic determination of larval samples (Table 1), thus proving the usefulness of metabarcoding for chironomid diversity assessment. Some of our retrieved species are indicators for high water quality and were previously detected in spring biotopes, such as Acricotopus lucens, Chironomus luridus, Dicrotendipes lobiger, Limnophyes minimus, Limnophyes pentaplastus, Psectrocladius limbatellus, Psectrotanypus varius and Tanytarsus usmaensis (Reiff et al. 2015), suggesting a general good water quality of our study sites. Only seven of those species can be determined based on larval morphology (Table 1). Even though various determination keys for larval and adult chironomids exist, not all taxa can be determined to species level even by experts (Kranzfelder et al. 2016). Especially larvae and female midges are almost impossible to determine morphologically, since often male genitals are necessary to distinguish species. Some chironomid species have a parthenogenetic life cycle (e.g. Paratanytarsus grimmii, Langton 1988), so only females occur especially in temporary wetlands (Dettinger-Klemm 2003). Without determination of all occurring chironomids, including females, around 27% of the species diversity could be lost (Ekrem et al. 2010). Even if morphology enables the determination down to genus level, an ecological interpretation is difficult since chironomid species of the same genus might have very dissimilar ecological preferences Milošević et al. (2013).

The saprobic index per site based on 14 chironomid species and their sequence abundancies ranged from 1.3 (control site) to 3.3 in one of the Bti treated sites (Table 2). Due to the demanding morphological chironomid species determination, it is common practise to exclude chironomids from bioassessment programs (Milošević et al. 2013). However, in standard water quality assessments in Germany sometimes only all red chironomid larvae are counted, summarized as Chironomus spec. and included in the saprobic index with a value around 3.5. Considering the high chironomid diversity and the range of saprobic values for chironomids between 0.8 (very good water quality) and 3.5 (bad water quality), the standard water quality assessment would have resulted in a severe underestimation of the studied water bodies due to the presence of C. riparius. In addition to the difficult morphological determination, small chironomid larvae (< 1 mm) from freshly hatched species can be easily overlooked by larvae picking. Thus, a metabarcoding approach based on water and homogenized sediment samples could be highly useful for future application in water quality assessments by increasing the chironomid diversity in a sample without specialised taxonomic expertise needed (Bista et al. 2017).

The advantages of metabarcoding over traditional monitoring for water quality assessments is gaining increasing attention. Since Haase et al. 2010) postulated the overlooking of many taxa in traditional stream monitoring programs, many studies proved that metabarcoding can provide higher numbers and more accurate taxonomic identifications than morphology-based methods for many freshwater macroinvertebrates (Hajibabaei et al. 2011, Carew et al. 2013, Elbrecht and Leese 2015, Elbrecht et al. 2017). Moreover, barcoding has increased for rapid biodiversity assessment and biomonitoring for many terrestrial taxa (Yu et al. 2012, Taberlet et al. 2012, Brehm et al. 2013). Ji et al. 2013) compared metabarcoded samples of arthropods and birds with standard biodiversity data sets, and found that the genetic data sets were taxonomically more comprehensive, quicker to produce and less reliant on taxonomic expertise. Cristescu (2014) raised the urge for a coordinated progression of species barcoding that integrates taxonomic expertise and genetic data. For the family Chironomidae an extended and reliably curated barcode database (analogous to the Trichoptera Barcode of Life Database, Zhou et al. 2016) would be highly useful for integrating chironomids in standard freshwater biomonitoring which enhance water quality assessments and might lead to better management of aquatic ecosystems.

Bti effects on chironomid community composition

In our study sites we could show that a considerable number of chironomids live in these wetlands subjected to mosquito control. The uptake and the mode of action of Bti is similar for mosquitos and chironomids (Ali et al. 1981). Regarding potential effects of mosquito control actions using the biocide Bti we expected an overall reduction in chironomid abundance in the Bti-treated sites as well as a reduction in species richness and resulting community composition changes.

The chironomid abundance until WAA 4 was significantly reduced by almost 65% in the Bti-treated sites compared to the Bti-untreated sites, including the never Bti-treated control site. At the control site this spring peak was especially pronounced, indicating that WAA 4 (here: begin of May) is a key time period for the overall chironomid emergence in this area. The observed abundance reduction can be explained by the recent Bti treatment, which killed not only the mosquito larvae but also the chironomid larvae, predominantly affecting freshly hatched larvae. Especially first instar larvae of C. riparius were shown to be highly affected by Bti in laboratory experiments while older larvae were less sensitive (Kästel et al. 2017). Until WAA 13 there was a non-significant trend towards more chironomids in the Bti-treated sites, which could be due to a reduced mosquito competition (Lundström et al. 2010) and subsequently chironomids with a second reproductive cycle in the same year had better conditions (more food resources available) to reproduce. Moreover, species have different egg laying and hatching times, and Bti does not affect eggs but only hatched individuals (Boisvert and Boisvert 2000). The species richness, however, was not significantly different between Bti-treated and untreated sites, neither on OTU-level, accounting for potential cryptic species diversity, nor on species level. The Bti-treatment thus seems to have a mainly quantitative effect on the abundance of the species present in the communities, which is stronger shortly after application during the main chironomid emergence peak in spring.

The adonis analysis corroborated the assumption that site and time (seasonality) had a dominant effect on chironomid communities (Table 3). The predominant effect of site can be explained by the study sites different vegetation (grassland, alder carr, oak carr and pine forest) which influences species composition by varying substrate availability, water chemistry, and the availability of nutritional resources (Van Den Brink and Van Der Velde 1991). Thus, it is not feasible to directly compare the chironomid species composition on the never Bti-treated control site (CL) and the Bti-treated sites (S, G and M) among each other since the vegetational surroundings are quite different and so is the species composition (see Table 2). Note, however, that it is hardly possible to find a “true” control site in the Upper Rhine Valley, i.e., a wetland which has never been treated with Bti next to wetlands subject to mosquito control. In this study we have therefore compared four different study sites, three of which have been subject to mosquito control with Bti.

Even though site and time influenced species composition the most, the first year of Bti intermittence significantly altered the chironomid community as well. This Bti effect was rather low (12%, Table 3). However, considering the more than 20 years of continued Bti application in the study area each spring (in some years even several applications per season) and the proven toxic effect of Bti on chironomid first instar larvae (Kästel et al. 2017), we can assume a more or less depleted community in terms of chironomid diversity. As dispersal for adult chironomids from the next Bti-untreated areas might take too long given the reduced flight capacities (Armitage et al. 1995), a recolonization of univoltine species would probably need longer than one season intermitting Bti treatment. Therefore, resilience in terms of significantly increased species richness may even only be expected after several seasons intermitting the Bti treatment. This highlights the importance for follow-up studies at the sites.

We showed the effectiveness of metabarcoding for chironomid diversity assessments, which led to a 70% increase in species determination compared to determination based on larval morphology. Thus, metabarcoding improves data quality by generating taxonomic resolution. Regarding the question of non-target effects of Bti on the chironomid community, our study found only minor significant effects even though Bti reduced the chironomid emergence by 65%. This could be due to a time lag of chironomid recolonization, since the study year was the first year of Bti intermittence after around 20 years of Bti application in the study area. A follow-up study after a few years of Bti intermittence could result in a more obvious recovery of the chironomid community composition in the Bti-untreated temporary wetlands.

This work has been financed by the Ministerium für Wissenschaft, Weiterbildung und Kultur Rheinland-Pfalz, Germany, in the frame of the programme "Research initiative", project AufLand.We are grateful to the representatives of the town Neustadt an der Weinstrasse and to the Stiftung Natur und Umwelt Rheinland-Pfalz for supporting this project ideally and financially. FL and VE are supported by the EU COST (European Cooperation in Science and Technology) Action DNAqua-Net (CA15219). We thank Martin Spies, Zoologische Staatssammlung München, for advices with questionable species and Christoph Leeb for generating the map of the study sites. Furthermore, we would like to thank Torbjørn Ekrem, Iliana Bista and Kristy Deiner for valuable comments substantially improving this manuscript.

This work has been financed by the Ministerium für Wissenschaft, Weiterbildung und Kultur Rheinland-Pfalz, Germany, in the frame of the programme "Research initiative", project AufLand.

Sampling: AK, SA; Laboratory work: AK, VE; Bioinformatic data analyses: VE; Biological data analyses: KT, AK, JM, SM; Multivariate statistics: SIS; Study design and supervision: KT; CB; FL; Manuscript writing: KT; AK, equal contribution.

All authors edited and commented on the manuscript draft.

The authors declare that there is no conflict of interest.