Reproductive Isolation and Ecological Niche Partition among Larvae of the Morphologically Cryptic Sister Species Chironomus riparius and C. piger

Background One of the central issues in ecology is the question what allows sympatric occurrence of closely related species in the same general area? The non-biting midges Chironomus riparius and C. piger, interbreeding in the laboratory, have been shown to coexist frequently despite of their close relatedness, similar ecology and high morphological similarity. Methodology/Principal Findings In order to investigate factors shaping niche partitioning of these cryptic sister species, we explored the actual degree of reproductive isolation in the field. Congruent results from nuclear microsatellite and mitochondrial haplotype analyses indicated complete absence of interspecific gene-flow. Autocorrelation analysis showed a non-random spatial distribution of the two species. Though not dispersal limited at the scale of the study area, the sister species occurred less often than expected at the same site, indicating past or present competition. Correlation and multiple regression analyses suggested the repartition of the available habitat along water chemistry gradients (nitrite, conductivity, CaCO3), ultimately governed by differences in summer precipitation regime. Conclusions We show that these morphologically cryptic sister species partition their niches due to a certain degree of ecological distinctness and total reproductive isolation in the field. The coexistence of these species provides a suitable model system for the investigation of factors shaping the distribution of closely related, cryptic species.


Introduction
Competition for resources will generally be most severe among closely related species, because they tend to have, due to their shared phylogenetic history, the most similar demands [1,2]. It is widely assumed that the sympatric coexistence of sibling or sister species requires some sort of resource partitioning under resourcelimited conditions [3][4][5], but see [6]. This ''limiting similarity'' concept [7] may not hold under certain, narrowly defined circumstances [8,9], but these instances are believed to be rather the exception from the rule [8] . Hence, testing for differences in the realised ecological niche will consequently be the first logical step in order to explain the coexistence of similar, closely related species. However, closely related species often tend to be morphologically similar for the same reason they are ecologically alike [10]. Therefore, proper species delimitation and unequivocal recognition in field studies are a necessary prerequisite, often requiring molecular methods [11] .
The dipteran midges Chironomus riparius Meigen 1804 (synonym C. thummi, respectively C. thummi thummi) and Chironomus piger Strenzke 1959 (synonym C. thummi piger) are sister taxa [12,13]. Larvae of both species are widely distributed in small streams, ditches, ponds and puddles throughout the holarctic [14]. The life cycle of C. riparius and C. piger consists of four larval stages, a short pupal stage and the adult midge. Adults form large mating swarms. A few days after hatching, female midges usually produce a single egg mass containing several hundreds of eggs. The larvae hatch after a few days, and the whole life cycle may be completed within four weeks. Depending on the water temperature, both species are usually multivoltin, with a first generation emerging early in spring and the final generation swarming around late autumn [15]. The overwintering generation consists solely of later larval stages (L3, L4). The species are often dominating the local Chironomus community [16]. In the study region they are frequently found together at the same sites [14], making the species pair an interesting model for the investigation of mechanisms enabling the sympatric coexistence of sibling species.
As the two sister taxa are morphologically cryptic, safe species discrimination was only possible only by analysis of polytene chromosomal structure in the past [17]. Despite their morphological similarity, genome size differs by 30%, mainly due to repetitive DNA [18]. Not only their taxonomic status regarding species or subspecies rank is unclear, also the reports on their degree of reproductive isolation are inconsistent. Some degree of prezygotic isolation in the field is warranted by differential swarming behaviour [19]. While some studies indicate that C. riparius and C. piger readily form viable and fully fertile interspecific hybrids in the laboratory [20], others estimate fertile hybrids in the wild to be effectively absent, due to fertility reductions caused by hybrid dysgenesis syndromes [21]. The actual degree of hybridisation and reproductive isolation in the field, however, has not yet been explored.
In this study we aimed to investigate distributional patterns of both species in an area where both species co-occur, and to reveal ecological factors that may have shaped the observed distribution. To this end, we investigated genetic differentiation between the species using mitochondrial and nuclear markers and related their relative abundance with environmental parameters. In particular, we answered successfully the following questions: -What is the degree of reproductive isolation among C. piger and C. riparius in the field, -is there a non-random spatial pattern of distribution and cooccurrence, and -can we identify ecological parameters potentially structuring the species distribution?

