Headwaters (stream orders 1–2) traditionally have been considered depauperate compared to mid-order streams (orders 3–4)—a conclusion that arises from a perception of streams as linear systems and emphasizes change in average α (local) diversity along streams. We hypothesized an opposite pattern for β (among-site) diversity and suggest that headwaters might account for a large degree of basin-scale biodiversity if considered within the more realistic framework of streams as branching networks. We assembled pre-existing biodiversity data from across the globe to test this hypothesis broadly at the population-genetic (mitochondrial haplotype diversity within species) and community (species/taxonomic diversity) levels, with a focus on macroinvertebrates. We standardized 18 (9 headwater and 9 mid-order) population-genetic and 16 (10 headwater and 6 mid-order) community-level ecoregional data sets from 5 global ecozones for robust comparisons of β-diversity estimates between the 2 stream-size categories. At the population-genetic level, we applied measures of among-site variation commonly used at both population-genetic (FST and ΦST) and community (Sørensen's dissimilarity with both presence/absence and abundance data) levels and developed a novel strategy to compare expected rates of loss of γ (regional) diversity as individual sites are eliminated sequentially from regions. At the community level, we limited analyses to Sørensen's presence/absence measures. We found that Sørensen's dissimilarity was significantly greater among headwaters than among mid-order streams at both population-genetic and community levels. We also showed that individual headwater reaches accounted for greater proportions of genetic γ diversity than did mid-order reaches. However, neither FST nor ΦST was significantly different between stream-size categories. These measures, which have been used traditionally for comparisons of population-genetic variation, measure proportions of total variation rather than solely among-site variation (i.e., they also are influenced by within-site variation). In contrast, Sørensen's dissimilarity measures only among-site variation and, therefore, is presumably more useful for reflecting general β diversity. Overall results suggest that, on average, headwaters probably contribute disproportionately to biodiversity at the network scale. This finding demands a shift in thinking about the biodiversity contributions of small headwaters and has strong conservation implications for imperiled headwater streams around the world.
Three decades ago, the influential River Continuum Concept (RCC; Vannote et al. 1980) modeled broad-scale spatial structure and function in streams in terms of gradual and predictable change from headwaters to mouth. One pattern predicted by the RCC was a unimodal distribution of biological diversity, with a peak in mid-order streams and lower values in headwaters and large rivers. Perhaps because small streams have been studied more extensively, the upper ½ of this predicted biodiversity pattern—an increase from headwaters to mid-orders—has been evaluated and supported in many stream types, particularly high-gradient mountain streams (Minshall et al. 1985, Ward 1986, Finn and Poff 2005, Sheldon and Warren 2009). These studies measured species or taxonomic diversity at the level of biological organization to which we refer henceforth as the community level. Of course, biodiversity includes multiple levels of biological organization (NRC 1999), and the pattern of increasing diversity with downstream distance from 1st to ∼4th- or 5th-order streams also is apparent at the level of intraspecific genetic diversity (henceforth, the population-genetic level). In other words, for many common species distributed from headwaters to mid-order streams, genetic diversity within species increases from upstream to downstream (Shaw et al. 1994, Crispo et al. 2006, Hänfling and Weetman 2006, Watanabe et al. 2008).
We make the above observations under the linear model of streams embodied by the RCC. Therefore, this linear view considers headwaters to be relatively depauperate (at either biological level), contributing little to whole-stream biodiversity. However, streams are not linear. Streams are branching networks, and this alternative to linearity has been conceptualized in various ways during the past decade (Fagan 2002, Benda et al. 2004, Grant et al. 2007, Morrissey and de Kerckhove 2009). Fisher (1997) warned that the long-influential linear view of streams is “at best incomplete and at worst incorrect”. With respect to stream biodiversity distribution, we expect the best—that the linear view is incomplete (although potentially misleading). Biodiversity across any landscape (or riverscape) can be decomposed into 3 fundamental components: local (α) diversity, regional (γ) diversity, and some description of diversity turnover or variation from locality to locality within a region (β diversity) (Whittaker 1960, Anderson et al. 2011). Patterns of α diversity in streams potentially can be accommodated within the linear conceptualization (as above), but this view ignores β and γ components of diversity.
If we consider a branching-network view of streams, all 3 components should be essential to describe diversity distributions. A brief examination of a network reveals that small streams are much more numerous than large streams. Headwaters (1st- and 2nd-order streams) are thought to comprise >75% of total stream length in most basins (Leopold et al. 1964, Benda et al. 2004). Their abundance alone suggests that headwaters could contribute substantially to γ diversity in streams. Furthermore, the relative spatial isolation of headwaters in the tips of stream networks is expected to enhance β diversity at both community and population-genetic levels (Finn and Poff 2005, 2011, Muneepeerakul et al. 2008, Hughes et al. 2009, Brown and Swan 2010). Headwaters also may exhibit greater among-stream differences in local habitat characteristics than do larger streams, a condition that is expected to promote β diversity (Lowe and Likens 2005, Meyer et al. 2007, Clarke et al. 2008).
A hypothesis of greater β diversity in headwaters than in mid-order streams has not been addressed directly, but some evidence supports the idea within limited biogeographical regions (Finn and Poff 2005, Hughes 2007). This pattern is opposite the pattern of α diversity expected along a linear gradient (Fig. 1), and it potentially exists at both community and population-genetic levels. If headwaters in general tend to be significantly more β-diverse than mid-order streams, a paradigm shift will become necessary in stream ecology, where headwaters have long been considered depauperate.
Fig.1.
Hypothesized inverse patterns of α and β diversity with increasing stream order from 1st to 4th. The α pattern is predicted by the River Continuum Concept (Vannote et al. 1980) and is found in many high-gradient streams. The β pattern is the hypothesis tested in our paper. Both relationships are presented as linear for simplicity, but the true nature of either relationship might be otherwise and probably varies among stream types.
The outcome hypothesized in Fig. 1 also would have significant practical implications. Greater β diversity means that each locality (e.g., each branch in a stream network) contributes a larger proportion of γ diversity. A major implication of this pattern for conservation is that each small branch lost results in a larger proportion of regional biodiversity lost. Many ongoing threats are arguably more insidious in headwaters than in larger streams. These threats include various effects of climate change (Brown et al. 2007, Finn et al. 2009, 2010) and loss of regulatory protection (Meyer et al. 2007, Leibowitz et al. 2008), which encourages more rapid degradation of small streams via filling, pumping, complete removal (Palmer et al. 2010), and other anthropogenic factors.
The hypothesis of greater β diversity in headwaters than in larger streams has been proposed (Clarke et al. 2008, Hughes et al. 2009), but it has not been tested for generality across a range of stream types. Furthermore, β diversity can be conceptualized and measured in many ways (see Anderson et al. 2011 for a comprehensive discussion), and how its interpretation might affect a test of this hypothesis is unclear. Throughout our paper, we approach β diversity as describing general variation in composition (genetic or taxonomic) among local sites within a given habitat class (i.e., stream size categories) (Anderson et al. 2011). We accumulated pre-existing data collected in stream networks from around the world to address the overarching question: Does general evidence exist across major ecozones and stream types that β diversity (quantified both at population-genetic and, to a lesser extent, community levels) is greater among headwaters than among mid-order streams? As an important subcomponent of this question, we applied some alternative common approaches from the fields of population genetics and community ecology to quantify β diversity. We also implemented a conservation-minded simulation that describes β diversity as an expected average rate of erosion of γ as local sites are lost from a region, and we compared these simulation results between the 2 stream-size categories.
Methods
General approach
For both population-genetic and community levels of diversity, we limited our analyses to the macroinvertebrate component of the stream biota. Invertebrates make up the bulk of macrofaunal diversity in streams, and they are well studied representatives of stream diversity. In many cases, the diversity patterns of invertebrates—particularly insects—are likely to represent the consequences of natural movement patterns, whereas other biota (e.g., fish) often include many introduced taxa. We also limited our analyses to nonamphidromous macroinvertebrates in an attempt to keep dispersal characteristics broadly similar among otherwise quite diverse taxa.
Across the globe, α and γ diversity in streams can be expected to vary considerably with biogeographic histories and differing major environmental pressures (Vinson and Hawkins 2003, Clarke et al. 2008, Bonada et al. 2009). Cross-region comparisons of α or γ are also troublesome because, at this scale, researchers typically must rely on data accumulated from several different sources that vary in methods, timing, intensity, level of taxonomic resolution, and spatial scales of resolution (Vinson and Hawkins 2003, Pecher et al. 2010). However, standardizing data sources to compare β diversity among regions is more feasible. Because β diversity in the general sense (Anderson et al. 2011) measures within-region variation, the main concern is to ensure similar sampling across all local sites within regions. Therefore, our strategy was to accumulate reliable data from headwaters or mid-order streams in many regions, standardize them, compare and evaluate the utility of some commonly applied β diversity measures, and ask whether overall (among-region) means of these estimates differed between the 2 stream-size classes.
We evaluated biodiversity data at both community and population-genetic levels, but we chose to emphasize the population-genetic level. Genetic diversity is an important component of landscape-scale biodiversity (Crandall et al. 2000, Manel et al. 2003, Sork et al. 2010) and is a good choice here for 2 main reasons. First, β diversity estimates generated from putatively neutral population-genetic data are expected to respond primarily to the effects of isolation among local sites in a region. At the community level, assemblages probably respond to a combination of isolation and habitat heterogeneity among streams (Thompson and Townsend 2006, Heino and Mykrä 2008, Brown and Swan 2010). Effects measured at the population-genetic level will address the isolation hypothesis more directly, without a potential confounding effect of habitat differences. Second, an international industry standard (GenBank) exists for filing genetic-sequence data, which makes them easy to search and retrieve (Benson et al. 2005). Such a standard does not exist for community-level diversity data. Therefore, the population-genetic level provides practical feasibility that the community level is lacking to gather comparable data easily among many regions. These data can be sorted readily into highly specific regions of the genome, thereby allowing a further degree of data standardization (e.g., by characterizing diversity at homologous loci across multiple species and regions).
Data specifications: both levels of organization
To filter the myriad pre-existing data sets potentially useful for our project, we formulated a series of rules to identify the population-genetic and community-level data most appropriate to answer our overarching question. Essentially, these strict guidelines served the purpose of maximizing region-to-region and data set-to-data set comparability in terms of sample size, spatial extent, and quality. At the coarsest level, we created 3 priority classes for data selection. The highest priority was for data generated within our own laboratory groups or groups with which we had been directly associated as a student or researcher over the preceding 5 y. We did not require data in this priority class to have been published, but at a minimum they had to have been part of a graduate-level thesis. The 2nd priority class consisted of data generated by our close research colleagues, outside of our own laboratory groups. Any data considered in this or the following class we required to be published in a peer-reviewed journal. The 3rd priority class included data generated by any of our colleagues participating in population or community sessions at the 2009 North American Benthological Society meeting in Grand Rapids, Michigan.
After accumulating potential data sets from the prioritized sources, we applied more-detailed filters (Appendix 1). We wanted to adhere to a strict definition of headwaters and mid-order streams that was consistent across regions (see Clarke et al. 2008), so a key requirement was that geographic coordinates were available for each local collection site in a potential data set, and that high-resolution maps (≥1∶25,000) were available to assign stream order accurately for each site. Following assignment of sites to stream orders, each regional data set was classified as either headwaters (1st–2nd order) or mid-order (3rd–4th order) or was split into 2 separate data sets according to stream size if sampling was sufficient (n ≥ 4 sampling sites/stream-size class) for both size categories. Any collection sites on stream orders >4 were removed from consideration. Henceforth, the term data set refers to a regional collection of samples, either all headwater or all mid-order in size, for which β diversity could be measured.
Another requirement imposed upon each data set was spatial independence of all sites (i.e., no 2 sites could be flow-connected; ver Hoef and Peterson 2010). After we sorted data sets into size classes, we culled any flow-connected sites, typically by randomly selecting one of the nonindependent sites for removal. In a few cases, culling was nonrandom to preserve the maximum number of sites in a data set. For example, if a 2nd-order collection site in a headwater data set was flow-connected to multiple 1st-order sites, the 2nd-order site was removed.
Increasing the geographic extent of sampling can increase β diversity estimates (Nekola and White 1999, Soininen et al. 2007), so we also imposed minimum and maximum spatial extents for each data set (Appendix 1), and we required that all sites in each data set were from a single terrestrial ecoregion (Olson et al. 2001). If a data set fell below the minimum size requirement, we removed it from further consideration. If a data set exceeded the maximum limit, we split it and assessed whether smaller subsets of sites could still meet the remaining requirements (Appendix 1). For example, data sets with sites from 2 ecoregions could be split by ecoregion and retained as 2 separate data sets. If a data set covered >50,000 km2 surface area within an ecoregion, we culled it to achieve a smaller spatial extent while retaining the maximum number of sites.
Diversity estimates generated in our paper are likely to differ to some degree from those reported in original publications given the reorganization of the original data sets described above. We also imposed additional guidelines to standardize population-genetic and community data sets (described below).
Data specifications: population-genetic level
Population-genetic data for stream macroinvertebrate species has accumulated rapidly over the past decade, and a large proportion of these data, worldwide, has been in the form of sequence data from the mitochondrial cytochrome c oxidase subunit I (COI) gene. To achieve robust cross-species and cross-region comparability while maintaining a reasonably large number of species and regions represented, we retained only population-genetic data sets composed of COI sequence data.
Single-species collections for population-genetic analyses typically contain information from orders-of-magnitude fewer individuals than do community collections. In some lines of inquiry (e.g., phylogeography; Pauls et al. 2006) collection of just 1 or 2 individuals/sampling locality can be suitable. However, a robust measure of variability among sites (i.e., β diversity), requires a reasonable estimate of α diversity for each site. Therefore, we required that enough individuals were sampled per site to provide an accurate representation of α diversity. We retained only those data sets with a mean of n ≥ 10 individuals/site, with no site having n < 7. For COI, evidence for stream insects suggests that 7 to 10 individuals are sufficient to capture most of the genetic diversity present in a local population (Monaghan et al. 2009). In our study, sites with insufficient n could be removed from some data sets and the smaller data set retained. Our final list for the population-genetic analyses consisted of 18 data sets: 9 headwaters and 9 mid-order (Table 1).
Table 1.
Population-genetic data sets used in our analyses. World Wildlife Fund (WWF) terrestrial ecoregions were taken from Olson et al. (2001). See Appendix 2 for index numbers. Number of populations is the number of local populations in our study (may be reduced from the number in the original studies to meet data-selection guidelines). Dispersal refers to the maximum capacity for terrestrial, among-stream dispersal. AA = Australasia, AT = Afrotropic, NA = Nearctic, NT = Neotropic, PA = Palearctic.
Analyses: population-genetic level
For each of the 18 data sets, we applied 4 different measures to characterize general β diversity, including 2 traditionally used in population genetics (FST and ΦST) and 2 traditionally used in community ecology (Sørensen's dissimilarity measures with either presence/absence or abundance information). We used the distribution of COI haplotypes as the raw input for each. The 2 categories provide 2 different perspectives on β diversity.
Sørensen's dissimilarity in its simplest form and with presence/absence data can be described as:
where W is the number of entities (haplotypes) shared between 2 samples, and A and B are the total number of haplotypes in each of the 2 samples. The maximum value is 1, which is achieved if the sites share no haplotypes. This metric can be converted easily for use with abundance data (McCune and Mefford 2006). We characterized presence/absence and abundance-based Sørensen's dissimilarity for each data set by calculating the average value across all pairs of sites.
The metrics FST and ΦST take a different approach to β diversity, and both require abundance data because they evaluate frequency differences. A simple way to consider either FST or ΦST for haplotype abundance data collected at 2 sites is as:
where Π represents the mean number of pairwise differences in haplotype (among individuals) either between or within the 2 sites. Like Sørensen's dissimilarity, the maximum value of either FST or ΦST = 1. However, the maximum value is more difficult to reach in this case because FST and ΦST are influenced by both among- and within-site variation (see also Jost 2008). Thus, to have an outcome of FST or ΦST = 1, not only must sites share no haplotypes (as with Sørensen's), but also each site must have no polymorphism (i.e., each site must have just a single haplotype). The difference between FST and ΦST is that the former accounts only for qualitative haplotype differences, whereas the latter incorporates information about evolutionary distance among haplotypes. For each of these 4 general measures of β diversity, we used t-tests to compare means from headwater vs mid-order data sets.
We also were interested in the effect of local habitat loss on γ diversity for each of the 2 stream-size categories. This issue is directly related to β diversity because each local habitat is expected to contribute a larger proportion of γ diversity in regions with greater β diversity. We developed a novel approach to visualize the rate at which γ would be lost if local habitats were removed sequentially from a region. We statistically compared these hypothetical loss rates for headwater vs mid-order data sets. We used regional haplotype richness as a measure of γ diversity.
Achieving this objective involved direct comparisons of γ diversity among diverse data sets, so it required 2 further standardization steps. First, we reduced all original COI sequence data to the length of the shortest segment in the original data sets (307 base pairs; Appendix 2). This step was necessary because longer sequences from the same genomic region yield more haplotypes identified/data set, rendering cross-species estimates of γ richness incomparable. Second, we reduced all data sets to n = 4 sites, the minimum size (Appendix 1), by splitting data sets with n > 4 sites (Table 1, number of populations) into geographic quadrants and retaining from each quadrant the site with most individuals sampled. This step was necessary because more densely sampled regions are more likely to have reached the actual γ value for the region, in a manner similar to reaching the asymptote in species-accumulation curves (Gaston and Blackburn 2000). Equalization of the number of local sites in each region was the most effective way to generate comparable γ-diversity values across the data sets.
For each of the reduced data sets, we simulated random, sequential loss of sites until a single site remained. We used linear regression to calculate the slope of the function relating number of sites remaining to total remaining regional richness (γ diversity). We ran this random-loss procedure 1000 to 10,000× and recalculated the slope term each time. Henceforth, we refer to the average slope across all 1000 to 10,000 random-loss simulations for each data set as the slope of loss. The slope of loss allows visualization of the effects of local habitat loss on regional diversity, and its magnitude is an approximation of β diversity. We used a t-test to assess differences in slope of loss between headwater and mid-order data sets. We also tested for differences between mid-order and headwater Trichoptera data sets separately because this group was the only major taxon represented by ≥2 species in both headwater (n = 4) and mid-order (n = 2) data sets.
Data specifications and analyses: community level
We assembled community-level data sets following General data specifications (above), and we included an additional requirement that data were sampled with consistent methods across all collection sites within a data set (Appendix 1). Consistent methods included samples collected within the same season(s), similar field collection methods (including types of microhabitat sampled), and the same taxonomic resolution applied for identification of specimens at all sites. We also ensured that no taxonomic redundancy existed in any data set (i.e., any taxa identified at multiple, potentially overlapping levels were either aggregated into higher taxonomic groupings or the less-resolved taxa were removed from analysis). Our final list for the community-level analyses consisted of 16 data sets: 10 headwater and 6 mid-order (Table 2).
Table 2.
Community-level data sets used in our analyses, descriptive information, mean Sørensen's dissimilarity, and original references for the data. See Appendix 3 for additional information. World Wildlife Fund (WWF) terrestrial ecoregions were taken from Olson et al. (2001). Number of sites = number of local sites (stream reach scale) in our study (may be reduced from original studies to meet data-selection guidelines and after further reduction of Iberian data sets; see text for details). Total number of taxa = taxon richness across the sites included in our study. AA = Australasia, NA = Nearctic, PA = Palearctic.
We allowed the number of sites/data set to vary to some degree. However, we had particularly extensive data sets from 5 different ecoregions on the Iberian Peninsula. We reduced the more-extensive Iberian data sets by randomly retaining 1 site from each 4th-order basin so that these data sets were more comparable to data sets from the rest of the world (which included 4–10 sites). This procedure resulted in sample sizes of 7 to 10 sites/reduced Iberian data set (Table 2). Therefore, all community data sets had final n = 4 to 10 sites.
We calculated mean Sørensen's dissimilarity (with presence/absence data only) for each community data set, and we used a t-test to assess whether means were different between headwater and mid-order data sets. Headwater and mid-order data sets from the same ecoregion were available for 4 ecoregions (1 Nearctic and 3 Palearctic; Table 2). For these data sets, we applied a paired t-test to test the hypothesis that β diversity (as Sørensen's dissimilarity) shifted in a consistent direction from headwaters to mid-order streams within ecoregions. The Colorado Rockies ecoregion had 2 headwater data sets (Table 2), so we used their mean Sørensen's dissimilarity as the headwater value for the paired t-test.
Results
Population-genetic level
The 18 data sets used for analyses at the population-genetic level spanned 5 terrestrial ecozones (Table 1). Australasia was represented in both headwater (n = 2) and mid-order (n = 6) groups, as were Nearctic (n = 3 and n = 1, respectively) and Palearctic (n = 3 and n = 1, respectively). The Afrotropic was represented by 1 headwater data set, and the Neotropic by 1 mid-order data set. All of the species represented had some capacity for terrestrial dispersal among streams. They generally could be classified into 1 of 2 dispersal subcategories, capable of flight or crawling (Table 1). Many (14 of 18) species were insects, only one of which lacked flight ability. The remaining 4 species were decapod crustaceans (1 crab and 3 shrimp), each capable of some degree of terrestrial crawling. The noninsects represented mid-order streams.
Comparisons of β-diversity measures between headwater and mid-order data sets had varying results (Fig. 2A–D). In general, the traditional population-genetic measures (FST and ΦST; Fig. 2A, B) revealed no difference and the traditional community measures revealed significant differences between headwater and mid-order data sets (Fig. 2C, D). FST (p = 0.67) and ΦST (p = 0.24) did not differ between headwater and mid-order data sets, although mean ΦST tended to indicate greater β diversity in headwaters. In contrast, β diversity measured as Sørensen's dissimilarity for presence/absence and abundance data was significantly greater (p = 0.004 and p = 0.007, respectively) in headwater than in mid-order streams.
Fig.2.
Mean (±1 SE) FST (t-test, p = 0.67) (A), ΦST (t-test, p = 0.24) (B), Sørensen's dissimilarity based on presence/absence (+/−) (t-test, p = 0.004) (C), and abundance (abund.) data (t-test, p = 0.007) (D) in headwater (HW) and mid-order (Mid) population-genetic data sets (Table 3).
We also found significant differences in the slope of loss of haplotypes between headwater and mid-order data sets (Fig. 3A, B). The slope of loss for headwater data sets ranged from 1.8 to 5.8 haplotypes/population lost (mean = 4.1, median = 4.3), and the slope of loss for the mid-order data sets ranged from 0.5 to 5.4 (mean = 1.8, median = 1.5) (Table 3). We obtained similar results when we restricted the analysis to include only Trichoptera (headwater mean = 3.7, mid-order mean = 1.2; t-test, p = 0.005). Mode of dispersal (crawl vs fly) had no effect on slope of loss among the species comprising our population-genetic data sets. The 4 mid-order, noninsect data sets had 4 of the 5 lowest slopes of loss (Table 3). The data sets having the top 6 slope-of-loss values each represented species with flight-capable adults, and 5 of these 6 were headwater data sets. For the mid-order data sets, mean slope of loss for the fliers (n = 5) was not different from that for the crawlers (n = 4; t-test, p = 0.36). Only 1 crawler species was represented in the headwater data sets, and its slope of loss equaled the mean for the group (4 haplotypes/population lost, Abedus herberti; Table 3).
Fig.3.
Mean slopes of loss for headwater (A) and mid-order (B) population-genetic data sets (t-test for difference in mean of slopes between stream-size categories, p = 0.004). Heavy lines show data sets representing species that disperse among streams only by crawling. Thin lines show data sets representing species capable of flight (Table 1).
Table 3.
Five estimates of among-population turnover for each of the 18 population-genetic data sets, listed by species name (see Table 1). +/− = based on presence/absence data.
Community level
Eight of the 16 community data sets (n = 5 headwaters and n = 3 mid-order) were collected via intensive sampling of the Iberian Peninsula (Palearctic) by NB and colleagues. Four of the remaining 5 headwater data sets were Nearctic, with the last headwater data set representing nonIberian Palearctic. Two of the remaining 3 mid-order data sets were Nearctic, with the final one from Australasia. All community data sets contained at least the full insect complement of the macrobenthos (minus Diptera for the Iberian data sets), and many contained all macroinvertebrate taxa (Appendix 3). Globally and within ecozones, the community data sets represented a wide array of stream types, from desert springs to alpine to forested lowland streams.
The final data sets included richness (presence/absence) data only, and mean Sørensen's dissimilarity was significantly greater in headwater than in mid-order data sets (p = 0.01; Fig. 4). Headwater means did not differ (at α = 0.05) from mid-order means within ecoregions in which both were represented (p = 0.10; Fig. 5), but means tended to be higher in headwaters than in mid-order streams. Mean Sørensen's values for the 3 Iberian mid-order data sets were exceptionally high, particularly compared to the other 3 mid-order data sets (Table 2).
Fig.4.
Mean (±1 SE) Sørensen's dissimilarity based on taxon presence/absence data (+/−) for headwater (HW) and mid-order (Mid) community data sets (t-test, p = 0.01; Table 2).
Fig.5.
Pairwise comparisons of mean Sørensen's dissimilarity based on presence/absence (+/−) data for headwater and mid-order stream communities in ecoregions having data sets for both stream sizes (paired t-test, p = 0.10). Two headwater data sets existed for the Colorado Rockies ecoregion: one from alpine (A) streams and one from subalpine (SA) streams. The mean of these 2 data sets was used in the paired t-test. Sc/sd = sclerophyllous/semideciduous forest, NE = northeastern, S = southern.
Discussion
Outcomes at both levels of organization
We hypothesized decreasing β diversity along a stream-size gradient from headwaters to mid-order streams, a pattern that is the inverse of that expected (and typically supported) for α diversity (Fig. 1). Results of our global-scale analysis of pre-existing diversity data at population-genetic and community levels supported our hypothesis when β diversity was characterized strictly as among-site variation (in haplotypes or taxa) within a region. Typical measures used to characterize population-genetic structure (FST and ΦST), which incorporate both within- and among-site variation, did not reveal significant differences between headwater and mid-order streams, although ΦST tended to be higher in headwaters than in mid-order streams, a result suggesting increased evolutionary distance among headwater populations (see further discussion below).
All stream types have in common a branching network structure, so this outcome of greater β diversity in headwaters in our global-scale study suggests that general rules could exist for predicting biodiversity distribution in streams (Rodriguez-Iturbe et al. 2009, Brown et al. 2011), potentially across multiple levels of biological organization. Hence, other major drivers like climate, flow regime, biogeographic history, or even assumed dispersal capacity appear to be outweighed by the simple influence of differences in position within a network (tips vs middle branches, in this case) on patterns of β diversity.
Networks are unique spatial structures in which tips are always more isolated from one another than are interior branches. Therefore, spatial isolation probably is a key mechanism underlying the observed patterns of β diversity. Isolation also might explain the greater effect size for mean Sørensen's dissimilarity in the population-genetic data sets than in the community-level data sets (Fig. 2C, D, Fig. 4; note differences in y-axes). Spatial patterns of putatively neutral population-genetic data are expected to reflect primarily gene flow, which is strongly influenced by spatial isolation (Hughes et al. 2009). Community-level β diversity also can be influenced by isolation (e.g., Thompson and Townsend 2006, Muneepeerakul et al. 2008, Bonada et al. 2009, Brown and Swan 2010), but community structure also responds strongly to habitat drivers (i.e., species sorting; Brown and Swan 2010, Brown et al. 2011). Indeed, habitat heterogeneity could be a secondary mechanism to increase β diversity in headwaters because among-stream heterogeneity might be greater in the tips than interior of stream networks (Gomi et al. 2002, Lowe and Likens 2005, Meyer et al. 2007). These effects should be detectable primarily at the community level. Various statistical approaches (e.g., Borcard et al. 1992, Cottenie 2005, Thompson and Townsend 2006) exist to disentangle the potentially interactive effects of habitat heterogeneity and spatial isolation.
The population-genetic level
Increased spatial isolation probably was the key mechanism behind our observation of increased β diversity among headwater populations, and at its extreme, genetic isolation of populations can lead to allopatric speciation. Indeed, one of the headwater data sets we evaluated (Elporia barnardi; Wishart and Hughes 2003) attained the maximum values for both Sørensen's measures (Table 3) because none of the independent populations that we evaluated shared any haplotypes. This pattern of maximum β diversity might indicate incipient speciation in this blepharicerid occupying isolated headwater streams in South African mountain ranges. Seidel et al. (2009) set out to quantify population-genetic patterns of isolated Gammarus pecos populations in Chihuahuan desert headwaters of western North America and instead found cryptic species in each drainage they sampled. Similarly, Myers et al. (2001) found strong population-genetic structure among headwater caddisfly populations across the Great Basin of North America and recommended that different subbasins should be regarded as different management units. Thus, β diversity is particularly pervasive among headwaters in harsh environments, such as deserts or high altitudes, which presumably impede among-stream gene flow (see also Schultheis et al. 2002, others in Table 1).
The differences between measures traditionally used in population-genetic and in community-level studies (Fig. 2A, B vs Fig. 2C, D) revealed the contrasting approaches used to quantify biological variation in the 2 disciplines. The population-genetic measures (FST and ΦST) did not show significant differences between headwaters and mid-order streams across the diverse data sets we evaluated. These measures account simultaneously for among- and within-site variation and, therefore, are more literally fixation indices rather than measures of among-population differentiation (Jost et al. 2010). Indeed, any putative measure of β diversity that is dependent on α diversity to some degree is not robust for comparing regions with differing α diversity (Jost 2007). Alternatively, Sørensen's dissimilarity measures only among-site variation and, therefore, can be compared directly among disparate regions. Thus, FST and ΦST generally are not suitable for comparing β diversity at either biological level of organization. Sørensen's distances are more appropriate, and other related alternatives have been developed specifically for genetic diversity data (e.g., Jost 2008).
From a conservation perspective, increased β diversity essentially means that individual localities each account for a greater proportion of regional-scale biodiversity. Hence, a prudent management choice might be to emphasize conservation in regions (or sections of stream networks) with high β diversity. Our slope-of-loss approach provided a way to visualize differences in the average contribution of local stream reaches to γ diversity between stream-size categories (Fig. 3A, B). The slope-of-loss values in our analysis can be interpreted as additive measures of β diversity (Lande 1996) because they are directly related to the difference between γ and mean α diversity. However, the key utility of the slope-of-loss approach lies in its ability to visualize the effect of sequential loss of local habitats. The results of the comparison of slope of loss between headwaters and mid-order data sets are qualitatively similar to the comparison of mean Sørensen's dissimilarity measures. The results suggest that effective conservation action in headwaters probably would retain valuable sources of biodiversity at the scale of whole stream basins.
Information on genetic diversity for species occupying both headwaters and mid-order streams within the same ecoregion was missing from our analysis. An important next step in testing whether β diversity is consistently greater in headwaters than in mid-order reaches will be to design studies to compare populations of the same species in mid-order reaches and their headwaters. Tests of the hypothesis (Fig. 1) with individual species occurring from 1st- to 4th-order streams would provide an added element of control that is missing from our current analysis. Several of the species used to represent headwaters in our study probably were headwater specialists (e.g., Wishart and Hughes 2003, Finn et al. 2006, 2007, Múrria and Hughes 2008), whereas many of the mid-order representatives were more widespread species. Decreased dispersal is thought to be a beneficial trait for habitat specialists (Hughes et al. 2009), particularly those that specialize on temporally stable habitat types (Roff 1990), and these differences could provide some explanation for our observations of higher β diversity among headwater data sets. Indeed, headwaters probably inherently harbor more habitat specialists (e.g., Gomi et al. 2002, Meyer et al. 2007) than mid-order streams. If this is the case, population-genetic β diversity would be expected to be greater in headwaters than in mid-order streams for the combined reasons of increased habitat specialization and increased spatial isolation.
The community level
The community-level analyses were more limited than the intensive population-genetic analyses. The number of data sets from headwaters and mid-order streams was comparable between the community- and population-genetic levels, but community data sets were analyzed based only on presence/absence data. Nevertheless, our overall observation of significant differences between headwaters and mid-order data sets (Fig. 4) followed the same pattern revealed at the genetic level.
In contrast, community-level β diversity did not differ significantly between pairs of headwater and mid-order data sets collected within the same ecoregion, although the headwaters tended to have greater values (p = 0.10; Fig. 5). The small sample size of paired data sets within ecoregions (n = 4) and the cluster of 3 of these ecoregions on the Iberian Peninsula resulted in too little power to yield a statistically significant result for this test. Evaluation of Fig. 5 reveals that the 3 Iberian data sets had relatively high values of β diversity in both mid-order and headwater streams that were comparable to values of β diversity in the headwater (but not mid-order) data sets for the Colorado Rockies. In general, Iberian streams are quite diverse (Bonada et al. 2009), and high β diversity might occur in both headwaters and mid-order streams because of an atypical biogeographic history (Soininen et al. 2007, Valladolid et al. 2007).
Our results, although relatively restricted, call for more research within and among other types of ecoregion before strong inference can be made regarding our hypothesis (Fig. 1) at the community level. The high variability among studies and regions in methods and quality of benthic invertebrate community data was daunting. However, our demonstration of an overall pattern of increased β diversity in headwaters vs in mid-order streams (Fig. 4) provides proof-of-concept, and the volume of existing community data suggests that further addressing our overarching question should be feasible in the near future.
Conclusions
Differences in stream size alone explained significant variation in β diversity at community and population-genetic levels in streams across the world, although additional variables, such as altitude, watershed topography, and habitat heterogeneity probably would further increase explanatory power. Therefore, stream network structure appears to influence population-genetic and community structure in similar ways, and future studies that evaluate both biological levels simultaneously should be instructive. Evaluation of species- and genetic-diversity patterns across the same landscape provides increased biogeographic understanding beyond that achieved from either individual level in streams (Bonada et al. 2009, Sei et al. 2009, Finn and Poff 2011) and in other environments (e.g., Cleary et al. 2006, Evanno et al. 2009).
We intentionally excluded large streams (≥5th order) from our analysis for the simple reason that not enough consistently collected macroinvertebrate data exist for analysis of β diversity at these stream sizes. The linear framework of the RCC predicts a decrease in α diversity in large streams relative to mid-order streams, and some evidence supports this prediction (Minshall et al. 1985, Watanabe et al. 2008). However, in their natural state, large streams often produce extensive floodplains, which are complex habitat mosaics that typically support high biodiversity (Ward 1998, Ward et al. 1999, Arscott et al. 2005). These large river systems probably contribute substantially to regional-scale biodiversity in streams. Development of a consistent means to incorporate diversity information from large streams would be useful for testing ideas about broad-scale distribution of diversity throughout entire stream networks.
A major implication of our analyses of β-diversity patterns in small-to-mid-order streams is that aquatic ecologists need to change the way we think about the contribution of headwaters to whole-stream biodiversity. The linear view of streams implied that headwater streams were depauperate and, therefore, potentially more expendable than higher-order streams. A more realistic view of streams as networks demonstrates the opposite pattern. Moreover, recent modeling exercises suggest that headwaters might play an important role in actively sustaining biodiversity across many stream sizes (Morrissey and de Kerckhove 2009). These findings demand that headwaters acquire a more pressing conservation status because each small branch lost represents the loss of unique diversity in a river network.
Acknowledgements
This project was a truly global-scale undertaking. We thank the Rosemary Mackay Fund Committee of the J-NABS editorial board and Pamela Silver for encouraging us to pursue the project and Stuart Bunn at the Australian Rivers Institute for funding. This paper is part of the outcome of the project RICHABUN (Patrones de abundancia y riqueza de insectos acuáticos en un amplio gradiente latitudinal: de las comunidades a las poblaciones, CGL2007-60163/BOS, Spanish Ministry of Education and Science). DSF especially thanks Chris Robinson, Dave Lytle, Stuart Bunn, and Sandy Milner for providing places (on 3 different continents) where she could think and write about ideas associated with this article over a span of 2 y. Comments from Will Clements and 2 anonymous referees substantially improved our paper. We also thank our many colleagues who were so gracious with their time, conversation, or data, especially GUADALMED (El estado ecológico de los ríos mediterráneos del nordeste de la Península Ibérica: Ríos en zonas fuertemente urbanizadas o agrícolas, HID98-0323-C05-01, Spanish Ministry of Education and Science) project members and these individuals: Javier Alba-Tercedor, Christine Albano, Fred Allendorf, Andrew Baker, Olivier Baggiano, Sofie Bernays, Mike Bogan, Chris Burridge, Ben Cook, Christine Engelhardt, Manuel Ferreras-Romero, Than Hitt, Joel Huey, Karen Kubow, Alison McLean, Michael Monaghan, José Luis Moreno, Marilyn Myers, Tim Page, Steffen Pauls, LeRoy Poff, Narcís Prat, Ana Previšiæ, Kathryn Real, Maria Rieradevall, Dan Schmidt, Alicia Schultheis, Makiri Sei, Rick Seidel (in memory), Lisa Shama, Carolina Solà, Ross Thompson, Kozo Watanabe, Marcus Wishart, Asako Yamamuro, and Carmen Zamora-Muñoz.
Literature Cited
Appendices
Appendix 1.
Guidelines for assembling population-genetic and community data sets. Refer to text for more details.Table3.
Appendix 2.
Addendum to Table 1: information on population-genetic data sets listed by species name. World Wildlife Fund (WWF) index number identifies the ecoregions listed by name in Table 1. COI = cytochrome c oxidase subunit I, bp = base pair.Table3.
Appendix 3.
Addendum to Table 2: information on community-level data sets (in the same order as in Table 2). World Wildlife Fund (WWF) index number identifies the ecoregions listed in Table 2. Taxonomic resolution: highest possible = species–family level, as per original publications (with no redundancy). Total no. taxa = taxon richness across the sites included in our study and at the level of resolution listed.Table3.
