BIODIVERSITY
RESEARCH
Long-term isolation and endemicity of
Neotropical aquatic insects limit the
community responses to recent
amphibian decline
Cesc M
urria1,2*, Amanda T. Rugenski3, Matt R. Whiles3 and
Alfried P. Vogler1,2
1Department of Life Sciences, Natural
History Museum, London SW7 5BD, UK,
2Department of Life Sciences, Imperial
College London, Silwood Park Campus,
Ascot SL5 7PY, UK, 3Department of Zoology
and Center for Ecology, Southern Illinois
University, Carbondale, IL, USA
*Correspondence: Cesc M
urria, D
epartement
de biologie, Universit
e de Sherbrooke,
Sherbrooke, Qu
ebec J1K 2R1, Canada.
E-mail: cmurria@ub.edu
ABSTRACT
Aim Neotropical highland streams have shown diminished ecosystem function-
ing after amphibian extirpation infected by the chytrid fungus Batrachochytrium
dendrobatidis. The loss of amphibians could affect communities of aquatic
insects co-occurring in these streams in various ways. We examined patterns of
species and genetic diversity of these communities and their evolutionary his-
tory along the chytrid expansion gradient to elucidate potential community
responses.
Location Six streams over a 320-km transect in Panama affected by chytrid
expansion from west to east for up to 14 years, and two apparently chytrid-free
streams in the east.
Methods Patterns of a- and b-diversity were investigated at three hierarchical
levels: genus, species and haplotypes. Genus identification was based on mor-
phology, and putative species were inferred by grouping the DNA barcodes
(749 cox1 sequences) with the GMYC method on all collected individuals of
Ephemeroptera, Trichoptera, Coleoptera and Plecoptera.
Results A total of 96 genera in 43 families (9 orders) of insects were encoun-
tered. Genus-level a-diversity was higher in the easternmost streams, possibly
due to a separate biogeographical history, whereas b-diversity was constant
along the chytrid expansion gradient. Community DNA barcoding resulted in
426 cox1 haplotypes and 154 putative species, most of them limited to single
sites. High b-diversity along the gradient at both species and haplotype levels
argues against community homogenization by migration in the wake of
amphibian declines. In contrast, phylo-b-diversity was low, indicating commu-
nity similarity at deep levels.
Main conclusions Aquatic insect communities in this region are influenced by
long-term limited dispersion that generated high endemicity. The pattern per-
sists mostly unperturbed after disease-driven amphibian declines; hence, if
indeed insects fill the niches vacated by tadpoles, they would originate from
local communities rather than immigration. Given the unique evolutionary his-
tory and physical isolation of local assemblages, the ecosystem deterioration
carries the risk of losing unique diversity.
Keywords
Batrachochytrium dendrobatidis, community DNA barcoding, genetic diversity,
SGDC, species diversity, stream ecology.
DOI: 10.1111/ddi.12343
ª 2015 John Wiley & Sons Ltd http://wileyonlinelibrary.com/journal/ddi 1
Diversity and Distributions, (Diversity Distrib.) (2015) 1–12
A
 J
ou
rn
al
 o
f 
Co
ns
er
va
ti
on
 B
io
ge
og
ra
ph
y
D
iv
er
si
ty
 a
nd
 D
