Chironomid species diversity at extremely high altitudes

Chironomid species diversity at extremely high altitudes

July 9, 2013

a Institute of Biodiversity and Ecosystem Dynamics (IBED), FNWI, University of Amsterdam. Science Park 904, 1098 XH, Amsterdam, the Netherlands

*Phone: +316 15135905; e-mail: roland.vandierendonck@student.uva.nl

Abstract

The high mountain ridges of the Andes may function as a geographical barrier, creating the conditions for allopatric speciation of aquatic life in high altitude streams. The aim of the present study was to compare chironomid species diversity in two separate streams. It was hypothesized that metal polluted waters downstream and mountain ridges separating the streams limit dispersal. Individual larvae were collected at 3000 m and 4000 m in the Quilcayhuanca catchment and at 4000 m in the Rúrec catchment, both in the Cordillera Blanca (Peru). The genetic structure of the chironomid sample was determined by sequencing the mitochondrial COI gene and the nuclear CAD gene and constructing Maximum Likelihood phylogenetic trees. Chironomid composition of the Quilcayhuanca 4000 m and Rúrec site were highly differentiated, but the sites had one phylogenetic species in common. Phylogenetic species composition of sites at 3000 and 4000 m in the Quilcayhuanca stream differed completely. Interestingly, the Rúrec site shared phylogenetic species with the 3000 m Quilcayhuanca site. It is concluded that dispersal of adult chironomids between the two streams is taking place. Nevertheless, the fact that the Quicayhuanca and Rúrec stream showed a different chironomid species composition could indicate that mountains are indeed a barrier for dispersal.

2

Introduction

Ever since the publication of Darwin’s On the Origin of Species (1859) the mechanisms that lead to the divergence of species have been topic of discussion. Dobzhansky (1935) reasoned that in the majority of cases, speciation is preceded by reproductive isolation. Over time, individuals from non-mixing populations can become incapable of interbreeding, the different populations consequently being classified as separate species. Contemplating on the cause of reproductive isolation, Dobzhansky (1937) proposed a variety of “isolating mechanisms” including geographical (physical distance), ecological (different habitat types), sexual (mating choice), mechanical (disparate genitalia) and genetic (hybrid sterility) isolation. The geography of speciation has been studied intensively and a range of speciation scenarios has been envisaged (Coyne & Orr, 2004). At one side of the spectrum there is sympatric speciation, or the divergence of species within a single geographic area. At the other side is allopatric speciation, which comprises geographical barriers that separate populations and thereby physically inhibit gene flow (Mayr, 1963). Allopatric speciation has long been the dominant explanatory principle in evolutionary biology, but various recent studies using a phylogenetic approach have showed that sympatric speciation is more common than was previously assumed (Losos & Glor, 2003; Johannesson, 2010). However, when related species are separated by physical barriers, such as seas, rivers or mountain ridges, allopatric speciation would seem to be the most probable interpretation (Losos & Glor, 2003).

One mountain range that is surely high enough to provide physical barriers is the Andes. Uplifted circa 12 million years ago, the Andes is a mountain range situated along the west coast of South America, with an altitude averaging 4000 m above sea level (Hoorn et al., 2010). In the Cordillera Blanca, part of the Peruvian Andes, glacial streams are geographically separated by high mountain ridges. Aquatic life in these high altitude streams is challenged by harsh conditions, including low water temperature, low oxygen saturation, low nutrient availability, high current velocity, high water discharge, and high ultraviolet-B radiation (Jacobsen et al., 2003; Jacobsen, 2008; Loayza-Muro et al., 2010). The latter is partly explained by a thinner ozone layer in tropical zones, a relatively short atmospheric sunlight path at higher altitudes, and increased exposure due to the transparency of freshwater streams with a low dissolved organic matter concentration (Villafañe et al., 1999; Kelly et al., 2001). In addition, biota in high alitude Andean streams is influenced by increased metal concentration, due to natural leaching by rock weathering or human mining activities (Loayza-Muro et al., 2010). Adverse effects of elevated metal