Species delimitation and identification
Two hundred sixty four individuals of C. riparius/piger were found at 34 sampling sites (Table 1). Microsatellite analysis detected a total of 76 alleles at the five loci (mean = 15.2, s.d. = 8.1). Factorial correspondence analysis on the microsatellite data revealed two distinct genotype clusters, termed A and B ( Figure 1A). Their distinctness was due to both private alleles and frequency differences at all loci ( Figure 1B). Identical results were obtained with other assignment methods like STRUCTURE (Pritchard et al., 1999) (results not shown). The statistical parsimony network revealed two major haplotype groups, linked by six mutational steps. Plotting the two nuclear genotypes on the haplotypes of the respective individuals revealed a complete congruence with these two haplogroups ( Figure 2). Polytene chromosome preparations identified genotype A (black symbols) consistently as C. riparius and genotype B (grey symbols) as C. piger.

Co-occurrence and population structure
At about half of the sampling sites containing C. riparius or C. piger, both species co-occurred in varying proportions ( Figure 3). However, individuals of the species occurred less often together than expected by chance (Fisher's exact test x 2 = 160, d.f. = 22, p,0.0001). The relative frequency of a species at a given site was not independent from their frequency at surrounding sites. We found a significant spatial autocorrelation of sampling sites up to 15 km apart ( Figure 4). A significant, albeit very weak genetic population structure within both species was detected. For C. piger, a W ST of 0.027 was calculated, while the estimate for C. riparius was 0.046 (Table 2), indicating a high amount of gene-flow among sampling sites among individuals of each species.

Discussion
Chironomus riparius and C. piger behave as good species in the field Microsatellite analysis showed the presence of two distinct genotype groups without intermediates, indicating complete reproductive isolation and the absence of putative hybrids ( Figure 1A). Most alleles were specific to one of the cluster with only few alleles shared in similar proportions by the two taxa ( Figure 1B). As only the lengths of PCR-fragments were scored, this partial overlap may be due to homoplasy or common ancestry. The results were so clear cut that the application of more sophisticated methods of hybrid detection (e.g. NewHybrids) was deemed unnecessary. The inference of reproductively isolated gene-pools is strengthened by the distinctness of the mitochondrial variation of the genotype groups ( Figure 2), indicating long lasting isolation with absence of both current and past hybridisation [22]. Even though the generation of hybrids in the laboratory is possible to varying degrees [17,21], both pre-and postzygotic isolation mechanisms [19,21] seemed to have maintained complete reproductive isolation in the wild, despite the opportunity to interbred. Therefore, the two taxa conform to several species concepts, including the biological [10] at least in the area investigated and they should be consequently regarded as good species.
Spatial repartition of C. riparius and C. piger along ecological gradients The unequivocal species assignment by molecular markers showed that larvae of C. riparius and C. piger occurred not only in the same general area, but in about half of the cases at the same site ( Figure 3). Still, they were found less often syntopically than expected from their overall abundance. The virtual absence of population structure on the spatial scale of the study (Table 2) suggests that this is not due to dispersal restrictions or geographical obstacles. It indicates rather competitive interaction, either present or past [23]. C. piger was dominant mainly in the west of the area, while C. riparius occurred more frequently in the east (Figure 3), as mirrored in the significant spatial autocorrelation of relative species frequency ( Figure 4). This pattern corresponds to the correlation of the species' relative frequencies with parameters measuring the average amount of precipitation during the summer months ( Table 3). The latter are climatic parameters which generally tend to be spatially autocorrelated. Other parameters that covaried significantly with relative species abundances were water chemistry variables (conductivity, nitrite and CaCO 3 , Table 3). Multiple regression retained only the precipitation variables (Table 4). This suggests that differential desiccation resistance, as observed in other Chironomus species [24], could have caused the observed patterns. However, the absolute differences in precipitation (see Appendix) are probably too small for a differential desiccation risk of the water bodies in the area. Therefore, the climate constitutes probably merely the ultimate cause for the differential species distribution. The amount of rain during the warm summer months with their increased evaporation determines the concentrations of nitrite and other ions in the shallow puddles and ditches both species inhabit. It has been shown that both high salinity and high nitrite concentrations impart larval development in C. riparius/piger [25,26]. Our results indicate that C. piger occurs in areas with less summer rain and tolerates higher nitrite concentrations and conductivity than C. riparius (Table 3). Therefore, the proximate cause for the observed correlation of the species frequencies with summer precipitation is more likely the gradient of water chemistry variables during the time of highest larval abundance [27].
Even though the inferred spatial repartition along ecological gradients is rather a hypothesis in need to be confirmed by subsequent experiments in the laboratory, it has become evident that C. riparius and C. piger are ecologically not completely equivalent. Although this has already been suspected before [14], our study is the first to demonstrate ecological partitioning among the species pair quantitatively in the field. Studies on the ecological differentiation of other Chironomus species have revealed a range of mechanisms that structure coexistence in sympatry. Dietary niche separation among two profundal species from the Chironomus plumosus-group has been suggested by stable isotope analysis [28]. The same species were found to differ in emergence time, suggesting also temporal niche separation [29]. Perhaps the most impressive example of interspecies competition avoidance is the spatial repartition of temporary rain water puddles by C. pulcher and C. imicola into shaded and sunny regions on a very small scale [30].
Despite the demonstrated spatial repartition along ecological gradients of C. riparius and C. piger, we found a substantial number of sites where both species co-occurred, indicating a substantial overlap in the realised ecological niche. Possible, not mutually exclusive explanations for this pattern include: i) a substantial stochasticity in the dispersal/colonisation of the sites. Even though the oviposition choice in another Chironomus species is influenced by nitrogenous compounds and conspecific larvae [31], a high degree of randomness regarding environmental conditions is generally assumed in the community assembly of chironomids [27]. Also which species arrived first at a yet unoccupied site may crucially influence the outcome of subsequent competition [32]. ii) Temporally fluctuating environmental conditions may also prevent complete competitive exclusion [33]. iii) Interaction with other species. Several other species of Chironomus are present at most of the investigated sites [16], as well as other mud dwelling taxa with similar requirements. iv) the abundance in the neighbourhood possibly also influences the local abundance of C. riparius and/or C. piger [33]. As this study documents, C. riparius and C. piger provide a promising model for the investigation of factors shaping the distribution of closely related, cryptic species. Currently ongoing experimental and ecological genomic studies on this emerging model system will help to gain a deeper understanding of the processes and factors that shape the realised niche of closely related species in sympatry. Understanding the internal factors and constraints shaping their distribution and coexistence will contribute to our mechanistical understanding of the processes shaping biodiversity in ecological communities.