is
tr
ib
ut
io
ns
INTRODUCTION
The accelerated rate of species extinction resulting from
human activities is closely associated with the loss of ecologi-
cal function (Hooper et al., 2005; Cardinale et al., 2006).
However, the interrelationships of species extinctions, altered
community composition and abiotic environmental changes
are complex and hamper our ability to predict the effects of
biodiversity loss (Vaughn, 2010). Current empirical
approaches have focused mainly on microcosms to examine
effects of diversity changes on system-scale functional pro-
cesses such as nutrient recycling, primary production and
organic matter dynamics (but see Atkinson et al., 2013;
Whiles et al., 2013). However, experimental studies con-
ducted over short periods on local communities may not
reflect broader evolutionary and ecological processes such as
the possible local loss of intraspecific genetic diversity and
the compensation by gene flow across metapopulations.
Moreover, the time delays between species losses and their
impacts on altered ecosystems may exceed predictions
based on small-scale experiments (Duffy, 2009). Similarly,
experimental manipulations generally do not account for
potential new colonizers with similar functional roles that
may compensate for the species loss.
Amphibian populations world-wide are experiencing mas-
sive declines, many of which are linked to Batrachochytrium
dendrobatidis (Bd), a pathogenic fungus that causes chytridi-
omycosis (Lips et al., 2006). Bd thrives in cool, humid con-
ditions, making pristine high-elevation sites in the tropics,
where endemism is high, most vulnerable. Pathogen preva-
lence is lower, and amphibian populations often remain sta-
ble at warmer lowlands (Bustamante et al., 2010; Becker &
Zamudio, 2011). Studies in Central America have docu-
mented a steady advance of Bd from west to east in Panama,
with massive die-offs of highly diverse, abundant amphibian
communities in high-altitude streams (Lips et al., 2005,
2006). For instance, the arrival of Bd coincided with the
extirpation of 40% of taxonomic diversity and 33% of the
phylogenetic diversity of frogs at one site in Panama (Craw-
ford et al., 2010). The ongoing Bd wave in Panama provides
a chronosequence of sites affected since 1996 to chytrid-free
streams (Fig. 1a). This situation represents a unique natural
(a)
(b) (c) (d)
Comm. GMYC sp. cox1 hapl.
Figure 1 (a) Distribution of sampling sites across Panama and diversity patterns after amphibian declines across the longitudinal
gradient for (b) richness of genera within communities, (c) GMYC species richness and (d) cox1 haplotypes richness. Blue circles are
total richness and red circles are richness after rarefaction based on Chuncanti, except r
ıo Blanco (black circle). Lines indicated the best
fit for a local polynomial regression for the rarefied richness.
2 Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd
C. M
urria et al.
experiment for testing fundamental questions about the
extirpation of a major trophic group and associated changes
in ecosystem function and their effects on unrelated lineages.
Larval amphibians feed on algae (grazers) and detritus and
may act as ‘ecosystem engineers’ driving nutrient recycling
because of their high abundance and consumption rates
(Flecker et al., 1999; Whiles et al., 2006). Their massive
declines have altered ecosystem function through changes in
nutrient uptake and cycling efficiency, resource availability
and overall biological activity (Whiles et al., 2006, 2013;
Col
on-Gaud et al., 2009; Rugenski et al., 2012). These
streams are also occupied by diverse communities of aquatic
insects, including some that could replace tadpoles ecologi-
cally because they have similar functional roles, such as vari-
ous Ephemeroptera, many of which are grazers, while other
species may also be affected in various direct and indirect
ways by the loss of amphibians (Col
on-Gaud et al., 2010a).
Previous local-scale studies examining the ecological effects
of amphibian extirpation on the macroinvertebrate commu-
nities comparing pre- and post-decline conditions found
subtle shifts in species composition, biomass and secondary
production (Col
on-Gaud et al., 2010a,b; Whiles et al., 2013),
whereas another study conducted over longer periods
showed a reduction in macroinvertebrate richness by ~40%
as time progressed (Rantala et al., 2015). However, all previ-
ous studies assessed changes in macroinvertebrate communi-
ties at local sites and used coarse identification to genus
because of the difficulties with species-level taxonomy. This
limits assessment of the precise species composition and its
variation among affected communities and therefore does
not permit the evaluation of local communities in the con-
text of larger-scale patterns and the evolutionary history of
the encountered lineages.
Here, we sampled previously studied sites together with
additional streams for a broader assessment of the variation
of aquatic insect communities affected by the chytrid disease
wave. The lack of an established taxonomy for insects form-
ing these local assemblages may be overcome using mito-
chondrial DNA (mtDNA) variation in the ‘barcode’ marker
cytochrome C oxidase I (cox1). Clusters of sequence varia-
tion in cox1 approximate the entities traditionally recognized
as Linnaean species in many invertebrates (Hebert et al.,
2003). This permits rapid assessment of taxonomically com-
plex species assemblages from multiple sites and various life
stages by grouping individuals using quantitative procedures
(T€anzler et al., 2012). Moreover, intraspecific variation of
cox1 is commonly used to determine genetic diversity and
assess gene flow among populations, while interspecific varia-
tion is used to infer their evolutionary histories (Avise,
2009). Thus, analysing the nucleotide variation of the cox1
gene for entire communities does not only serve to identify
species, but also integrates this information with their evolu-
tionary differentiation and current gene flow, linking pro-
cesses at different temporal scales (Marske et al., 2013; Joly
et al., 2014). Community DNA barcoding therefore has great
potential for inferring population dynamics, phylogeography,
macroecological patterns and evolutionary history of diverse
groups (Hajibabaei et al., 2007; T€anzler et al., 2012) and for
predicting large-scale biodiversity patterns using local genetic
inventories (Papadopoulou et al., 2011).
The simultaneous analysis of taxonomy, population genet-
ics and evolution of species assemblages may lead to a more
complete understanding of processes affecting the persistence
and turnover at various time-scales and the potential effects
of tadpole loss on aquatic insects. Moreover, understanding
the evolutionary mechanisms and current population dynam-
ics should be valuable for assessing the fragility of tropical
macroinvertebrate communities to perturbation and their
response to a rapidly changing environment. Species compo-
sition and their evolutionary and biogeographical history are
likely influenced by barriers to dispersal among streams and
by the vagility of species, along with the particular functional
roles of each species. The biogeography of aquatic insects in
Central America remains relatively unexplored, but high
diversity, endemicity and turnover among mountain ranges
are expected based on the stability hypothesis of habitats,
which posits that weak climatic fluctuations during the Pleis-
tocene in the tropics promoted the persistence of isolated
populations and allopatric speciation (Fjeldsa et al., 1999;
Cadena et al., 2011; Garc
ıa-L
opez et al., 2013). Similarly, in
temperate regions, the population persistence in permanent
streams has led to higher endemism and smaller geographical
species ranges for running (long-term stable habitat) than
standing (instable habitat) water species (Ribera et al., 2003;
Hof et al., 2006). Given the complexity of factors affecting
the distributions of aquatic insects and their potential
responses to amphibian declines, we assessed the diversity
and turnover using a similar analytical framework at three
hierarchical levels from the genus level to intraspecific
genetic differentiation, which provides a wider view of the
importance of the historical and contemporary factors deter-
mining current diversity patterns (Bonada et al., 2009;
Dexter et al., 2012; Marske et al., 2013).
Considering the rapidity and severity of amphibian
declines and consequent changes in ecosystem function
(Col
on-Gaud et al., 2010a,b; Whiles et al., 2013), two scenar-
ios may be considered for the effects on insect assemblages.
First, if long-term habitat stability and limited dispersal drive
the evolutionary history of aquatic insects in the region, we
expect substantial isolation among local communities even at
small spatial scales. This scenario would predict that the
post-decline communities are composed of species already
present locally, and even if the ecological roles of tadpoles
could be assumed by functionally similar aquatic insect spe-
cies, the arrival of new immigrants is unlikely. Therefore, the
local species would change in abundance after amphibian
declines, but not in composition (or some ecological func-
tions remain unfilled), and b-diversity should be high and
stay high. In that case, we predict a positive correlation
between species and genetic diversity because limited migra-
tion among discrete, undisturbed sites over extended time-
scale would have parallel effects on local diversity at these
Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd 3
Aquatic insects in highland Neotropical streams
two levels (Vellend & Geber, 2005; Vellend et al., 2014). An
alternative scenario is expected if there is historical and
recent long-distance dispersal. In this case, novel widespread
aquatic insect immigrant species may occupy available niche
space in the altered ecosystem, and we expect increasing
homogenization at both interspecific and intraspecific levels
at sites with the longest time since amphibian decline, possi-
bly at the expense of the original residents, which may
reduce b-diversity when studied across multiple sites. This
could also lead to a distortion of diversity at species and
genetic levels because rapid range expansions, possibly driven
by selection only at species level, would undermine the pro-
cesses (such as stochastic dispersal) needed for the multi-
hierarchical patterns to form (Vellend & Geber, 2005).
METHODS
Study area and sampling methods
A total of 8 streams were sampled (Fig. 1a). These included
six pristine post-decline streams located across four moun-
tain ranges, spanning from sites at Fortuna in the west to
Cerro Azul in the east and representing a chronosequence of
decline dates from 1996 to 2009. We also sampled one chy-
trid-free site in the Chucanti Nature Preserve and one low-
land pristine stream (r
ıo Frijolitos) in the Soberan
ıa National
Park (Table 1). All streams, except r
ıo Frijolitos located in
the lowlands, were sampled at similar altitudes between 700
and 900 m, or at 1200 m in the westernmost sites (Chorro,
Aleman), and were located in five unconnected rain forests
and surrounded by agricultural areas that are expanded at
the expense of rain forest. Only Aleman and Chorro were
connected by water (3.48 km apart in different subcatch-
ments), while the remaining sites were unconnected and sep-
arated by mountains. Geographical distance between
sites ranged from 2.11 km (Guabal-Blanco) to 415.85 km
(Aleman-Chucanti), with distances among most sites
between 50 and 180 km (Appendix S1). Except r
ıo Blanco
that was smaller, all streams had similar habitats and sub-
strate which ranged from large boulders and cobbles in rif-
fles, to silt and sand in pools. For example, between 2008
and 2011 for Chorro, Guabal and Chucanti, respectively, the
average width was 4.5 m, 3.7 m and 3 m; the temperature
was 19 
 2°C, 21 
 1.5°C and 21 
 2°C; and discharge was