3 levels include metal toxicity and increased acidity. Furthermore, the smothering of the streambed with layers of metal precipitates reduces habitat and food availability and quality for benthic fauna (Courtney & Clements, 2002).

The composition of macroinvertrebrate communities in Peruvian high altitude streams is affected both by altitude and metal related factors. Loayza-Muro et al. (2010) observed that the macroinvertrebrate community in unpolluted streams consists of a relatively large percentage of mayflies, stoneflies and caddisflies, whereas in polluted streams, dipterans and coleopterans were abundant. One taxon notable for its resilience to pollution is the dipteran family Chironomidae. The benthic larvae of these insects are among the most abundant life forms in Andean freshwater streams. Comparing polluted with reference high altitude streams at 3000 m (‘low altitude’) and 4000 m (‘high altitude’) in the Peruvian Andes, Loayza-Muro et al. (submitted) found that polluted streams both at low and high altitude harboured only one single chironomid species, whereas unpolluted streams displayed greater species richness. Furthermore, the chironomid species living in polluted streams appeared at both high and low altitude, indicating that metal concentration is a stronger predictor of species diversity than altitude. The 3000 m and 4000 m unpolluted streams, however, had no species in common.

High mountain peaks could form a barrier for adult chironomids to disperse between high altitude Andean streams. Furthermore, metal pollution downstream could pose a challenge to the survival of drifting larvae. The question whether these factors prevent gene flow between invertebrate communities in separate streams, thereby creating conditions for allopatric speciation, remains unanswered. A tree displaying the evolutionary relationships of individual larvae could aid inferences on past speciation events. In a follow-up to De Baat and Loayza-Muro’s study in 2012, the present study compares the chironomid species diversity in the Quilcayhuanca and Rúrec catchment, two geographically separated unpolluted freshwater streams in the Peruvian Andes. If migration and movement via water is not possible and if mountains prevent adult dispersal, a different species diversity would be expected at both sites.

Although the morphological characters of chironomids have been thoroughly studied, it remains difficult to determine individual chironomid larvae to the species level using these traits – in fact, chironomid larvae display little morphological diversity within subfamilies and genera. Furthermore, there is no taxonomic key for Andean Chironomidae, as past studies primarily focused on North-American and European chironomids (Cranston et al., 1995). Therefore, in this study, genetic barcoding is used as an alternative way to determine species richness (Hebert et al., 2003). By comparing genetic sequences of the mitochondrial COI gene and the nuclear CAD-complex gene, individuals can be

4 classified as different genotypes, adopting a phylogenetic species concept that takes the relative difference in genetic code as a measure of relationship (Coyne & Orr, 2004). Using this data, a phylogenetic tree was build. Given that the Quilcayhuanca and Rúrec streams, though only 35 km apart, are being separated by high mountains possibly preventing dispersal of adults, and that metal polluted waters downstream possibly counteract dispersal of larvae, it was expected that the sites have few phylogenetic species in common.

Figure 1: Map of the area. Chironomid

larvae used in this study were collected at four sites in the Cordillera Blanca (Peru), here shown as black dots on the map. The streams are described as Quilcuayhanca Reference High (QRH), Quilcayhuanca Reference Low (QRL), Quilcayhuanca Polluted High (QPH) and Rúrec Reference High (RRH), where High is 4000 m and Low is 3000 m. Triangles represent the mountain range separating the two catchments. The study sites are located near the city of Ancash (below). Figure is modified from Loayza-Muro et al. (2013).

5 Materials and Methods

Sample sites

The Cordillera Negra and Cordillera Blanca are two parallel mountain ranges in the mid-western part of the Peruvian Andes, situated at opposite sides along the Santa river. In the Cordillera Blanca, the most eastern of the two mountain ranges, samples were taken at 4000 m in freshwater streams in both the Quilcayhuanca catchment and the Rúrec catchment, about 35 km apart. Additionally, samples were collected at a polluted site at 4000 m and a 3000 m unpolluted site, both in the Quilcayhuanca catchment. All sites are situated in the Huaraz area (Figure 1).