Sampling
The sampling area lies in the middle of the upper Rhine valley in a rectangle of roughly 40 by 60 km between 49u099-49u339N and 8u109-8u139E. It comprises the Rhine valley plain, in the west limited by the mountains of the Pfä lzer Wald and in the east by the rising hills of the Odenwald range. The area is hydrologically characterised by the presence of many drainage ditches, slowly flowing small streams, temporary puddles, the oxbows and the main stream of the river Rhine.
The sampling took place from mid September to November 2004, thus sampling the over-wintering generation of Chironomus larvae [34]. The sampling period was scheduled in autumn in order to avoid the large fluctuations in abundance among species throughout summer. Moreover, sampling the hibernating larvae assemblage that will foster next years first generation presents the result of competition processes during the growth season [35].
Sampling took place as described in [12]. Briefly, potential Chironomus habitats were considered opportunistically within the study region, but we mainly focused on typical Chironomus riparius/  piger habitats (small streams, creeks, and ditches with fine, muddy sediment). An area of 161 m was sampled with a 30640 cm net of 0.5 mm mesh size. Due to the small size and low depth of most water bodies, we did not consider different areas within a water body during sampling. All Chironomus larvae instar stage found (instar stage 3 and 4), as identified by the presence of ventral tubuli, were brought alive into the laboratory. For the present study, we chose all thirtyfour sampling sites where C. riparius and/or C. piger which had been identified earlier using a COI barcoding approach [16].