72 
 8 l s1, 60 
 8 l s1 and 55 
 12 l s1 (Rugenski,
2013).
Sampling was conducted over a 4-week period in Febru-
ary–March of 2011 during the dry season. Stream flow had
been stable during the dry season, except r
ıo Blanco that
experienced a flood event 3 weeks prior to sampling. Macro-
invertebrates were sampled at all available habitats from rif-
fles and pools along a 150-m reach using a 250-lm mesh
kick net. Standard sampling took place for at least 150 min
until no new taxa were captured. Samples were preserved in
absolute ethanol, and macroinvertebrates were identified to
genus in the laboratory (Rold
an, 1988; Merritt et al., 2008)
prior to molecular analysis. For each genus, functional feed-
ing group assignments were based on Merritt et al. (2008).
Molecular analysis and estimates of molecular
species identities
The three most abundant orders (58 of 96 genera) encom-
passing all functional feeding groups were selected for com-
munity DNA barcoding (see functional feeding groups in
Appendix S2), including Ephemeroptera classified as grazers
and collectors (17 genera); Trichoptera, of which the major-
ity are classified as filterers and shredders, but also including
a few grazers and collectors (Glossossomatidae, Hydroptili-
dae) and predators (Polycentropodidae) (18 genera); and
Coleoptera, with the majority classified as collectors, but also
including predators (Dytiscidae, Gyrinidae) and grazers
(Psephenidae) (23 genera). Additional individuals belonging
to Plecoptera (predators) were sequenced because they
included the only genus (Anacroneuria) distributed at all
sites. The cox1 gene was sequenced for all available specimens
up to 20 individuals per genus and site (see Appendix S3 for
details of DNA sequencing).
Putative species were delimited separately for each order
using the Generalized Mixed Yule-Coalescent (GMYC) model
Table 1 Site description and diversity at genus level for all orders of aquatic insects. Year decl., year of decline; Ind., number of
individuals collected for the sampling period; Gen. Rich., number of genera per site, Chironomidae and Tipulidae were not included;
Gen.Rich.R., rarefied Gen.Rich
Site Code Area Altitude X_UTM Y_UTM Year decl. Ind. Gen. Rich. Gen. Rich.R
Chorro CO R.F. Fortuna, Chiriqu
ı 1218 364,119 963,913 1996 306 33 32
Aleman AL R.F. Fortuna, Chiriqu
ı 1206 365,157 967,239 1996 299 39 29
Blanco BL P.N. Omar Torrijos, Cocl
e 771 543,390 957,542 2004 147 19 –
Guabal GU P.N. Omar Torrijos, Cocl
e 684 545,100 958,785 2004 395 39 32
Capira CP P. N. Altos de Campana, Panam
a 808 618,314 960,341 2007 321 33 30
Frijolitos FR Pipeline Road, Panam
a 104 640,012 1,011,706 – 236 45 45
Cerro Azul CA P. N. Chagres, Panam
a 860 676,275 1,021,286 2009 291 44 41
Chucant
ı CU R. Chucanti, Panam
a 937 779,864 973,300 – 248 50 50
Total 2243 96 –
4 Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd
C. M
urria et al.
(Fujisawa & Barraclough, 2013) from the cox1 haplotype data
after all identical sequences were collapsed. To increase the
robustness and the accuracy of the GMYC model, DNA
sequences of outgroups and all available cox1 data for the
genera encountered were downloaded from GenBank (120
sequences, see Appendix S4 for GenBank accession numbers
used) and added to the final molecular data (Fig. 2). The
GMYC analysis was conducted on a maximum-likelihood
phylogenetic tree obtained with RAxML 7.0.4 (Stamatakis,
2006). The tree was made ultrametric using penalized likeli-
hood as implemented in r8s v. 1.7 (Sanderson, 2003). The
GMYC analysis was conducted using the R package Splits
(http://r-forge.r-project.org/projects/splits//) with the ‘single
threshold’ option.
(a) (b)
(c)
(d)
Figure 2 Maximum-likelihood phylogenetic trees of unique cox1 haplotypes for all four lineages: (a) Ephemeroptera, (b) Plecoptera,
(c) Coleoptera, (d) Trichoptera. Branches coloured in red correspond to variation within the inferred GMYC groups. Each haplotype
(terminal on trees) is coloured in accordance with its locality in Fig. 1. Note that haplotypes located in more than one site are coloured
with different colours. Black bars represent sequences downloaded from GenBank (see details in Appendix S4). Family names were
abbreviated to the first three letters: (a) LEH, Leptohyphidae; BAE, Baetidae; LEP, Leptophlebiidae; EUT, Euthyplociidae; CAE,
Caenidae; HEP, Heptageniidae; SIP, Siphlonurus. (b) PER, Perlidae. (c) HYO, Hydrophiloidea; HYA, Hydraenidae; GYR, Gyrinidae;
DYT, Dytiscidae; DRY, Dryopidae; PSE, Psephenidae; PTI, Ptilodactylidae; ELM, Elmidae. (d) CAL, Calamoceratidae; LEP, Leptoceridae;
HYD, Hydropsychidae; POL, Polycentropodidae; PHI, Philopotamidae.
Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd 5
Aquatic insects in highland Neotropical streams
Analysis of community diversity at genus, inter- and
intraspecific levels
To compare differences in a-diversity among sites, richness
was calculated at three levels: (1) individuals identified to
genus for the entire community, (2) individuals sequenced
and clustered into GMYC groups and (3) genetic variation
within GMYC groups (number of haplotypes and nucleotide
diversity). Measures of a-diversity were rarefied for the site
that had the lowest sample size. R
ıo Blanco exhibited signifi-
cantly lower abundance and diversity than other sites (see
below), possibly as a result of its smaller size and the recent
flood, and hence was removed from comparisons among
sites. The trend of variation in rarefied a-diversity along the
chytrid expansion was explored using a local polynomial
regression (LOESS) fitting simple models to variation of
richness site-by-site to build up a function that described the
deterministic part of the variation.
To assess turnover in communities along the chytrid disease
path, b-diversity at the three levels and phylogenetic distance
were measured for each pair of sites. Taxonomic b-diversity
was measured using the rarefied matrix of genus composition
and the Bray–Curtis dissimilarity index, which takes into
account abundance (Anderson et al., 2011). b-diversities of
GMYC groups and haplotypes were calculated using the
Sørensen index because of its analogy with the PhyloSor index
(Bryant et al., 2008), which was used for measuring the phylo-
b-diversity. The PhyloSor index measures the overlapping
fraction of phylogenetic branch lengths connecting haplotypes
between communities (Bryant et al., 2008). Values of all these
similarity indexes range from 0 (complete dissimilarity) to 1
(complete similarity) for a pair of communities. Indices of
b-diversity were also used to examine the community-similar-
ity distance decay by fitting exponential decay curves using
nonlinear regression (Soininen et al., 2007).
Nucleotide diversity p (i.e. the average number of nucleo-
tide differences per site between two sequences; Nei, 1987)
was measured for GMYC groups that had at least 5 individu-
als sequenced per site. To test for differences in mean p
among sites, a nonparametric Kruskal–Wallis test was per-
formed. Parallel evolutionary processes at species and genetic
levels (i.e. species and genetic diversity correlation (SGDC);
Vellend, 2003) were tested by correlating separately the num-
ber of genera and GMYC groups per site with the mean p
per site using Pearson’s correlation. All statistical analyses
were carried out using the R packages stats (R Development
Core Team, 2011), Vegan (Oksanen et al., 2011) and Ape
(Paradis et al., 2004).
RESULTS
At the genus level, we identified 2243 individuals belonging
to 96 genera in 43 families and 9 orders (Appendix S2). The
highest a-diversity was at Chucanti with 50 genera followed
by Frijolitos (45) and Cerro Azul (41). These three streams
were located at the easternmost side of the gradient, and
among these sites, only Cerro Azul was infected with chytrid
in 2009. The remaining four sites had similar a-diversity
ranging between 29 and 32 (Fig. 1b, Table 1). Many of the
identified genera were located across all sites; however, sev-
eral genera within the families Gomphidae (Odonata), Dyti-
scidae, Elmidae and Hydraenidae (Coleoptera) were located
only at the easternmost sites from Capira, which led to
higher richness. None of the genera exclusive to the east had
the same functional role as tadpoles (Gomphidae and
Dytiscidae: predators; Elmidae and Hydraenidae: collectors).
b-diversity at the genus level was similar along the geograph-
ical gradient and ranged from 0.53 to 0.3 and did not show
a geographical distance decay of similarity (R2 = 0.068;
P = 0.25, Fig. 3a).
(a) (b)
Figure 3 (a) Exponential distance decay of similarity fitted at the genus level (white squares; Bray–Curtis index; r2 = 0.08, P = 0.21),
species level (black circles; Sørensen index; r2 = 0.51, P < 0.001) and haplotype level (grey circles; Sørensen index; r2 = 0.85, P < 0.001).
(b) Lineal geographical distance decay of similarity (PhyloSor index) fitted for Ephemeroptera (white squares; r2 = 0.423, P = 0.001),
Trichoptera (grey circles; r2 = 0.0003, P = 0.45), Coleoptera (black triangles; r2 = 0.012, P < 0.28) and Plecoptera (black circles and
black dash line; r2 = 0.057, P = 0.09).
6 Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd
C. M
urria et al.
DNA sequences were obtained for 749 individuals repre-
senting 426 unique cox1 haplotypes (GenBank accession
numbers KR134410-KR134835) (Table 2, Appendix S5–S6).
The maximum-likelihood tree including external DNA
sequences was consistent with our genus-level identifications
in all cases (Fig. 2). GMYC analysis grouped these haplotypes
into 154 putative species (Appendix S6), Ephemeroptera had
the highest richness, followed by Trichoptera, Coleoptera and
Plecoptera (Table 2). Rarefied a-diversity at species (around
20 species per site) and genetic (around 40 haplotypes per
site) levels were similar among sites except for r
ıo Frijolitos
and Cerro Azul, which had higher diversity. Accordingly,
LOESS analysis of rarefied a-diversity at species and genetic
levels did not show a trend along the chytrid disease expan-
sion (Fig. 1C–D).
The mtDNA genealogies were characterized by high intra-
specific and interspecific genetic structure at local sites
(Appendix S5-S7). Almost all haplotypes were located only at
one site (mono-coloured bar in Fig. 2), with the exception
of 19 haplotypes that were located in two close sites, two
haplotypes located in three sites, and only one haplotype
(Anchytarsus, Coleoptera) found across the five westernmost
sites. GMYC groups mirrored the structure of haplotypes,
and generally, each GMYC group was located exclusively at
one site. Rarely, GMYC groups were distributed in two
nearby sites or, in the coleopteran genera Psephenus, Anchy-
tarsus and Heterelmis, at four or five sites.
b-diversity at the species level ranged from 0.39 to 0 and,
at the haplotype level, ranged from 0.22 to 0 (Fig. 3a). For
both levels, the highest values of b-diversity were between
the closely adjacent r
ıo Chorro and Alem
an. At the species
level, communities separated by more than 300 km showed
complete turnover, and this distance was reduced to
~100 km at the intraspecific level. Nonlinear regressions of
community-similarity distance decay fit an exponential decay
for species (r2 = 0.51, P < 0.001) and genetic (r2 = 0.85,
P = 0.001) levels. In contrast, phylo-b-diversity for Trichop-
tera, Coleoptera and Plecoptera ranged from 0.6 to 0.2 and
did not show any geographical distance decay of similarity
(P > 0.05), while Ephemeroptera showed a moderate linear
decay (r2 = 0.423, P = 0.001) (Fig. 3b). The largely constant
values of phylo-b-diversity indicated that turnover of vari-
ability was at the terminal levels (species or haplotypes), but
the same major lineages were represented at all sites and
diversified locally.
Genetic diversity was measured for 60 GMYC groups that
ranged from 5 to 21 individuals sequenced per site
(mean = 8.3). Nucleotide diversity p ranged from 0 to 0.013,
and these values were independent of location of the site
along the gradient (Fig. 4a) and the number of individuals
sequenced per site or order (Appendix S7). For instance, the
highest values of p were found for Farrodes (Ephemeroptera)
at Guabal (p = 0.0135) and Cerro Azul (p = 0.0116) that
had 10 and 13 individuals sequenced, respectively, whereas
the 9 individuals sequenced for Farrodes from Aleman
showed p = 0.00. At Capira, the genus Farrodes had two
GMYC groups that showed dissimilar values of p = 0.0105
and p = 0.0021 for 6 and 7 individuals sequenced, respec-
tively. The only possible generalization concerns the genus
Anacroneuria (Plecoptera) that showed consistently low p
(0.00–0.0033). The highest values of p were found in
Ephemeroptera, whereas Coleoptera and Trichoptera showed
intermediate values. Cerro Azul showed the highest mean p
per site; however, there were no significant differences
among sites (v2 = 10.64, P = 0.1). Species and genetic diver-
sities at local sites were positively correlated (r2 = 0.72,
P = 0.007), whereas neither was significantly correlated at
the genus level (r2 = 0.27, P = 0.18) (Fig. 4b).
DISCUSSION
Our study reveals high turnover among streams simulta-
neously at the species and genetic levels, which together sug-
gests that dispersal between the sampled areas has been limited
over extended time-scales. We also find that a-diversity at spe-
cies and genetic levels did not significantly change along the
chytrid disease expansion; hence, the time since amphibian
extirpation (from up to 14 years to chytrid-free sites) was not
Table 2 Number of species (GMYC) and cox1 haplotypes (hap) for the four selected insect orders. Ind., number of individuals
sequenced; NucDiv, mean nucleotide diversity p per site for the GMYC groups with at least 5 individuals sequenced
Site
Ephemeroptera Trichoptera Coleoptera Plecoptera Total
Ind. GMYC Hap Ind. GMYC Hap Ind. GMYC Hap Ind. GMYC Hap Ind. GMYC Hap NucDiv*
CO 48 14 25 19 2 12 14 3 10 5 2 3 86 21 50 0.00254
AL 46 14 26 22 7 16 2 1 2 6 3 6 76 23 50 0.00403
BL 0 0 0 23 6 12 1 1 1 12 1 2 36 8 15 0.00077
GU 62 14 35 31 6 14 45 9 32 11 2 5 149 31 86 0.00392
CP 22 11 19 53 5 13 20 6 13 8 2 2 103 23 47 0.00318
FR 64 15 40 34 9 22 16 6 12 9 3 5 123 33 79 0.00391
CA 68 17 51 14 7 12 28 6 20 7 2 3 117 35 86 0.00755
CU 15 4 14 20 8 14 14 7 7 10 4 7 59 22 42 0.00266
Total 325 69 202 216 42 104 140 29 90 68 14 30 749 154 426
*See Appendix S5–S6 for details.
Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd 7
Aquatic insects in highland Neotropical streams
correlated with the local richness of aquatic insects. The high
b-diversity and rapid exponential distance decay of similarity
detected at both organizational levels are evidence against any
large-scale processes sweeping through these insect communi-
ties and the homogenization through widespread species that
might mirror the chytrid wave. This suggests that the effect of
amphibian decline on the pattern of species and intraspecific
diversity at large spatial scale is minimal.
Existing ecological studies have established changes in eco-
system function (Whiles et al., 2006, 2013; Col
on-Gaud
et al., 2009; Rugenski et al., 2012) and declines of insect
diversity several years after the loss of amphibians (Rantala
et al., 2015). These assessments were generally limited to par-
ticular local sites, hence lacking the regional perspective, and
they were based on the unrefined identification to genus level
only. Our phylogenetic analysis from DNA data extended
across the entire region provides a novel, evolutionarily
explicit framework to discriminate between different possible
scenarios for the response of aquatic insect communities.
The significant changes in ecosystem function following the
tadpole loss at local sites (e.g. increase in community respira-
tion and N flux to the epilithon, or decrease in N uptake
length and N flux to fine particulate organic matter; Whiles
et al., 2013) are apparently not compensated for by other
components of the ecosystem such as aquatic macroinverte-
brates (Rantala et al., 2015). Our results are consistent with
this scenario, as the post-decline macroinvertebrate commu-
nities show very localized distributions that suggest they are
not composed of recent immigrants that would occupy the
vacant niche space after tadpole extirpations. With our data,
we cannot exclude that species in post-decline sites might
have been drawn from a local pool of in situ evolving
lineages, which may have been restricted to pockets of the
local sites, or existed at very low abundance when tadpoles
were present and now increased in abundance (which we did
not measure). While insects do not profit from the disap-
pearance of tadpoles, it remains to be seen whether and how
these profound ecological changes linked to shifts in basal
resources associated with the loss of tadpoles will have direct
(e.g. competition for food resources) or indirect (e.g. biotur-
bation, altering algal community structure) effects on macro-
invertebrate communities (Rantala et al., 2015), which may
ultimately result in the disappearance of isolated aquatic
macroinvertebrate communities.
The only feature unique to the chytrid-free sites was the
composition of communities in the genus-level study, which
detected several taxa limited to the easternmost side of the
geographical gradient and a-diversity increased abruptly
from Capira eastwards. The absence of these specific taxa at
the westernmost sites was confirmed by species composition
datasets from 2008 to 2011 from Chorro and Guabal
(Col
on-Gaud et al., 2010a,b; Rugenski, 2013). Given that
these taxa use resources different from tadpoles and are in
different functional feeding groups, and furthermore, Cerro
Azul experienced declines two years before our study but still
harboured these groups, it is not likely that the lower rich-
ness at the western side was associated with amphibian extir-
pation. Instead, the pattern of a-diversity may be explained
by the historical disjunction of the North and South Ameri-
can fauna around the Isthmus of Panama.
The limited sampling effort may also affect the estimates
of b-diversity, given our snapshot sampling over a 1-day
(a) (b)
Figure 4 (a) Distribution of nucleotide diversity (p) across the longitudinal gradient for all GMYC groups that had at least 5
individuals sequenced: Ephemeroptera (white squares), Trichoptera (grey circles), Coleoptera (black triangles) and Plecoptera (black
circles). (b) Species and genetic diversity correlation (SGDC). Species diversity refers the number of genera per site (grey dots, r2 = 0.27,
P = 0.18) and GMYC groups per site (black dots, r2 = 0.72, P = 0.007). For the calculations of the local genetic diversity, the mean
nucleotide diversity was established for each GMYC group at a site for those represented by at least 5 sequenced individuals (Appendix
S6).
8 Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd
C. M
urria et al.
period that may be incomplete and therefore artificially
increase the inferred species turnover. However, it is unlikely
that a pattern of high endemicity would be created simply
from undersampling, as it would favour the detection of
widespread and common species, which we do not find. At
the genus level, numerous taxa were recovered consistently
among sites and their composition was comparable with
other studies in this region (Col
on-Gaud et al., 2009, 2010a,
b), and in the case of the locality at Chucanti, our repeated
sampling over 3 years revealed a very similar community
composition (Rugenski, 2013). This suggests that our sample
is a good representation of the local biota, and therefore, the
turnover within genera at the species and genetic level is
likely to be an adequate reflection of the true change among
sites. The detected turnover also may be a consequence of
other factors such as seasonality and patchiness of macroin-
vertebrate communities associated with flow disturbances
(e.g. the low diversity found at r
ıo Blanco was presumably
affected by a recent flooding event) (Lake, 2000; R
ıos-Touma
et al., 2011). However, even with limited sampling, the rare-
fied species and genetic diversities were remarkably similar at
the 8 study sites, with the greatest proportion of species
belonging to Ephemeroptera, approximately half as many
species of Coleoptera and Trichoptera, and only a few species
to Plecoptera. The largely uniform number of mayfly species,
all of which are considered to be part of the grazer and col-
lector guild, is a strong indication of their ubiquity and eco-
logical significance, while at the species and genetic levels,
there is strong turnover among sites.
The high interspecific and intraspecific turnover among
sites supports the notion of long-term habitat stability in
tropical highlands that resulted in limited dispersion and an
increase in genetic differentiation among mountain sites
(Cadena et al., 2011; Garc
ıa-L
opez et al., 2013). Similarly,
high community turnover at the levels of haplotypes, species
and higher clades was found in terrestrial scarab beetles
across several Neotropical highland forests, which was attrib-
uted to the long-term ecological stability of local populations
that may track suitable habitats to accommodate climatic
changes without long-distance dispersal (Garc
ıa-L
opez et al.,
2013). The aquatic insects examined here showed even
higher genetic structure among communities at a similar
geographical scale, possibly because dispersal among streams,
especially among headwaters, is even more limited than
among the largely contiguous terrestrial rain forest habitats
inhabited by scarabs (Bohonak & Jenkins, 2003; Finn et al.,
2011). Long-term habitat stability may be surprising given
the undoubted effects of Pleistocene climatic fluctuations on
species movements, especially in the Northern Hemisphere
(Hewitt, 2000), but in the Tropics, the resulting shifts in
habitat distributions are comparatively small and in case of
highland habitats can be tracked by populations up and
down mountain slopes with minimal movements (Rull,
2005; Graham et al., 2006; M
urria & Hughes, 2008). Popula-
tion persistence in unstable habitats requires high dispersal
rate to ensure survival, and accordingly, species tend to
exhibit lower genetic endemism and larger ranges than those
in long-term stable habitats (Ribera et al., 2003; Papadopou-
lou et al., 2009). Habitat stability would amplify patterns of
high b-diversity over geographical distance and the size of
regional species pools in tropical rain forests (Graham et al.,
2006; Carnaval et al., 2009) and also preserve biogeographi-
cal signatures at higher hierarchical levels.
Nucleotide diversity p was highly variable among and
within lineages (except Plecoptera that had the lowest p) and
among sites (Fig. S7). Distribution of p in mtDNA among
species has been correlated with differences in relevant traits
such as long-lived or low fecundity (Romiguier et al., 2014),
recurrent adaptive evolution in mtDNA genes (Bazin et al.,
2006) or geographical range and dispersal abilities of species
(Fujisawa et al., 2015). In our study, the disparity of p
within and among lineages is not predicted from differences
in geographical ranges, dispersal abilities and the correlated
population size (Fujisawa et al., 2015) because almost all spe-
cies have narrow, isolated geographical distributions. The
high disparity of p shown by closely related species that
share functional traits and inhabit the same isolated stream
(e.g. the two species Farrodes at Capira) may refute the
hypothesis of differences in adaptive traits and the explana-
tions based on habitat stability because these species presum-
ably have similar population history associated with local
abiotic factors such as flood events or biogeography. It is
unlikely that these observations are explained by sampling
artefacts, as the number of individuals sequenced per GMYC
groups and p were not correlated (Fig. S7). The high dispar-
ity of p remains unexplained but may be the results of con-
tingent factors affecting the long-term effective population
size of each species, as species abundance or frequent bottle-
necks affect comparatively small populations in a fluctuating
environment (Fujisawa et al., 2015), or may reflect the time
since the last occurrence of a selective sweep for each species
(Bazin et al., 2006). Regardless of the determinant of the
detected high variability of p, there was no significant shift
among the pre- and post-decline sites, indicating that dispar-
ity in p is not an initial indication of declining populations
following amphibian extirpation.
We did, however, find a correlation between number of
GMYC groups and intraspecific diversity of the local species
present at these sites (confirming the SGDC of Vellend,
2003). This correlation is usually explained by the action of
uniform processes that generate patterns at both hierarchical
levels, such as dispersal among sites that affect the total spe-
cies counts as the product of immigration–extinction dynam-
ics and affect genetic diversity as the result of genetic drift.
The detected high turnover among sites at both hierarchical
levels suggests that indeed extended isolation of local assem-
blages (i.e. a neutral process) is the main mechanism causing
the observed macroecological pattern, rather than selection
for particular species traits (e.g. those relevant to niches
vacated by amphibians). This is further support for an evolu-
tionary history of isolation and restricted migration that lim-
its species ranges and the potential responses of aquatic
Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd 9
Aquatic insects in highland Neotropical streams
insects to amphibian declines. Just as the loss of endemic
amphibian assemblages resulted in numerous species extinc-
tions in the wake of chytrid infection (Crawford et al.,
2010), each Neotropical highland stream harbours a unique
diversity of aquatic insects that is highly erodible and irre-
placeable because of long-term isolation. Local extinctions
caused by anthropogenic disturbances, for example from
agriculture, logging, mining and wind farms, are not likely to
be compensated for by recolonization because of limited dis-
persion among mountain ranges.
ACKNOWLEDGEMENTS
This research was supported by National Science Foundation
Grant DEB #0717741 and Leverhulme Trust Grant F/00696/
P. We thank The Smithsonian Tropical Research Institute
and Autoridad Nacional del Ambiente for logistical support
and sampling permits, and Guido Berguido for his assistance
with sampling in the Chuncanti nature preserve. All research
complies with the current laws of the Republic of Panama,
as stated in the scientific permits SE/A-25-10 and SE ⁄ A-40-
10. CM was supported by a Beatriu de Pin
os postdoctoral
fellowship (BP-DGR-2011) from Agencia de Gesti
o d’Ajuts
Universitaris i de Recerca, Catalunya.
REFERENCES
Anderson, M.J., Crist, T.O., Chase, J.M., Vellend, M., Inouye,
B.D., Freestone, A.L., Sanders, N.J., Cornell, H.V., Comita,
L.S., Davies, K.F., Harrison, S.P., Kraft, N.J.B., Stegen, J.C.
& Swenson, N.G. (2011) Navigating the multiple meanings
of b diversity: a roadmap for the practising ecologist. Ecol-
ogy Letters, 14, 19–28.
Atkinson, C., Vaughn, C.C., Forshay, K.J. & Cooper, J.T.
(2013) Aggregated filter-feeding consumers alter nutrient
limitation: consequences for ecosystem and community
dynamics. Ecology, 94, 1359–1369.
Avise, J.C. (2009) Phylogeography: retrospect and prospect.
Journal of Biogeography, 36, 3–15.
Bazin, E., Gl
emin, S. & Galtier, N. (2006) Population size
does not influence mitochondrial genetic diversity in ani-
mals. Science, 312, 570–572.
Becker, G.C. & Zamudio, K.R. (2011) Tropical amphibian
populations experience higher disease risk in natural habi-
tats. Proceedings of the National Academy of Sciences USA,
108, 9893–9898.
Bohonak, A.J. & Jenkins, D.G. (2003) Ecological and evolu-
tionary significance of dispersal by freshwater invertebrates.
Ecology Letters, 6, 783–796.
Bonada, N., M
urria, C., Zamora-Mu~noz, C., El Alami, M.,
Poquet, J.M., Punt
ı, T., Moreno, J.L., Bennas, N.,
Alba-Tercedor, J., Ribera, C. & Prat, N. (2009) Using com-
munity and population approaches to understand how
contemporary and historical factors have shaped species
distribution in river ecosystems. Global Ecology and Bioge-
ography, 18, 202–213.
Bryant, M.J., Lamanna, C., Morlon, H., Kerkhoff, A.J., En-
quist, B.J. & Green, J.L. (2008) Microbes on moutainsides:
contrasting elevational patterns of bacteria and plant diver-
sity. Proceedings of the National Academy of Sciences USA,
105, 11505–11511.
Bustamante, H.M., Livo, L.J. & Carey, C. (2010) Effects of
temperature and hydric environment on survival of the
Panamanian Golden Frog infected with a pathogenic chy-
trid fungus. Integrative Zoology, 5, 143–153.
Cadena, C.D., Kozak, K.H., G
omez, J.P., Parra, J.L., Mccain,
C.M., Bowie, R.C., Carnaval, A.C., Moritz, C., Rahbek, C.,
Roberts, T.E., Sanders, N.J., Schneider, C.J., Vanderwal, J.,
Zamudio, K.R. & Graham, C.H. (2011) Latitude, elevation
climatic zonation and speciation in New World vertebrates.
Proceedings of the Royal Society B, 279, 194–201.
Cardinale, B.J., Srivastava, D.S., Duffy, J.E., Wright, J.P.,
Downing, A.L., Sankaran, M. & Jouseau, C. (2006) Effects
of biodiversity on the functioning of trophic groups and
ecosystems. Nature, 443, 989–992.
Carnaval, A.C., Hickerson, M.J., Haddad, C.F.B., Rodrigues,
M.T. & Moritz, C. (2009) Stability Predicts Genetic diversity
in the Brazilian Atlantic Forest hotspot. Science, 323, 785–789.
Col
on-Gaud, C., Whiles, M.R., Kilham, S.S., Lips, K.R., Prin-
gle, C.M., Connelly, S. & Peterson, S.D. (2009) Assessing
ecological responses to catastrophic amphibians declines:
patterns of macroinvertebrates production and food web
structure in upland Panamanian streams. Limnology and
Oceanography, 54, 331–343.
Col
on-Gaud, C., Whiles, M.R., Brenes, R., Kilham, S.S., Lips,
K.R., Pringle, C.M., Connelly, S. & Peterson, S.D. (2010a)
Potential functional redundancy and resource facilitation
between tadpoles and insect grazers. Freshwater Biology, 55,
2077–2088.
Col
on-Gaud, C., Whiles, M.R., Lips, K.R., Pringle, C.M., Kil-
ham, S.S., Connelly, S., Brenes, R. & Peterson, S.D.
(2010b) Stream invertebrate responses to a catastrophic
decline in consumer diversity. Journal of North American
Benthological Society, 29, 1185–1198.
Crawford, A.J., Lips, K.R. & Bermingham, E. (2010) Epi-
demic disease decimates amphibian abundance, species
diversity, and evolutionary history in the highlands of cen-
tral Panama. Proceedings of the National Academy of Sci-
ences USA, 107, 13777–13782.
Dexter, K.G., Terborgh, J.W. & Cunningham, C.W. (2012)
Historical effects on beta diversity and community assem-
bly in Amazonian trees. Proceedings of the National Acad-
emy of Sciences USA, 109, 7787–7792.
Duffy, J.E. (2009) Why biodiversity is important to the func-
tioning of real-world ecosystems. Frontiers in Ecology and
the Environment, 7, 437–444.
Finn, D.S., Bonada, N., M
urria, C. & Hughes, J.M. (2011)
Small but mighty: headwaters are vital to stream network
biodiversity at two levels of organization. Journal of North
American Benthological Society, 30, 963–980.
Fjeldsa, J., Lambin, E. & Mertens, B. (1999) Correlation
between endemism and local ecoclimatic stability
10 Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd
C. M
urria et al.
documented by comparing Andean bird distribution and
remotely sensed land surface data. Ecography, 22, 63–78.
Flecker, A.S., Feifarek, B.P. & Taylor, B.W. (1999) Ecosystem
engineering by a tropical tadpole: density-dependent effects
on habitat structure and larval growth rates. Copeia, 2,
495–500.
Fujisawa, T. & Barraclough, T.G. (2013) Delimiting species
using single-locus data and the Generalized Mixed Yule
Coalescent approach: a revised method and evaluation on
simulated data sets. Systematic Biology, 62, 707–724.
Fujisawa, T., Vogler, A.P. & Barraclough, T.G. (2015) Ecol-
ogy has contrasting effects on genetic variation within spe-
cies versus rates of molecular evolution across species in
water beetles. Proceedings of the Royal Society B, 282,
20142476.
Garc
ıa-L
opez, A., Mic
o, E., M
urria, C., Galante, E. & Vogler,
A.P. (2013) Beta-diversity at multiple hierarchical levels:
explaining the high diversity of scarab beetles in tropical
montane forests. Journal of Biogeography, 40, 2134–2145.
Graham, C.H., Moritz, C. & Williams, S.E. (2006) Habitat
history improves prediction of biodiversity in rainforest
fauna. Proceedings of the National Academy of Sciences USA,
103, 632–636.
Hajibabaei, M., Singer, G.A., Hebert, P.D.N. & Hickey, D.A.
(2007) DNA barcoding: how it complements taxonomy,
molecular phylogenetics and population genetics. Trends in
Genetics, 24, 167–172.
Hebert, P.D.N., Cywinska, A., Ball, S.L. & deWaard, J.R.
(2003) Biological identifications through DNA barcodes.
Proceedings of the Royal Society B, 270, 313–321.
Hewitt, G.M. (2000) The genetic legancy of the Quaternary
ice ages. Nature, 68, 87–112.
Hof, C., Br€andle, M. & Brandl, R. (2006) Lentic odonates
have larger and more northern range than lotic species.
Journal of Biogeography, 33, 63–70.
Hooper, D.U., Chapin, F.S., Ewel, J.J., Hector, A., Inchausti,
P., Lavorel, S., Lawton, J.H., Lodge, D.M., Loreau, M.,
Naeem, S., Schmid, B., Setala, H., Symstad, A.J., Vander-
meer, J. & Wardle, D.A. (2005) Effects of biodiversity on
ecosystem functioning: a consensus of current knowledge.
Ecological monographs, 75, 3–35.
Joly, S., Davies, T.J., Archambault, A., Bruneau, A., Alison,
D., Kembel, S.W., Peres-Neto, P.R., Vamosi, J.C. &
Wheeler, T.A. (2014) Ecology in the age of DNA barcod-
ing: the resource, the promise and the challenges ahead.
Molecular Ecology Resources, 14, 221–232.
Lake, P.S. (2000) Disturbance, patchiness, and diversity in
streams. Journal of North American Benthological Society,
19, 573–592.
Lips, K.R., Brem, F., Brenes, R., Reeve, J.D., Alford, R.A.,
Voyles, J., Carey, C., Livo, L., Pessier, A.P. & Collins, J.P.
(2006) Emerging infectious disease and the loss of biodi-
versity in a Neotropical amphibian community. Proceedings
of the National Academy of Sciences USA, 103, 3165–3170.
Lips, K.R., Burrowes, P.A., Mendelson, J.R. & Parra-Olea, G.
(2005) Amphibian declines in Latin America: widespread
population declines, extinctions, and impacts. Biotropica,
37, 163–165.
Marske, K.A., Rahbek, C. & Nogu
es-Bravo, D. (2013) Phy-
logeography: spanning the ecology-evolution continuum.
Ecography, 36, 1169–1181.
Merritt, R.W., Cummins, K.W. & Berg, M.B. (2008) An
introduction to aquatic insects of North America, 4th edn.
Kendall/Hunt Publishing Company, Dubuque, IA.
M
urria, C. & Hughes, J.M. (2008) Cyclic habitat displace-
ments during Pleistocene glaciations have induced inde-
pendent evolution of Tasimia palpata populations
(Trichoptera:Tasimiidae) in isolated subtropical rain forest
patches. Journal of Biogeography, 35, 1727–1737.
Nei, M. (1987) Molecular Evolutionary Genetics. Columbia
University Press, New York.
Oksanen, J., Blanchet, F.G., Kindt, R., Legendre, P., Minchin,
P.R., O’hara, R.B., Simpson, G.L., Solymos, P., Stevens,
M.H.H. & Wagner, H. (2011) Vegan: Community Ecology
Package. R package version 2.0-2.
Papadopoulou, A., Anastasiou, I., Keskins, B. & Vogler, A.P.
(2009) Comparative phylogeography of tenebrionid beetles
in the Aegean archipelago: the effects of dispersal ability
and habitat preference. Molecular Ecology, 18, 2503–2517.
Papadopoulou, A., Anastasiou, I., Spagopoulou, F., Stalime-
rou, M., Terzopoulou, S., Legakis, A. & Vogler, A.P. (2011)
Testing the Species-Genetic Diversity Correlations in the
Aegean Archipelago: toward a haplotype-based macroecol-
ogy? The American Naturalist, 178, 241–255.
Paradis, E., Claude, J. & Strimmer, K. (2004) APE: analyses
of phylogenetics and evolution in R language. Bioinformat-
ics, 20, 289–290.
R Development Core Team (2011) R: a language and envi-
ronment for statistical computing. R Foundation for Statisti-
cal Computing, Vienna, Austria. Available at: http://
www.R-project.org/.
Rantala, H.M., Nelson, A.M., Fulgoni, J.N., Whiles, M.R.,
Hall, R.O., Dodds, W.K., Verburg, P., Huryn, A.D., Prin-
gle, C.M., Kilham, S.S., Lips, K.R., Colon-Gaud, C., Rugen-
ski, A.T., Peterson, S.D., Fritz, K., McLeran, K.E. &
Connelly, S. (2015) Long-term changes in structure and
function of a tropical headwater stream following a
disease-driven amphibian decline. Freshwater Biology, 60,
575–589.
Ribera, I., Foster, G.N. & Vogler, A.P. (2003) Does habitat
use explain large scale species richness patterns of aquatic
beetles in Europe? Ecography, 26, 145–152.
R
ıos-Touma, B., Encalada, A.C. & Prat, N. (2011) Macroin-
vertebrate assemblages of an andean high-altitude tropical
stream: the importance of season and flow. International
Review of Hydrobiology, 96, 667–685.
Rold
an, G. (1988) Gu
ıa para el estudio de los macroinverte-
brados acu
aticos del Departamento de Antioquia, edn. Col-
ciencias-Fondo FEN, Bogot
a.
Romiguier, J., Gayral, P., Ballenghien, M. et al. (2014) Com-
parative population genomics in animals uncovers the
determinants of genetic diversity. Nature, 515, 261–263.
Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd 11
Aquatic insects in highland Neotropical streams
Rugenski, A.T., M
urria, C. & Whiles, M.R. (2012) Tadpoles
enhance microbial activity and leaf decomposition in a
neotropical headwater stream. Freshwater Biology, 57,
1904–1913.
Rugenski, A.T. (2013). Influence of disease-driven amphibian
declines on ecosystem structure and function in Panamanian
headwater streams. PhD thesis. Southern Illinois University,
USA.
Rull, V. (2005) Biotic diversification in the Guayana
Highlands: a proposal. Journal of Biogeography, 32, 921–
927.
Sanderson, M.J. (2003) r8s: Inferring absolute rates of molec-
ular evolution and divergence times in the absence of a
molecular clock. Bioinformatics, 19, 301–302.
Soininen, J., Mcdonald, R. & Hillebrand, H. (2007) The dis-
tance decay of similarity in ecological communities. Ecogra-
phy, 30, 3–12.
Stamatakis, A. (2006) RAxML-VI-HPC: maximum likeli-
hood-based phylogenetic analyses with thousands of taxa
and mixed models. Bioinformatics, 22, 2688–2690.
T€anzler, R., Sagata, K., Surbakti, S., Balke, M. & Riedel, A.
(2012) DNA barcoding for community ecology - how to
tackle a hyperdiverse, mostly undescribed Melanesian
fauna. PLoS ONE, 7, e28832.
Vaughn, C.C. (2010) Biodiversity losses and ecosystem func-
tion in freshwaters: emerging conclusions and research
directions. BioScience, 60, 25–35.
Vellend, M. (2003) Island biogeography of genes and species.
The American Naturalist, 162, 358–365.
Vellend, M. & Geber, M.A. (2005) Connections between spe-
cies diversity and genetic diversity. Ecology letters, 8, 767–
781.
Vellend, M., Lajoie, G., Bourret, A., M
urria, C., Kembel,
S.W. & Garant, D. (2014) Drawing ecological inferences
from coincident patterns of population- and community-
level biodiversity. Molecular Ecology, 23, 2890–2901.
Whiles, M.R., Hall, R.O. Jr, Dodds, W.K., Verburg, P.,
Huryn, A.D., Pringle, C.M., Lips, K.R., Kilham, S.S.,
Col
on-Gaud, C., Rugenski, A.T., Peterson, S.D. & Connel-
ly, S. (2013) Disease-driven amphibian declines alter
ecosystem processes in a tropical stream. Ecosystems, 16,
146–157.
Whiles, M.R., Lips, K.R., Pringle, C.M., Kilham, S.S., Bixby,
R.J., Brenes, R., Connelly, S., Colon-Gaud, J.C., Hunte-
Brown, M., Huryn, A.D., Montgomery, C. & Peterson, S.
(2006) The effects of amphibian population declines on the
structure and function of Neotropical stream ecosystems.
Frontiers in Ecology and Evolution, 4, 27–34.
SUPPORTING INFORMATION
Additional Supporting Information may be found in the
online version of this article:
Appendix S1 Geographical distances among sites.
Appendix S2 Taxonomical composition per site.
Appendix S3 Methods of molecular analysis.
Appendix S4 Cox1 sequences downloaded from GenBank.
Appendix S5 Individual identification (BMNH number) of
all cox1 sequenced.
Appendix S6 Taxonomical information of each GMYC
group and haplotype, and their location.
Appendix S7 Table S7: Number of individuals sequenced,
Nucleotide Diversity (p) and genetic diversity (GD) per site;
Figure S7: Relationship between number of individuals
sequenced and Nucleotide Diversity (p).
BIOSKETCHES
Cesc M
urria is a postdoctoral researcher at Universit
e de
Sherbrooke (Quebec) where he examines macroecological
patterns of diversity at the population and community level
and their correlation.
Amanda Rugenski is a postdoctoral researcher at Arizona
State University where she studies how species traits and
consumer-driven nutrient recycling influence structure and
function in aquatic ecosystems.
Matt Whiles’ laboratory at Southern Illinois University
focuses on quantifying roles of consumers in freshwater eco-
system function.
Alfried Vogler’s laboratory at the Natural History Museum,
London, uses molecular phylogenetic and evolutionary analy-
ses to study species diversity and genomics of insects.
Author contributions: All authors conceived the study, C.M.
and A.T.R. did the field work, taxonomic identifications and
molecular analysis, C.M. performed all statistical and phylo-
genetic analysis, and figures. All authors discussed the results,
C.M. and A.P.V. wrote a draft of the paper, and all authors
contributed substantially to the revision.
Editor: Jacqueline Beggs
12 Diversity and Distributions, 1–12, ª 2015 John Wiley & Sons Ltd
C. M
urria et al.