Chironomid sampling

Individual chironomid larvae were collected by Loayza-Muro and De Baat in January and February 2012, using sieves and plastic trays to attain them from gravel-pebble sediments and stones on the stream banks. Individual larvae were separated in 1.5 mL screw cap tubes, and 100 µL of TRIzol (Invitrogen, USA) was added. Samples were stored at -80 °C, shipped to the University of Amsterdam (The Netherlands), and stored at -80 °C. A total number of 58 chironomid larvae – 16 from the 4000 m Quilcayhuanca reference site and 42 from the 4000 m Rúrec reference site – were selected for DNA isolation and amplification.

DNA isolation

For DNA isolation from larvae in TRIzol, a protocol was followed based on that of Loayza-Muro et al. (submitted). Samples were taken from -80 °C and thawed at room temperature (RT). Subsequently, 20 µL of chloroform was added and samples were incubated at RT for 2-3 minutes. Next, samples were centrifuged at 12000 x g for 15 minutes at RT (samples 14-23 and 294-317) or 4 °C (samples 148-153 and 318-335). The aqueous phase, containing RNA, was removed and stored at -80 °C. To the remaining organic phase and interphase, 30 µL of pre-warmed TNES-6U buffer (10 mM Tris; pH 7.5, 125 mM NaCL, 10 mM EDTA-2Na; pH 7.5, 1 % SDS and 6 M Urea) was added. After 10 minutes of incubation at RT, samples were spun at 18000 x g for 15 minutes at RT (samples 14-23, 294-317 ) or 4 °C (samples 148-153, 318-335). The aqueous phase, containing DNA, was transferred to a clean tube, whereas the interphase, rich in protein, was discarded. An equal volume of ice cold 2-propanol was added to the aqueous phase and samples were incubated at -80 °C overnight. Samples were spun at 18,000 x g for 30 minutes at 4 °C. The supernatant was discarded, the pellet washed three times with 100 µL ice-cold 80 % ethanol and dried in a vacuum centrifuge. DNA was redissolved by adding 20 µL of T10E1 buffer (10 mM

6 Tris, 1 mM EDTA) to the tubes, which were placed in a 37 °C heating block for 10 minutes and vigorously shaken. Finally, DNA concentrations were measured using a spectrophotometer (NanoDrop, USA).

DNA amplification and sequencing

A ±700 base pair (bp) fragment of the mitochondrial COI gene and a ±1000 bp fragment of the nuclear CAD gene were amplified. Primers 911 (forward) and 912 (reverse) and 54 F (forward) and 405 R (reverse) were used to amplify gene fragments of respectively COI and CAD (table 1).

Table 1: Primers used in PCR amplification and sequencing

Locus Name Length (nt) Sequence Reference

COI 911 25 TTT CTA CAA ATC ATA AAG ATA TTG G (Folmer et al., 1994)

912 26 TAA ACT TCA GGG TGA CCA AAA AAT CA (Folmer et al., 1994) CAD 54 F 23 GTN GTN TTY CAR ACN GGN ATG GT (Moulton &

Wiegmann, 2004) 405 R 23 GCN GTR TGY TCN GGR TGR AAY TG (Moulton &

Wiegmann, 2004) R= G A (Purine)

Y = C T (Pyrimidine)

For DNA amplification, to each tube a 20 µL reaction volume consisting of 8.2 µL H2O, 4.0 µL 5 X PCR buffer, 4.0 µL 1 mM dNTP’s, 0.6 µL 10 mg/ml bovine serum albumin, 0.4 µL 10 µM primers (forward and reverse), 0.4 µL 5U/µL “Phyre Hot Start II” polymerase (Finnzymes, Finland) and 2.0 µL template DNA were added.