DNA isolation and microsatellite analyses
Larvae were kept in the laboratory for at least 5 days without feeding, in order to remove potential PCR inhibiting substances from the gut [36]. Head and first body segments were removed for polytene chromosome analysis as described in [17]. Briefly, salivary glands were prepared from fresh larval tissue and fixed in 50% acetic acid. Chromosomes were stained in 2% orcein acetic acid for 15 min and fixed on glass slides for microscopical analysis. Remaining caudal tissue was homogenized in 700 ml standard CTAB buffer containing 0.1 mg/ml proteinase K. After digestion for at least 1 h at 62u C, chloroform/isoamyl alcohol 24:1 treatment was performed followed by 1 h precipitation at 220u C. DNA pellets were washed twice with ethanol 70% and resolved in 30 ml water.
Genetic structure, mitochondrial haplotype phylogeny and species identification Factorial correspondence analysis (FCA) was applied on multilocus genotypes to explore the distribution of genetic variation graphically (GENETIX 4.04 software, [38]. Genetic population structure was assessed using the AMOVA approach [39] as implemented in the Excel add-in GenAlEx [40]. For this analysis only sampling sites with at least seven conspecific individuals were taken into account. C. riparius/piger COI haplotypes were identified from [16] (GenBank Accession numbers DQ910547-DQ910729). The phylogeny of the COI haplotypes was inferred using statistical parsimony (SP) [41]. The SP network was constructed with the program TCS v. 1.21 [42]. Nesting of clades followed the rules given in [43] and [44]. Inferred reproductively isolated entities were taxonomically identified using polytene chromosome preparations of a subset of individuals [17].

Co-occurrence
We used a Fisher's exact test (10 6 permutations) to investigate whether the co-occurrence of the identified taxa was random. Spatial patterns of the relative frequency of C. riparius and relevant environmental parameters at the sample sites with at least seven individuals found were inferred with spatial autocorrelation analysis. Seven mutually exclusive lag classes of 5000 m width were used to compute Moran's I spatial correlation coefficient for each class. Statistical significance of Moran's I was assessed with 999 Monte Carlo permutations. The Excel Add-in RookCase version 0.99 [45] was used for the calculations.

Physico-chemical and climatic characterisation of sampling sites
Thirty-eight ecological parameters were recorded in order to characterize abiotic habitat conditions at the respective sampling size. These parameters were chosen to cover a wide range of ecological conditions known to influence freshwater communities, and the distribution of chironomid species in particular. Recorded characteristics include physicochemical parameters [46,47], sediment composition [48], climatic conditions [49], and structural habitat characteristics (e.g., size and depth of water body).  For the determination of sediment organic content, measured as loss on ignition, approximately 30 g of sediment sample were dried at 60uC for three days and weighed subsequently. Samples were then muffled at 550uC for 4 h, followed by determination of percentage weight loss.. For the identification of relative particle size composition of the samples, 150 g of homogenised sediment were washed through six sieves with decreasing mesh size and the content of each sieve was dried and weighted.
Conductivity, pH, water temperature and O 2 saturation were measured with a WTW Multi 340i multimeter at each sampling site. Ammonium, nitrite and phosphate concentrations were calorimetrically determined using AquamerkH quicktests. Chloride, CaCO 3 and nitrate concentrations were measured with colour tests (MerkoquantH). The stream velocity was measured using an AMR ALMEMOH device.
Nineteen biologically meaningful climatic parameters were extracted for each sampling site from the BIOCLIM environmental layers with a spatial resolution of 0.5 min, implemented in the computer program DIVA-GIS version 4.2 [50]. Mean, median, standard deviation, minimum and maximum values for the recorded parameters are given in Appendix S1.

Statistical analysis
Means, standard deviations, median, minimum and maximum values for all 38 variables taken into account are given in the Appendix. All data with the exception of pH were either log 10 (x+1; continuous variables) or arcsin (percentages) transformed to conform to the underlying assumptions of normality and heteroscedasticity in subsequent analyses. We calculated Pearson's correlation coefficients (r) between the relative C. riparius frequencies and all respective variables. Due to the multitude of comparisons, we calculated a q value for each test to estimate the minimum false discovery rate which is incurred when calling that test significant. Variables with p values,0.05 and q values,0.10 in correlation analysis were retained for a multiple regression (forward selection).

Supporting Information
Appendix S1 Recorded environmental parameters.