For COI amplification, an initial denaturation step at 98 °C for 10 s was followed by an annealing step at 55 °C for 30 sec. A cycle of denaturation at 98 °C for 10 s, annealing at 55 °C for 10 s and elongation at 72 °C for 20 s was repeated 35 times. The PCR was concluded by an elongation step at 72 °C for 5 min and a cooling step at 4 °C for 5 min.

For CAD amplification, a 3 min denaturation step at 98 °C was followed by 4 cycles of denaturation at 94 °C for 30 s, annealing at 52 °C for 30 s and one min elongation at 68 °C; 6 cycles of denaturation at 94 °C for 30 s, annealing at 51 °C for one min and elongation at 68 °C for one min; 36 cycles of denaturation at 94 °C for 30 seconds, annealing at 45 °C (samples 16, 17, 19, 20, 295 and 297)

7 or 49 °C (other samples) for 20 seconds and one and a half minute of elongation at 68 °C; 10 minutes of extension at 68 °C and 5 minutes of cooling at 4 °C.

Finally, 3 µL of PCR products was loaded on an 1 % agarose gel in order to control whether primer dimers had formed and whether the amplified gene fragments were of the required size. CAD PCR products occasionally showed multiple bands on the gel, indicating contamination or non-specific amplification. In these cases, the PCR product was purified. The remaining PCR product was run on an 1 % agarose gel. Bands were cut out under UV-light and cleaned following the Invisorb Spin DNA Extraction kit protocol. PCR products were send for sequencing to MacroGen (Amsterdam).

Of the COI PCR products, 48 were sequenced, whereas 40 CAD PCR products were sequenced. COI sequences were generally of high quality, but the majority of CAD sequences was not included in the final dataset because of small contig size and large stretches of poor quality.

Data analysis

A BLASTn search in GenBank was executed to confirm that sequences were indeed fragments of COI or CAD and of chironomids. Furthermore, sequences were translated into amino acids to ensure no stop codons were present in the amplified DNA fragments, using either an invertebrate mitochondrial DNA translation table or a general nuclear DNA translation table. Sequences were analysed and modified in CodonCode Aligner (CodonCode Inc.). Reviewing the chromatogram data in the Trace window, the 5’ and 3’ 20-200 bases of each sequence were deleted because of poor sequencing quality. Forward and reverse sequences were assembled in contigs. Ambiguous positions in the CAD sequences – possibly a display of heterozygosity at these bases – were marked with nucleotide ambiguity codes (IUPAC). High quality sequences were imported in MEGA 5.2.1, where they were aligned using the multiple sequence alignment programme ClustalW and trimmed to equal lengths (Tamura et al.,2011). A model test was run to find the best DNA substitution models, rates and patterns. Using the Maximum Likelihood method, trees were constructed using MEGA software, applying 500 bootstrap replications to test the robustness of the phylogeny reconstruction. Gaps and missing data were completely removed from the analysis. The heuristic method of nearest-neighbour-interchange was applied to infer trees. Initial trees were made automatically using the default NJ/BIOJN algorithms, including 1st+2nd+3rd codon positions. A tree based on COI sequences and a tree based on CAD sequences were constructed. In order to root the trees, the first hits when BLASTing one of the sequences (89% maximum-identity) were downloaded from GenBank.

8 represented in respectively the COI tree and CAD tree. Furthermore, 53 COI sequences and 28 CAD sequences of chironomid DNA previously isolated, amplified and sequenced by Milo de Baat were added to create a total dataset of 89 COI sequences and 44 CAD sequences. These included DNA fragments of nine individuals collected in an 3000 m unpolluted stream and three individuals collected in a polluted stream at 4000 m, both in the Quilcayhuanca area.

Results

As shown in figures 2 and 3, both the COI and CAD sequence Maximum Likelihood tree are well-resolved, with only the tree based on COI fragments displaying a single polytomy. Furthermore, both trees have high bootstrap values at the final clades, indicating evidence for the shown clustering of individuals is solid. In figure 3, monophyletic clusters in the CAD tree are defined by small roman numerals. Corresponding clusters in the COI tree in figure 2 are noted by the same numerals. Incongruities between the mitochondrial COI and the nuclear CAD tree could indicate a case of hybridization. However, since for all individuals embedded in both the COI and the CAD tree, the monophyletic groups matched, the likelihood of monophyletic groups in the tree representing real species increased. Because there is little variation within clades, it could be assumed that each branch represents a distinct haplotype. The number of individuals in a cluster vary from 1 to 35. Since bootstrap values at earlier branchings are on average low, ancestral reconstructions are weak.

Interpreting the clusters in the phylogenetic reconstructions, a 2 % species concept was adopted, meaning that a 2 % difference between codes defines individuals as being of different phylogenetic species (Nixon & Wheeler, 1990). Following this reasoning, as in both trees branches are longer than the 2 % distance, all colored boxes represent phylogenetic species (figures 2 and 3). Some individuals are represented in the CAD tree but not in the COI tree and vice versa. Because the number of usable COI sequences was higher than the number of CAD sequences, it was decided to use the genetic structure of COI sequences for further analysis. Relative abundances of phylogenetic species at the different sites were visualized in pie charts (figure 4).

At the Quilcayhuanca 4000 m unpolluted site, four phylogenetic species were distinguished in 36 samples. Thus, one species is added to the inventory of Loayza-Muro et al. (submitted), who identified three separate species at this site. The majority of individuals (86 %) clustered in one monophyletic group. Other phylogenetic species appear to be much less abundant, as only a few individuals (3,1 and 1, respectively) belonging to these groups were collected. Consequently, it appears that the observed species diversity is a function of sampling size, as the more individuals one collects, the more species one encounters. The Quilcayhuanca site at 3000 m harbors at least five species in 8 samples, neither of which

9

Figure 2: COI Tree. Maximum

Likelihood tree of partial COI sequences from chironomid larvae collected on different sites in the Quilcayhuanca and Rúrec catchment (legend) with the highest log likelihood (-3476.84). The best fitting substitution model was the General Time Reverse model + Gamma-distributed variation + invariant sites (GTR + G + I). A 544 bp fragment of the COI sequence of 90 individuals – 89 from the study side and one from a Costa Rican sample that was downloaded from GenBank as an outgroup to root the tree. The bootstrap values, percentages of replicates in which the nodes were similarly placed are shown above the bars. Scale = distance, 0.02 corresponding to the 2.0 % species distance. Clusters that correspond with a haplotype in the CAD tree, are marked with the small roman numerals used to distinguish haplotypes in the CAD tree (figure 3).

is also present at the 4000 m Quilcayhuanca site. This was already shown by Loayza-Muro et al. (submitted) (figure 2).

Of 27 individuals collected in the Rúrec 4000 m reference site, 15 % clustered with a large monophyletic group of Quilcayhuanca 4000 m specimen (figure 3). Another individual clustered with two individuals from Quilcayhuanca 3000 m.

10

Figure 3: CAD Tree. Maximum

Likelihood tree of partial CAD sequences from chironomid larvae collected on different sites in the Quilcayhuanca and Rúrec catchment (legend) with the highest log likelihood (-16310.60). The best fitting substitution model was the Tamura-Nei model + Gamma-distributed variation (TN +G) A 785 bp fragment of the CAD sequence of 45 individuals – 44 from the study side and one from an sample from Thailand that was downloaded from GenBank to root the tree. Bootstrap values are shown above the bars. Scale = distance, 0.02 corresponding to the 2.0 % species distance. All eight haplotypes are marked with small roman numerals.

11 Discussion

The unpolluted sites at 4000 m in the Quilcayhuanca and Rúrec catchment are characterized by similar environmental values (Loayza-Muro et al., 2013). However, since the sites are physically separated by a mountain ridge possibly inhibiting dispersal of adult chironomids and since metal pollution downstream prevents larval dispersal, it was expected that the genetic diversity of chironomids at the two sites would be dissimilar. Comparing the phylogenetic species composition at both 4000 m sites, this inference seems valid. The diversity at both 4000 m streams was mainly composed of one very common species, with 80 % of the sampled individuals being classified as belonging to that taxon, and three less abundant species (figure 4). One species was found at both sites, though being far more abundant in the Quilcayhuanca stream than in the Rúrec. This indicates that there has been dispersal of chironomids across the mountain range between Quilcayhuanca and Rúrec catchments at high altitudes. Nevertheless, the high differentiation between sites could indicate that dispersal is rare.

Comparing physical and chemical variables at 3000 m and 4000 m in the Quilcayhuanca catchment, Loayza-Muro et al. (submitted) measured a 50 % increase in UV-B radiation with elevation. Conversely, temperature decreased by almost 50 % from 3000 m to 4000 m. As was previously found, the two unpolluted Quilcayhuanca sites had no species in common. Interestingly, one individual from the Rúrec 4000 m site was identified as belonging to the same phylogenetic species as two individuals from Quilcayhuanca, 4000 m (n=36) Quilcayhuanca , 3000 m (n=8) Rúrec, 4000 m (n=27)

86% 8% 3% 3% 31 3 1 1 25% 25% 25% 12% 13% 2 2 2 1 1 78% 15% 3% 4% 21 4 1 1

Figure 4: Relative abundances. Pie charts of relative abundances of phylogenetic species collected at the

Quilcayhuanca 4000 m and 3000 m and Rúrec 4000 m unpolluted sites. Numbers were derived from the Maximum Likelihood tree of COI sequences. Different colors correspond to different species, relative abundances of species are written inside or outside the slices. At the boxes right next to the charts the absolute number of individuals are shown.

12 the Quilcayhuanca 3000 m site. However, the low number of larvae involved in this particular result may be a reason to withhold from drawing strong conclusions.

Phylogenetic species versus real species

One of the underlying assumptions of the phylogenetic species concept applied in this study is that intraspecific variation is smaller than 2 %. Compared to other insect families, which display highly variant and often overlapping intra- an interspecific divergences, chironomids have relatively stable COI genes (Cognato, 2006). For example, between species within the Tanytarsini tribe in the subfamily Chironominae an average intra- and interspecific distance of respectively 0,87% and 14,7% was measured (Ekrem et al., 2006). Furthermore, because the COI tree and CAD tree were congruent, and variation within clades was observed to be low, the likelihood that phylogenetic species actually correspond to real species is increased.

Individuals could be classified into different phylogenetic species, but it appeared to be impossible in this study to infer to which genus, order or even which subfamily within the Chironomidae these species belong. This is partly explained by the fact that South-American chironomid species are poorly represented in GenBank. When BLASTing CAD sequences, even species from families such Cucilidae or Drosophilidae had a similar distance as did Chironomidae. Another reason why inferring the higher taxonomic levels was undoable is that COI sequences may in general have poor power when it comes to identification of unknown species, because the monophyletic groups in COI trees may not be fully congruent with species trees (Ekrem et al., 2007). Furthermore, in order to sketch scenarios of the evolutionary past of the encountered species, the inner nodes of the tree could potentially depict the relations. However, the history of descent of the fourteen phylogenetic species in the COI tree is vague since the bootstrap values indicating the percentage of bootstrapped trees that have the shown configuration are low. The uplift of the Andes, 12 Ma could not be used to add a time axis to the phylogenetic trees, as it is advised not to rely solely on biogeographical events, but use an external calibration such as fossils of ancestral taxa (Cranston, 2012).

Generalizing samples

When attempting to extrapolate the results from this study, several problems were encountered. First of all, the dataset is limited to a number of 89 individuals. It cannot be excluded that if more individuals are sampled at the exact same locations, a different genetic structure is found. Generally, in every sample of individuals, a small number of species is very abundant and the majority of species is rare, meaning the more individuals are collected, the more species are found (Hubbell, 2001). A second problem is of

13 phonological nature. Individuals were collected in January and February, it could be that the presence and abundance of species correlates with the time of the year.

The first of the two problems can be handled by testing whether the number of individuals is a statistically sufficient number to draw conclusions. It would be very interesting to investigate whether the found species abundances were to be expected based on the sample size and the number of sampling sites, and also whether the species seem randomly distributed across the area, or if there is some statistical evidence for migrational barriers. In order to investigate these questions, a formula for neutral biodiversity could be used (Etienne, 2005). Theories of neutral biodiversity assume that relative species abundances of species at the same trophic level are predominantly explained by random fluctuations and population, as opposed to niche related factors (Hubbell, 2001). However, the number of individuals in this study was too little to give enough power to the formulae. Hence, future studies should collect more larvae so stronger inferences can be made about the genetic diversity of chironomids in Andean high altitude streams. It would also be wise to increase the number of sample sites.

In conclusion, because there are phylogenetic species that occur on both sites of the mountain range, it can be deduced that dispersal of adult chironomids takes place between the Quilcayhuanca and Rúrec stream. Nevertheless, both sites showed a differentiation of chironomid diversity. This could indicate that the mountains indeed are barrier for dispersal, creating the right conditions for allopatric speciation to occur.

Acknowledgments

I would like to thank Hans Breeuwer for his enthusiastic support in the preparation of and during practical work and his critic advice on draft versions of the report and presentation. Michiel Kraak is warmly acknowledged for sharing his views on the written work. I also want to thank Peter Kuperus and Betsy Voetdijk for their help during lab work and their everlasting patience when explaining regulations and protocols. Furthermore, Milo de Baat and Raúl Loayza-Muro are acknowleged for collecting the chironomids in Peru, this study would not have been possible without their effort.

14 References

COGNATO, A.I., 2006. Standard percent DNA sequence difference for insects does not predict species boundaries. Journal of Economic Entomology, 99(4), pp. 1037-1045.

COURTNEY, L.A. and CLEMENTS, W.H., 2002. Assessing the influence of water and substratum quality on benthic macroinvertebrate communities in a metal-polluted stream: an experimental approach.

Freshwater Biology, 47(9), pp. 1766-1778.

COYNE, J.A. and ORR, H.A., eds, 2004. Speciation. Sunderland: Sinauer Associates. CRANSTON, P., 1995. Chironomids: From Genes to Ecosystems. Melbourne: CSIRO.

CRANSTON, P.S., HARDY, N.B. and MORSE, G.E., 2012. A dated molecular phylogeny for the Chironomidae (Diptera). Systematic Entomology, 37(1), pp. 172-188.

DARWIN, C., 1859. On the Origin of Species by Means of Natural Selection. London: Murray.

DOBZHANSKY, T., 1937. Genetic Nature of Species Differences. The American Naturalist, 71(735), pp. 404-420.

DOBZHANSKY, T., 1935. A Critique of the Species Concept in Biology. Philosophy of Science, 2(3), pp. 344. EKREM, T., WILLASSEN, E. and STUR, E., 2007. A comprehensive DNA sequence library is essential for identification with DNA barcodes. Molecular Phylogenetics and Evolution, 43(2), pp. 530-542.

ETIENNE, R.S, 2005. A new sampling formula for neutral biodiversity. Ecology letters, 8(3), pp. 253-260. FOLMER, O., BLACK, M., HOEH, W., LUTZ, R. and VRIJENHOEK, R., 1994. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Molecular Marine

Biology and Biotechnology, 3(5), pp. 294-299.

HEBERT, P.D., CYWINSKA, A. and BALL, S.L., 2003. Biological identifications through DNA barcodes.

15 HOORN, C., WESSELINGH, F.P., TER STEEGE, H., BERMUDEZ, M.A., MORA, A., SEVINK, J., SANMARTÍN, I., SANCHEZ-MESEGUER, A., ANDERSON, C.L., FIGUEIREDO, J.P., JARAMILLO, C., RIFF, D., NEGRI, F.R., HOOGHIEMSTRA, H., LUNDBERG, J., STADLER, T., SÄRKINEN, T. and ANTONELLI, A., 2010. Amazonia through time: Andean uplift, climate change, landscape evolution, and biodiversity. Science, 330(6006), pp. 927-931.

HUBBELL, S.P., 2001. The unified neutral theory of biodiversity and biogeography. Princeton: Princeton

University Press.

JACOBSEN, D., 2008. Low oxygen pressure as a driving factor for the altitudinal decline in taxon richness of stream macroinvertebrates. Oecologia, 154(4), pp. 795-807.

JACOBSEN, D., ROSTGAARD, S. and VÁSCONEZ, J.J., 2003. Are macroinvertebrates in high altitude streams affected by oxygen deficiency? Freshwater Biology, 48(11), pp. 2025-2032.

JOHANNESSON, K., 2010. Are we analyzing speciation without prejudice? Annals of the New York

Academy of Sciences, 12061, pp. 143-149.

KELLY, D.J., CLARE, J.J. and BOTHWELL, M.L., 2001. Attenuation of Solar Ultraviolet Radiation by Dissolved Organic Matter Alters Benthic Colonization Patterns in Streams. Journal of the North American

Benthological Society, 20(1), pp. 96-108.

LOAYZA-MURO, R.A., MARTICORENA-RUIZ, J., PALOMINO, E.J., MERRITT, C., DE BAAT, M.L., VAN GEMERT, M., VERWEIJ, R.A., KRAAK, M.H.S. and ADMIRAAL, W., 2012. Persistence of Chironomids in Metal Polluted Andean High Altitude Streams: Does Melanin Play a Role? Environmental Science &

Technology, 47(1), pp. 601-607.

LOAYZA-MURO, R.A., ELÍAS‐LETTS, R., MARTICORENA‐RUÍZ, J.K., PALOMINO, E.J., DUIVENVOORDEN, J.F., KRAAK, M.H.S. and ADMIRAAL, W., 2010. Metal-induced shifts in benthic macroinvertebrate community composition in Andean high altitude streams. Environmental Toxicology and Chemistry, 29, pp. 2761-2768.

LOAYZA-MURO, R.A., DE BAAT, M.L., PALOMINO, E.J., KUPERUS, P., KRAAK, M.H.S., ADMIRAAL, W. and BREEUWER, J.A.J., Metal and altitude drive genetic diversity of chironomids in Andean streams. Submitted for publication

16 LOSOS, J.B. and GLOR, R.E., 2003. Phylogenetic comparative methods and the geography of speciation.

Trends in Ecology & Evolution, 18(5), pp. 220-227.

MAYR, E., 1963. Animal Species and Evolution. Cambridge: Belknap Press.

MOULTON, J.K. and WIEGMANN, B.M., 2004. Evolution and phylogenetic utility of CAD (rudimentary) among Mesozoic-aged Eremoneuran Diptera (Insecta). Molecular Phylogenetics and Evolution, 31(1), pp. 363-378.

NIXON, K.C. and WHEELER, Quentin D., 1990. An amplification of the phylogenetic species concept.

Cladistics, 6(3), pp. 211-223.

TAMURA, K., PETERSON, D., PETERSON, N., STECHER, G., NEI, M. and KUMAR, S., 2011. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Molecular Biology and Evolution, 28(10), pp. 2731-2739.

VILLAFAÑE, V.E., ANDRADE, M., LAIRANA, V., ZARATTI, F. and HELBLING, E.W., 1999. Inhibition of phytoplankton photosynthesis by solar ultraviolet radiation: studies in Lake Titicaca, Bolivia. Freshwater