University of Southern California Dissertations and Theses

The development of novel measures of landscape diversity in assessing the biotic integrity of lotic communities

(USC Thesis Other)

The development of novel measures of landscape diversity in assessing the biotic integrity of lotic communities

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content The development of novel measures of landscape diversity in assessing the biotic integrity of
lotic communities
By
Ariel Levi Simons
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
Biology (Marine Biology and Biological Oceanography)
December 2020
ii
Dedication
This thesis is dedicated to Stacie Dahl. A person with a deep interest in the world. From their
love and support I was able to complete my doctorate.

iii
Acknowledgements
I would like to acknowledge the love and support of the following people:
My parents, for giving me freedom and the skills to handle it. Stacie Dahl, for always being there
for me. My students, for teaching me. And all of my friends as they help make life worth living.

iv
Table of contents
Dedication ................................................................................................................................................... ii
Acknowledgements ................................................................................................................................. iii
List of Tables ............................................................................................................................................. vii
List of Figures .............................................................................................................................................. x
Abstract ..................................................................................................................................................... xiv
Introduction ................................................................................................................................................. 1
References ......................................................................................................................................................... 5
Chapter 1 ................................................................................................................................................... 10
Abstract ............................................................................................................................................................ 10
Introduction .................................................................................................................................................... 10
Materials and Methods ................................................................................................................................ 15
Scope of data ....................................................................................................................................................................................... 15
Sample group selection ................................................................................................................................................................... 16
Sample acquisitions and classifications .................................................................................................................................. 17
Calculating the CSCI .......................................................................................................................................................................... 17
Land use ................................................................................................................................................................................................ 18
Statistical analysis ............................................................................................................................................................................. 18
Results ............................................................................................................................................................... 20
Trends within diversity measures ............................................................................................................................................. 20
Linear models of diversity and health ...................................................................................................................................... 23
Discussion ........................................................................................................................................................ 27
Modeling watershed health and the importance of single and multi-sample measures of diversity ........... 28
Relationships between measures of watershed health and the physical environment ..................................... 29
𝞯 Diversity decay and evidence of niche differentiation .................................................................................................. 31
Conclusion ....................................................................................................................................................... 32
Acknowledgements ...................................................................................................................................... 33
References ....................................................................................................................................................... 33
Chapter 2 ................................................................................................................................................... 47
Abstract ............................................................................................................................................................ 47
Introduction .................................................................................................................................................... 48
Ecological stress and co-occurrence network topology ................................................................................................... 52
Network size .................................................................................................................................................................................. 52
Connectance ................................................................................................................................................................................... 52
Mean co-occurrence strength ................................................................................................................................................. 53
Modularity ....................................................................................................................................................................................... 53
v
Degree heterogeneity ................................................................................................................................................................. 55
Materials and methods ................................................................................................................................ 56
Sample scope ....................................................................................................................................................................................... 56
Sample acquisitions and classifications .................................................................................................................................. 57
Calculating the CSCI .......................................................................................................................................................................... 57
Land use ................................................................................................................................................................................................ 58
Sample group selection ................................................................................................................................................................... 58
Network construction ..................................................................................................................................................................... 59
Topological measures ...................................................................................................................................................................... 59
Modeled biotic integrity index .................................................................................................................................................... 60
Results ............................................................................................................................................................... 60
Trends in co-occurrence network topology .......................................................................................................................... 61
Linear models of watershed biotic integrity using co-occurrence network topology ........................................ 64
Discussion ........................................................................................................................................................ 67
Synthesis and future directions .................................................................................................................................................. 70
Acknowledgements ...................................................................................................................................... 71
References ....................................................................................................................................................... 71
Chapter 3 ................................................................................................................................................... 94
Abstract ............................................................................................................................................................ 94
Introduction .................................................................................................................................................... 95
Materials and Methods ............................................................................................................................. 101
Sample scope .................................................................................................................................................................................... 101
Acquiring and calculating data for the CSCI ....................................................................................................................... 101
Acquiring samples for metabarcoding .................................................................................................................................. 102
DNA metabarcoding and bioinformatics .............................................................................................................................. 103
Land use ............................................................................................................................................................................................. 105
Sample group selection ................................................................................................................................................................ 105
Calculations of 𝞯 diversity .......................................................................................................................................................... 106
Modeled indices of biotic integrity ......................................................................................................................................... 106
Results ............................................................................................................................................................ 107
Sequencing statistics for metabarcoded amplicons ........................................................................................................ 109
Linear models of diversity and biotic integrity ................................................................................................................. 110
𝞯 diversity patterns and the environment .......................................................................................................................... 111
𝞯 diversity and land use ............................................................................................................................................................... 111
𝞯 diversity and altitude ................................................................................................................................................................ 113
𝞯 diversity and distance .............................................................................................................................................................. 114
𝞯 diversity and models of community assembly .............................................................................................................. 115
The relative importance of diversity measures in our modeled indices ................................................................ 116
The relative importance of environmental measures in our modeled indices ................................................... 117
Discussion ..................................................................................................................................................... 119
Metabarcoding and limitations with taxonomic assignments .................................................................................... 120
vi
The importance of local and landscape diversity in modeling regional measures of stream community
health ................................................................................................................................................................................................... 121
Measures of regional stream health and environmental variables .......................................................................... 122
𝞯 diversity decay and models of community assembly ................................................................................................. 128
Conclusion .................................................................................................................................................... 129
References .................................................................................................................................................... 131
Conclusion ............................................................................................................................................... 156
Appendix A: supplementary material ............................................................................................ 159
Chapter 1: supplementary material ..................................................................................................... 160

vii
List of Tables
Chapter Table
number
Table caption
1 1 The relative importance of various diversity parameters in our
modeled stream health index.
1 2 The relative importance of altitude, land use, and distance in
describing linear models of the mean value of the California
Streams Condition Index (CSCI) and modeled index per HUC 8
watershed.
2 1 Topological measures of co‐occurrence networks, their ecological
relevance, and predicted relationship with an increase in stress due
to upstream land use: network size, connectance, mean co‐
occurrence strength, modularity, and degree heterogeneity.
2 2 The relative importance of the topological measures used in our
modeled stream health indices (p < 10
–4
).
2 3 Coefficients of sample site altitude and land use in linear models
describing linear models of the percent of genera of BMIs per
sample site per functional feeding group (All p < 10
–4
unless
otherwise noted).
2 4 The relative importance of altitude, standard deviation on altitude,
land use, standard deviation of land use, and distance in describing
viii
our linear models of the mean value of the CSCI and modeled index
per HUC 8 watershed (All p < 10
–4
unless otherwise noted).
3 1 Number of sequences used from each primer for each community.
3 2 10-fold cross-validated Pearson correlations between mean and
modeled CSCI using communities aggregated to various taxonomic
levels (p < 10
-4
).
3 3 AIC scores of Exponential (Power-law) models of zeta diversity
decline with sampling order (p < 10
-4
).
3 4 The relative importance of zeta diversity parameters in describing
linear models of the mean value of the CSCI with communities
aggregated to various taxonomic levels.
3 5 The relative importance of land use, altitude, and distance in
describing linear models for modeled and (mean) CSCI values using
communities aggregated to various taxonomic levels.

Appendix A S1 Linear models of taxonomic diversity measures with land use,
altitude, and distance.
Appendix A S2 Linear models of functional diversity measures with land use,
altitude, and distance.
Appendix A S3 Linear models of taxonomic β and ζ10 diversity measures with land
use and altitude.
Appendix A S4 Linear models of taxonomic β and ζ10 diversity measures with land
use and altitude.
ix

x
List of Figures
Chapter Figure
number
Figure caption
1 1 An illustration of the first 3 orders of ζ diversity. The mean number
of unique categories of organisms per sample, a measure of α
diversity, is represented by the value ζ1. In comparing two or more
samples the average number of unique categories held in common
between any two samples is represented by the value ζ2. The mean
value then of β diversity, as described by Jaccard distance (Jaccard,
1900), for sets of two or more samples is ζ2/(2 . ζ1 - ζ2) (Hui and
McGeoch, 2014). For sets of three or more samples the value of ζ3
represents all of the unique categories held in common between
three samples. This process can be extended to N samples, allowing
for a determination of the values ζ1 through ζN.
1 2 The mean values per HUC 8 watershed of αtaxa (a), αfunctional (b), βtaxa
(c), βfunctional (d), ζ10,taxa (e), ζ10,functional (f), CSCI (g), and modeled
health index (h).
1 3 Trends in the mean values per HUC 8 watershed of αtaxa, βtaxa, and
ζ10,taxa versus changes in altitude, land use, and distance. Note that a
rise in β or a decline in ζ10 diversity indicates a greater level of
dissimilarity between samples within a group.
xi
1 4 A comparison of the mean CSCI and modeled index with altitude,
land use, and distance with linear best fit lines. Both indices show
similar behavior in response to changes in altitude, land use, and
distance.
2 1 An example of a stressor reducing both the taxonomic richness of 3
communities, from an initial state (𝞪1, 𝞪2, 𝞪3) to a degraded state
(𝞪1’, 𝞪2’, 𝞪3’), as well as the number of unique categories of taxa
held in common between communities.
2 2 A comparison of the first modeled CSCI and mean CSCI colored by
land use (r = .81, p < 10
–4
). CSCI, California Stream Condition
Index.
2 3 A comparison of the second modeled CSCI and mean CSCI colored
by land use (r = .61, p < 10
–4
). CSCI, California Stream Condition
Index.
3 1 An illustration of the first four orders of zeta diversity. The mean
number of unique categories of organisms per sample, a measure of
alpha diversity, is represented by the value 𝞯1. In comparing two or
more samples the average number of unique categories of organism
held in common between any two samples is represented by the
value 𝞯2. This process can be extended to N samples, allowing for a
determination of the values 𝞯1 through 𝞯N.
xii
3 2 (A) The fraction of ASVs from each primer and community, and the
number of ASVs resolved to various taxonomic levels. (B) The
number of unique groups resolved to various taxonomic levels.
3 3 Pearson correlations between environmental variables, measures of
zeta diversity, mean CSCI scores, and modeled indices using
communities of diatoms classified using the rbcL primer and
aggregated to various taxonomic levels (p < 10
-4
).
3 4 Pearson correlations between environmental variables, measures of
zeta diversity, mean CSCI scores, and modeled indices using BMI
communities classified using the 18S-V9 primer and aggregated to
various taxonomic levels (p < 10
-4
).
3 5 Pearson correlations between environmental variables, measures of
zeta diversity, mean CSCI scores, and modeled indices using
prokaryotic communities classified using the 16S 515F-926R primer
and aggregated to various taxonomic levels (p < 10
-4
).
3 6 Pearson correlations between environmental variables, measures of
zeta diversity, mean CSCI scores, and modeled indices using diatom
communities classified using the 18S-V9 primer and aggregated to
various taxonomic levels (p < 10
-4
).
3 7 Pearson correlations between environmental variables, measures of
zeta diversity, mean CSCI scores, and modeled indices using soft
algal communities classified using the 18S-V9 primer and
aggregated to various taxonomic levels (p < 10
-4
).
xiii
3 8 Pearson correlations between environmental variables, measures of
zeta diversity, mean CSCI scores, and modeled indices using
prokaryotic communities classified using the 16S-V4 primer and
aggregated to various taxonomic levels (p < 10
-4
).
Appendix A S1 Trends in the mean values per HUC 8 watershed of αfunctional,
βfunctional, and ζ10,functional versus changes in altitude, land use, and
distance. Note that a rise in β or a decline in ζ10 diversity indicates a
greater level of dissimilarity between samples within a group.

xiv
Abstract
Freshwater streams are of great ecological importance, helping to regulate a variety of
biogeochemical cycles. While the ecosystem services provided by streams help to support a
number of human activities, they are frequently degraded by a number of anthropogenic
stressors. Efforts to monitor the ecological health of stream ecosystems have frequently focused
on measurements based on the local diversity of benthic macroinvertebrates (BMIs) sorted by
morphology. However, individual BMI communities are effectively linked across a landscape
through processes such as dispersal, and as such there is a need to measure the health of multiple
communities on a regional basis. Here we investigate the potential role of three techniques in
development of such regional assessments: zeta diversity, topological measures of co-occurrence
networks, and environmental DNA (eDNA). Through our investigations we were able to
demonstrate the following key findings: (1) patterns of landscape diversity, as described by
either zeta diversity or topological measures of co-occurrence networks, can account for a
significant proportion in the regional health of watersheds as described by the California Streams
Condition Index (CSCI), (2) eDNA is a viable alternative to classifying and organizing
biological communities found in freshwater streams. Our results indicate the potential for the
development of regional scale assessments of the health of freshwater streams which can be
derived novel patterns of landscape diversity in biological communities classified via eDNA.
1
Introduction
Globally, as human populations and their utilization of natural capital has grown so has
the level of anthropogenic stress upon a diverse array of ecosystems (Best, 2019; He and
Silliman, 2019). Lotic ecosystems, that is the biological communities within flowing bodies of
water such as rivers and streams, are particularly sensitive to the presence of these stressors
ranging from changes in land use (Rodríguez-Romero et al., 2018) to climate (Höök et al., 2019).
In spite of these anthropogenic stresses lotic ecosystems still help to mediate a diverse array of
ecosystem services which support human activities, ranging from maintenance of water quality
(Sweeney et al., 2004) to sustaining healthy fisheries (Colvin et al., 2019). Their importance has
led then to the development of a number of tools for assessing the biological integrity of lotic
communities, both as a means of assessing current human impacts as well as motivating efforts
at ecological restoration (Palmer et al., 2014). Monitoring of the biological condition of these
ecosystems has traditionally focused on measurements of the composition of particular
communities, such as benthic macroinvertebrates (BMIs) or diatoms, at a site and comparing it to
their relatively undisturbed counterparts (Fetscher et al., 2014; Beck et al., 2019).
This approach to the assessment of lotic ecosystems, while in widespread use, has to
contend with two significant limitations. First, community compositions are determined by
sorting sampled taxa by their morphology. While such methods are well established they are
both labor intensive and often limited in their level of taxonomic resolution (Jones, 2008;
Manoylov, 2014). Second, bioassessments focus on the composition of individual communities.
This assumes that the compositions of lotic communities are due solely to changes in their
physical environment, although evidence suggests a strong role for dispersal in shaping the
2
distribution of unique organisms across a landscape (Leibold et al., 2004; Brown et al., 2011).
Addressing these limitations has involved the search for alternative methods of classifying the
compositions of lotic communities, as well as the development of bioassessment frameworks
which take into account ecological processes which tie individual communities into regional
scale metacommunities.
One promising avenue for overcoming the limitations of morphologically based methods
of determining the composition of lotic communities involves the use of environmental DNA
(eDNA). Rapid developments in DNA sequencing and computational technology have greatly
reduced the costs and efforts associated with this method. Use of eDNA can enable finer levels
of taxonomic resolution in biological communities than those derived from morphologically
based methods (Pawlowski et al., 2019), and it can allow for the structure of ecologically
important communities of organisms such as bacteria and archaea which would otherwise be
overlooked (Thomsen and Willerslev, 2015).
With use of eDNA bioassessments of lotic communities can ‘zoom’ in on individual
communities in greater detail, potentially capturing the compositions of all of constituent groups
from microbes to metazoa. To ‘zoom’ out and assess lotic communities as components of a
regional whole we turned two novel methods in landscape ecology used to study
metacommunities: zeta diversity and co-occurrence networks.
Zeta diversity was developed as a generalized extension of established measures of
landscape diversity, which could allow comparisons to be made be large numbers of
communities simultaneously (Hui and McGeoch, 2014). This was motivated as beta diversity, a
well-established tool in landscape ecology for measuring pairwise differences in the composition
of communities, has been found to be heavily biased towards the presence or absence of rare
3
taxonomic groups in a landscape (Latombe et al., 2017). By quantifying the compositional
overlaps across any number of communities, zeta diversity could then allow for measurements of
the gain or loss of relatively rare or common taxonomic groups by simultaneously comparing
small or large groups of communities. In the course of developing zeta diversity, it was also
found to enable the investigation of a variety of other ecological questions ranging from
processes of community assembly (McGeoch et al., 2019), dispersal limitation (Leihy et al.,
2018), and the effects of environmental gradients on the structure of metacommunities (Hui et
al., 2018).
While zeta diversity can help investigate a variety of ecological questions, using patterns
of similarity between communities in a landscape, it is somewhat limited in inferring patterns of
co-occurrence between unique taxonomic groups. In landscape ecology patterns of co-
occurrence, specifically spatial co-occurrence, have been used to test a variety of hypotheses
ranging from determining the existence of ecological guilds (Heino and Grönroos, 2013) to
quantifying changes in species interactions following an environmental disturbance (Kafley et
al., 2019). As a whole, these patterns of co-occurrence in a metacommunity can be organized
into networks. The topology of these networks, that is the structural arrangement of their
constituent co-occurrences, have been used to quantify changes in biological communities across
a variety of environmental gradients (Mandakovic et al., 2018; Tavella and Cagnolo, 2019).
Using these tools, we sought to develop new methods in bioassessment for lotic
communities, and in doing so address the following questions:
1. Are there patterns in structure of morphologically classified lotic BMI communities
which can be used as reliable indicators of their biotic integrity on a regional scale?
4
2. Are there patterns in structure of various types of lotic communities, classified from
eDNA, which can be used as reliable indicators of their biotic integrity on a regional
scale?
Chapter 1, titled “Using alpha, beta, and zeta diversity in describing the health of stream‐
based benthic macroinvertebrate communities”, discusses the development of a regional scale
assessment of biotic integrity for BMIs using linear models of various measures of zeta diversity.
To build our models we used approximately 4000 morphologically sorted communities of lotic
BMIs, sampled over decades across California, and compared our modeled measure of biotic
integrity against the established measure of the California Streams Condition Index (CSCI). We
found that, on a regional scale, our model could account for a majority of the observed variation
in CSCI scores. We also used zeta diversity to find strong evidence for the assembly of lotic BMI
communities to be driven by a process of niche differentiation. This suggested that changes in
landscape diversity alone has the potential to be used as a tool for investigating ecological
patterns in lotic BMI communities, as well as providing for a means of assessing their biotic
integriy.
Chapter 2, titled “Using co‐occurrence network topology in assessing ecological stress in
benthic macroinvertebrate communities”, builds on the theme of investigating the use of novel
measures of landscape diversity for bioassessments. Starting with the same data used in Chapter
1, we generated co-occurrence networks for regional groups of communities and calculated
various measures of their topology (e.g. the density of edges or how evenly there are distributed
across the network). To test how applicable these topological parameters could be in the
bioassessment of lotic BMI communities we constructed a model of CSCI scores those values.
We found a simple linear model of just five topological parameters could account for the
5
majority of the observed variation in CSCI scores. Our results provided further support for the
idea that landscape ecology patterns alone could provide for reliable measures of regional scale
biotic integrity.
Chapter 3, titled “Zeta diversity patterns in metabarcoded lotic communities as a tool for
bioassessment”, returns to the investigation of zeta diversity as tool for bioassessment but with a
focus on multiple types of biological communities found in lotic environments. The size of the
group of communities used in this investigation was far smaller, approximately 150, but in using
metabarcoding instead of morphology to classify organisms they were far more taxonomically
rich. We found support for an approach in synthesizing applications of eDNA and the landscape
ecology of metacommunities in performing bioassessments, but only for particular taxonomic
groups (e.g. BMIs and diatoms).
References
Beck, M. W., Mazor, R. D., Theroux, S., & Schiff, K. C. (2019). The Stream Quality Index: A
multi-indicator tool for enhancing environmental management. Environmental and Sustainability
Indicators, 1, 100004.

Best, J. (2019). Anthropogenic stresses on the world’s big rivers. Nature Geoscience, 12(1), 7-
21.

Brown, B. L., Swan, C. M., Auerbach, D. A., Campbell Grant, E. H., Hitt, N. P., Maloney, K. O.,
& Patrick, C. (2011). Metacommunity theory as a multispecies, multiscale framework for
6
studying the influence of river network structure on riverine communities and
ecosystems. Journal of the North American Benthological Society, 30(1), 310-327.

Colvin, S. A., Sullivan, S. M. P., Shirey, P. D., Colvin, R. W., Winemiller, K. O., Hughes, R. M.,
... & Danehy, R. J. (2019). Headwater streams and wetlands are critical for sustaining fish,
fisheries, and ecosystem services. Fisheries, 44(2), 73-91.

Fetscher, A. E., Stancheva, R., Kociolek, J. P., Sheath, R. G., Stein, E. D., Mazor, R. D., ... &
Busse, L. B. (2014). Development and comparison of stream indices of biotic integrity using
diatoms vs. non-diatom algae vs. a combination. Journal of applied phycology, 26(1), 433-450.

He, Q., & Silliman, B. R. (2019). Climate change, human impacts, and coastal ecosystems in the
Anthropocene. Current Biology, 29(19), R1021-R1035.

Heino, J., & Grönroos, M. (2013). Does environmental heterogeneity affect species co‐
occurrence in ecological guilds across stream macroinvertebrate
metacommunities?. Ecography, 36(8), 926-936.

Höök, T. O., Foley, C. J., Collingsworth, P., Dorworth, L., Fisher, B., Hoverman, J. T., ... &
Tank, J. (2019). An assessment of the potential impacts of climate change on freshwater habitats
and biota of Indiana, USA. Climatic Change, 1-20.

7
Hui, C., & McGeoch, M. A. (2014). Zeta diversity as a concept and metric that unifies incidence-
based biodiversity patterns. The American Naturalist, 184(5), 684-694.

Hui, C., Vermeulen, W., & Durrheim, G. (2018). Quantifying multiple-site compositional
turnover in an Afrotemperate forest, using zeta diversity. Forest Ecosystems, 5(1), 1-9.

Jones, F. C. (2008). Taxonomic sufficiency: the influence of taxonomic resolution on freshwater
bioassessments using benthic macroinvertebrates. Environmental Reviews, 16(NA), 45-69.

Kafley, H., Lamichhane, B. R., Maharjan, R., Khadka, M., Bhattarai, N., & Gompper, M. E.
(2019). Tiger and leopard co-occurrence: intraguild interactions in response to human and
livestock disturbance. Basic and Applied Ecology, 40, 78-89.

Latombe, G., Hui, C., & McGeoch, M. A. (2017). Multi‐site generalised dissimilarity modelling:
using zeta diversity to differentiate drivers of turnover in rare and widespread species. Methods
in Ecology and Evolution, 8(4), 431-442.

Leibold, M. A., Holyoak, M., Mouquet, N., Amarasekare, P., Chase, J. M., Hoopes, M. F., ... &
Loreau, M. (2004). The metacommunity concept: a framework for multi‐scale community
ecology. Ecology letters, 7(7), 601-613.

8
Leihy, R. I., Duffy, G. A., & Chown, S. L. (2018). Species richness and turnover among
indigenous and introduced plants and insects of the Southern Ocean Islands. Ecosphere, 9(7),
e02358.

Mandakovic, D., Rojas, C., Maldonado, J., Latorre, M., Travisany, D., Delage, E., ... & Cabrera,
P. (2018). Structure and co-occurrence patterns in microbial communities under acute
environmental stress reveal ecological factors fostering resilience. Scientific reports, 8(1), 5875.

Manoylov, K. M. (2014). Taxonomic identification of algae (morphological and molecular):
species concepts, methodologies, and their implications for ecological bioassessment. Journal of
phycology, 50(3), 409-424.

McGeoch, M. A., Latombe, G., Andrew, N. R., Nakagawa, S., Nipperess, D. A., Roigé, M., ... &
Steinberg, P. D. (2019). Measuring continuous compositional change using decline and decay in
zeta diversity. Ecology, 100(11), e02832.

Palmer, M. A., Hondula, K. L., & Koch, B. J. (2014). Ecological restoration of streams and
rivers: shifting strategies and shifting goals. Annual Review of Ecology, Evolution, and
Systematics, 45, 247-269.

Pawlowski, J., Kelly-Quinn, M., Altermatt, F., Apothéloz-Perret-Gentil, L., Beja, P., Boggero,
A., ... & Feio, M. J. (2018). The future of biotic indices in the ecogenomic era: Integrating (e)
9
DNA metabarcoding in biological assessment of aquatic ecosystems. Science of the Total
Environment, 637, 1295-1310.

Rodríguez-Romero, A. J., Rico-Sánchez, A. E., Mendoza-Martínez, E., Gómez-Ruiz, A.,
Sedeño-Díaz, J. E., & López-López, E. (2018). Impact of changes of land use on water quality,
from tropical forest to anthropogenic occupation: A multivariate approach. Water, 10(11), 1518.

Sweeney, B. W., Bott, T. L., Jackson, J. K., Kaplan, L. A., Newbold, J. D., Standley, L. J., ... &
Horwitz, R. J. (2004). Riparian deforestation, stream narrowing, and loss of stream ecosystem
services. Proceedings of the National Academy of Sciences, 101(39), 14132-14137.

Tavella, J., & Cagnolo, L. (2019). Does fire disturbance affect ant community structure? Insights
from spatial co-occurrence networks. Oecologia, 189(2), 475-486.

Thomsen, P. F., & Willerslev, E. (2015). Environmental DNA–An emerging tool in conservation
for monitoring past and present biodiversity. Biological conservation, 183, 4-18.

10
Chapter 1
Using alpha, beta, and zeta diversity in describing the health of stream-based benthic
macroinvertebrate communities
Abstract
Ecological monitoring of streams has frequently focused on measures describing the taxonomic,
and sometimes functional, α diversity of benthic macroinvertebrates (BMIs) within a single
sampled community. However, as many ecological processes effectively link BMI stream
communities there is a need to describe groups of communities using measures of regional
diversity. Here we demonstrate a role for incorporating both a traditional pairwise measure of
community turnover, β diversity, in assessing community health as well as ζ diversity, a more
generalized framework for describing similarity between multiple communities. Using 4395
samples of BMI stream communities in California, we constructed a model using measures of α,
β, and ζ diversity which accounted for 71.7% of among-watershed variation in the mean health
of communities, as described by the California Streams Condition Index (CSCI). We also
investigated the use of ζ diversity in assessing models of stochastic versus niche assembly across
communities of BMIs within watersheds, with the niche assembly model found to be the likelier
of the two.
Introduction
Stream ecosystems provide a number of ecological services for human activities, ranging
from supporting populations of native species to moderating the delivery of nutrients to coastal
11
ecosystems (Anderson et al., 2002; Dudgeon et al., 2006). It has also been recognized that
monitoring the health of aquatic stream communities is of importance not just to the functioning
of the immediate environment, but to downstream ecosystems and communities impacted by
events such as harmful algal blooms (Carpenter et al., 1998; Paerl et al., 2016). However, given
the close proximity of many riparian systems to agricultural and urbanized land, such
communities are frequently stressed by human activity in the surrounding watershed (Urban et
al., 2006; Tonkin et al., 2016). Such stress can take a variety of forms depending on the
geographic scope and intensity of land use, such as increased nutrient levels from agricultural
runoff or increases in sediment load from upstream erosion (Sala et al., 2000; Allan, 2004;
Dudgeon et al., 2006, Abell et al., 2007, Strayer and Dudgeon, 2010). Ecological stress has been
found to be associated with a decline in local (α) stream community diversity (Hendrickx et al.,
2007; Petrin et al., 2008; Heino, 2009). Local stream communities are not isolated, being
connected through mechanisms such as dispersal, and may more accurately be considered part of
a regional-scale metacommunity (Grönroos et al., 2013; Socolar et al., 2016). There is
importance then in considering both local and regional measures of biodiversity (Socolar et al.,
2016), as well as regional gradients in environmental and anthropogenic activity, in shaping
community structure across watersheds (Shade et al., 2008; Rawi et al., 2013; Brendonck et al.,
2015; Tonkin et al., 2016).
One approach to assessing the impacts of environmental gradients, both natural and
anthropogenic, across a set of stream communities is to consider the regional-scale (γ) diversity
of the communities in question (Heino, 2009). γ diversity has the advantage of being
straightforward to calculate, being the product of the average local-scale diversity (α) and
pairwise similarity between local communities (β) (Whittaker, 1972). Measures of diversity
12
within (i.e., α), and between or among (i.e., β and γ) communities, have been commonly used in
describing stream communities (Vander Vorste et al., 2017). However, such measures are limited
in fully describing diversity in groups of more than two samples (Chao et al. 2008). For example,
estimates of β diversity are biased towards the turnover of rare types of organisms (Latombe et
al., 2017b).
To address this limitation, we incorporated a new measure, zeta (ζ) diversity, in our study
of a set of California stream communities. ζ diversity describes the overlap in categories shared
between multiple samples (Figure 1) (Hui and McGeoch, 2014). In our study these categories are
defined either by membership of unique taxa or of a particular functional feeding group. Using
this framework, one can describe the mean degree of overlap in unique categories among N
samples (ζN) across a given set of samples. From this metric one can then define the mean α
diversity, as defined by categorical richness, across a sample set as ζ1, while the mean β
diversity, as defined by Jaccard’s similarity index (Jaccard, 1900), can then be defined as ζ2/(2 •
ζ1 - ζ2) (Hui and McGeoch, 2014). What is novel is that the ζ diversity framework enables
measurement of diversity involving three or more sampled communities which cannot be
determined using solely α and β diversity. For example, the loss of categories of taxa found
across most samples in an environment can be a useful indicator of widespread environmental
degradation (Gaston and Fuller, 2008; Pond, 2012).
13

Here we propose the utility of ζ diversity in describing variations in the composition of
BMI communities in streams over the entire state of California. We demonstrate the role of a
number of factors in shaping ζ diversity, such as variation in the natural landscape (sample site
altitude), as well as levels of human impact (land use upstream of sample sites). Using samples
(Figure 2) collected by the State of California’s Surface Water Ambient Monitoring Program
(SWAMP), we also investigated trends in how ζ diversity declines with the number of sampled
communities (namely, ζ order). How ζ diversity declines with ζ order can allow for a test of
likelihood for two models of community assembly, either a stochastic process (Munoz et al.,
2008) or one driven by a niche differentiation process (Scheiner et al., 2011). Given prior
observations about how local environment (Siqueira et al., 2012; Astorga et al., 2014) and
predation (Chase et al., 2009) have shaped BMI community structures our expectation was to
Sample A

Sample A
Sample B

Sample A
Sample B
Sample C

Figure 1: An illustration of the first 3 orders of ζ diversity. The mean number of unique categories of organisms per
sample, a measure of α diversity, is represented by the value ζ 1. In comparing two or more samples the average
number of unique categories held in common between any two samples is represented by the value ζ 2. The mean
value then of β diversity, as described by Jaccard distance (Jaccard, 1900), for sets of two or more samples is ζ 2/(2 .
ζ 1 - ζ 2) (Hui and McGeoch, 2014). For sets of three or more samples the value of ζ 3 represents all of the unique
categories held in common between three samples. This process can be extended to N samples, allowing for a
determination of the values ζ 1 through ζ N.
14
observe dominance of niche differentiation processes over stochastic ones. We then
demonstrated how various measures of ζ diversity can be used to describe the health of
communities of BMIs as described by the California Streams Condition Index (CSCI) (Rehn et
al., 2015; Mazor et al., 2016). ζ diversity, as a generalized extension of α and β diversity, has the
potential to both better illustrate, at the watershed scale, trends in turnover as well as a more
general framework for assessing the health of groups of BMI communities.
15

Materials and
Methods
Scope of data
The initial scope of
data covered in this
analysis consists of 4984
stream samples from 2997
unique geographic
locations across the state of
California, constituting a
23-year period (1994-
2016). Every sample
contains the following
data: BMIs enumerated and
sorted to a standardized
level (generally a genus-
level identification except
chironomids which were
identified to subfamily;
Richards and Rogers,
(a) (b)
(c) (d)
(e) (f)
(g) (h)
Figure 2: The mean values per HUC 8 watershed of α taxa (a), α functional (b), β taxa (c),
β functional (d), ζ 10,taxa (e), ζ 10,functional (f), CSCI (g), and modeled health index (h).
16
2011), a bioassessment index score based on a composite of taxonomic and functional diversity
within BMI assemblages known as the California Stream Condition Index (CSCI), sample site
altitude in meters, and the percent developed land use (agricultural, urban, and managed
landscapes) within a 5 km clipped buffer of the watershed upstream of the sampling site. We
used the CSCI as our measure of community health as it has been extensively validated across
the state of California, and which has been used as the primary bioassessment tool used to assess
beneficial uses in California and regulatory compliance with National Pollutant Discharge
Elimination System (NPDES) programs (Mazor et al., 2016). Sites were assigned to hydrologic
regions (hereafter referred to as “watersheds”) at different regional scales defined by the U.S.
Geological Survey Hydrologic Unit Code (HUC), a standardized watershed classification system
developed by the U.S. Geological Survey that organizes watersheds in a nested hierarchy by size
(Seaber et al., 1987). Each HUC is assigned a geographically unique numerical ID, with
geographic scale decreasing in increments of two digits from a regional (HUC 2) to
subwatershed (HUC 12) scale. Our initial data contained 118 HUC 8, 22 HUC 6, 14 HUC 4, and
3 HUC 2 watersheds. We designated HUC 8 watersheds as the minimum geographic scale for
investigating patterns of diversity.
Sample group selection
To investigate patterns in our measures of α, β, and ζ diversity we first constructed
equally sized groups of samples from within each HUC 8 watershed. We chose watersheds with
at least 25 unique samples as having sufficient data density for meaningful analysis. This
filtering reduced our overall data set from 4984 to 4395 unique samples. We then randomly
subsampled 25 stream samples within each of the remaining 56 HUC 8 watersheds. For analysis
of trends in sample diversity this process was repeated 100 times. For each group of 25 samples
17
we calculated the mean sample site altitude (altitude), the mean percent developed upstream land
use (land use), and for mean geographic separation distance in meters between samples
(distance) we used the distm function within the R package geosphere (Hijmans et al., 2017).
Sample acquisitions and classifications
Slightly more than half (55%) of the BMI samples were collected following the reach-
wide protocol of (Peck et al., 2006); the other samples were collected with targeted riffle
protocols, which produce comparable data (Gerth and Herlihy, 2006; Herbst and Silldorff, 2006;
Rehn et al., 2007). For most of the available data, taxa were identified to genus, although this
level of effort (as well as the total number of organisms per sample) varied among samples,
necessitating standardization of BMI data. Membership of a particular taxon to a functional
feeding group was done using CAMLnet (Ode, 2003). Using these methods, we assigned all
identified BMIs to one of 334 unique taxonomic groups and 8 unique functional feeding groups.
Calculating the CSCI
The CSCI is a predictive measure of community health at a given site which compares
observed taxa and metrics to values expected under reference conditions based on site-specific
landscape-scale environmental variables, such as watershed area, geology, and climate (Mazor et
al., 2016). This index comprises two sets of measurements using a standardized taxonomy for
BMI communities (Richards and Rogers, 2011): the first being a ratio of observed-to-expected
taxa (O/E), and the second a predictive multi-metric index (pMMI) made of 6 metrics related to
ecological structure and function of the BMI assemblage describing the site’s community. The
CSCI, and its components, were designed to have minimal influence from major natural
18
gradients. Therefore, they can be used as a measure of biological conditions with a consistent
meaning in different environmental settings (Reynoldson et al. 1997, Hawkins et al. 2010).
Land use
The type and geographic extent of land use in the upstream vicinity of each sampling site
data is derived from the National Land Cover Data set (Homer et al., 2007). The type of land use
attributed to human activity is measured by the total percent of land cover in a designated area
dedicated to agriculture, urbanization, or otherwise managed vegetative landscapes such as golf
courses. The sum of all of these measures is applied to each sample site using a 5 km watershed-
clipped buffer upstream of a stream sampling site using ArcGIS tools (version 10.3;
Environmental Systems Research Institute, Redlands, California) (Mazor et al., 2016). Land use
data are derived for all sample sites using satellite measurements acquired in the year 2000. It
should be noted that the samples in our study were located in areas where the percent developed
land use did not change significantly over the range of time covered in this study (r
2
= -0.02, p =
n.s.).
Statistical analysis
The diversity measures we used in this study, separately calculated for samples
categorized both by taxonomic groups as well as functional feeding groups, were mean α
diversity, β diversity, the ζ diversity decline power-law exponent (b), and the mean number of
groups held in common across 10 samples (ζ10). Both measures of α and ζ diversity were
calculated using the R package ZETADIV (Latombe et al., 2017a), with the mean α diversity per
group of samples calculated as ζ1. The mean value of β diversity, as determined from the Jaccard
distance, was calculated per group of samples using the vegdist function in the R package vegan
19
(Oksanen et al., 2013). To then determine the contribution of each of our experimental factors to
our observed variation in β diversity we performed a 1000 permutation PERMANOVA using the
adonis function in the R package vegan (Oksanen et al., 2013).
To investigate how similar on average our samples were within groups, as described by ζ
diversity, we then tested two models, a power-law of the form ζN = ζ1 . N
-b
and an exponential of
the form ζN = ζ1 . e
b.(N-1)
, in describing how overlaps between samples declines with ζ order. The
ζ diversity decline exponent, b, is a measure of the dissimilarity within a sample set with a
greater value of b indicating a greater overall dissimilarity between constituent samples across a
multi-site assemblage. We fitted these models to the first 10 orders of ζ within a given set of 25
stream samples. The first 10 orders of ζ diversity, ζ10, were chosen as this was sufficient to show
an asymptotic ‘floor’ for the average number of unique taxa or functional feeding group held in
common across a sample set. Support for each model was determined using an Akaike
Information Criterion (AIC) score within ZETADIV.
For each set of 25 samples the AIC score was found to be significantly less for the
power-law model than the exponential decline model for community compositions described
either by taxonomic (p < 10
-4
) or functional feeding groups (p < 10
-4
), as measured by a
Wilxocon signed-rank test. As a further check on the support for a power-law model describing ζ
diversity decline, we repeated this analysis using the same sample groups, but with count data
randomly re-assigned to the taxonomic and functional groups present in those same sample
groups. Repeating a comparison of the AIC scores for both models for this randomized data
produced ζ diversity curves more likely to be described by an exponential decay model for both
taxonomic (p < 10
-4
) and functional diversity (p < 10
-4
), which supports the assertion of taxa
being distributed via stochastic process producing an exponential model of ζ diversity decay (Hui
20
and McGeoch, 2014). Subsequently, the power-law model was then selected as the likelier of the
two in describing the decay in ζ diversity with ζ order.
Our measure of α diversity was compared against altitude, land use, and distance using
the lm function in the stats R package (v3.5.1 Team, 2018). This same package was used to
construct a best-fit linear model to predict the mean CSCI score of a set of samples given their
diversity measures. To determine the relative importance of each factor in a linear model the
function calc.relimp was used within the relaimpo R package (Grömping, 2006).
To construct a linear model to best predict the ecological health of a set of samples, as measured
by their mean CSCI value, we first applied a backwards elimination method in selecting
significant diversity parameters (Pearman, 1997; Snodgrass, 1997). In comparing the AIC scores
of each linear model after the removal of a parameter we dropped two parameters, btaxa and
bfunctional, as they made no significant contribution in predicting mean CSCI score. The resulting
linear involves six parameters: αtaxa, αfunctional, βtaxa, βfunctional, ζ10,taxa, and ζ10, functional.
Results
Trends within diversity measures
From our analysis of 100 permutations on 56 within-watershed groups of samples we
observed four general trends between our measures of diversity (αtaxa, αfunctional, βtaxa, βfunctional,
ζ10,taxa, and ζ10, functional) and environmental parameters (altitude, land use, and distance). Both the
mean taxonomic and functional richness per sample have a negative correlation with
environmental stress, described either by a rise in land use or a decline in altitude (Figure 3a and
d; Appendix S2: Figure S1a and d). There is negative correlation between the mean level of
similarity within a group between samples and the mean level of environmental stress (Figure
21
3b, c, e, and f; Figure S1b, c, e, and f). There is also a negative correlation between the mean
level of similarity within a group between samples and the mean geographic distance between
samples (Figure 3 g-I; Appendix 2: Figure S1 g-i). On a watershed scale, coefficients describing
the strength and significance of the relationships between measures of community diversity and
turnover with environmental metrics indicate distance appears not to underlie much of the
variation in our measures of diversity as compared to either altitude or land use (Appendix 1:
Table S1-2).
We then further investigated the contributions of our experimental factors to the observed
variations in β diversity using a PERMANOVA with 1000 permutations. We found altitude
(F(1,4315 = 58.8, p < 10
-4
) and upstream land use (F(1,4315) = 32.8, p < 10
-4
) accounted for
more of the variation in β diversity than either HUC 8 watershed (F(41,4315) = 7.3, p < 10
-4
) or
sampling year (F(22,4315) = 8.4, p < 10
-4
). As a result, our focus was on the role of land use,
distance, and altitude within HUC 8 watersheds in describing patterns of diversity. We also
determined, in line with prior observations (Vander Vorste et al., 2017), that variations in β
diversity tended to increase with geographic scale, with it increasing from (F(6,4315) = 11.9, p <
22
10
-4
) with HUC 6 watersheds, to (F(7,4315) = 38.9, p < 10
-4
) at HUC 4 watersheds, and
(F(1,4315) = 133.4, p < 10
-4
) at HUC 2 watersheds.

Figure 3: Trends in the mean values per HUC 8 watershed of α taxa, β taxa, and ζ 10,taxa versus changes in altitude, land
use, and distance. Note that a rise in β or a decline in ζ 10 diversity indicates a greater level of dissimilarity between
samples within a group.
23
Linear models of diversity and health
Using six parameters, αtaxa, αfunctional, βtaxa, βfunctional, ζ10,taxa, and ζ10, functional, a linear model
was constructed to best predict the mean value of the CSCI score for a set of samples (Table 1).
The functional relationship between the mean CSCI score per sample group and the remaining
diversity parameters are as follows:

MeanCSCI = -0.3 + 6.2 x 10
-3
x αtaxa + 0.2 x αfunctional – 0.6 x βtaxa + 1.0 x βfunctional – 6.5 x 10
-3
x
ζ10,taxa + 1.7 x 10
-2
x ζ10, functional

Table 1: The relative importance of various diversity parameters in our modeled stream health index.
Parameters F(1,5568) P-value Relative importance
(%)
αtaxa 1.2 x 10
4
< 10
-4
32.5
αfunctional 1151 < 10
-4
36.7
βtaxa 21.9 < 10
-4
6.3
βfunctional 690.6 < 10
-4
10.3
ζ10,taxa 19.4 < 10
-4
3.1
ζ10, functional 50.5 < 10
-4
11.2

24
This model was found to explain 71.7% of the observed variation in the mean value of
the CSCI score for a set of samples. This modeled health index was found to vary in accordance
with altitude, land use, and distance for a set of samples in a similar fashion as the mean CSCI
score (Table 2, Figure 4).

25

26

Figure 4: A comparison of the mean CSCI and modeled index with altitude, land use, and distance with linear best fit
lines. Both indices show similar behavior in response to changes in altitude, land use, and distance.
27
Table 2: The relative importance of altitude, land use, and distance in describing linear models of the mean value of
the CSCI and modeled index per HUC 8 watershed.
Health index CSCI Modeled index
Proportion of variation due to altitude F(1,5571) = 9243 F(1,5596) = 3162
P-value, variation due to altitude < 10
-4
< 10
-4

Proportion of variation due to land use F(1,5571) =1.1 x 10
4
F(1,5596) = 6890
P-value, variation due to land use < 10
-4
< 10
-4

Proportion of variation due to distance F(1,5571) = 45.1 F(1,5596) = 12.5
P-value, variation due to distance < 10
-4
< 10
-3
Relative importance of altitude (%) 24.0 15.6
Relative importance of land use (%) 73.1 78.7
Relative importance of distance (%) 2.8 5.7
Proportion of variance explained by model
(%)
78.3 64.3

Discussion
We found 𝞯 diversity, as an extension of established measures of diversity, to play a
potential in role in both the investigation of BMI community assembly processes as well as
assessing the health of groups of BMI communities. In determining a greater likelihood for a
28
power-law versus exponential decay model for how 𝞯 diversity declines with 𝞯 order, we found
BMI communities in California streams are more likely to assembled via a niche differentiation
rather than a stochastic process. While both measures of 𝞪 and 𝞫 diversity were found to
explain sizable proportions of variation in the mean CSCI, our measure of BMI community
health, within each of the 56 watersheds across California we found the addition of measures of
𝞯10 diversity to improve the ability of our modeled index to account for variation in our measure
of BMI community health. As our modeled index uses a more generalized framework composed
of presence/absence diversity measures, rather than relying on a set of undisturbed reference
communities, we believe there is the potential to utilize the 𝞯 diversity framework in creating
community health assessment indices applicable in environments beyond the current geographic
scope of the CSCI.
Modeling watershed health and the importance of single and multi-sample
measures of diversity
We found both the mean watershed values of αtaxa and αfunctional to have the greatest
relative importance in our model to describe the variation in the mean CSCI on a per watershed
basis. Measures of diversity across multiple samples in a watershed (β and ζ) further improved
our model’s predictive powers. It is also noteworthy that both β and ζ10 diversity have a
significant role in our model, although they are potentially capturing two different patterns
describing multiple sampled communities. In general measures of diversity which incorporate
the overlap between large numbers of communities, such as ζ10 diversity, are more sensitive to
the presence of common groups of organisms than those which only involve pairwise
comparisons, such as β diversity, which are inherently more biased to the presence of rare groups
29
(Latombe et al., 2017b). As we observed that both ζ10 and β diversity have a significant role in
our model of the health of BMI communities across a watershed (Table 1) we propose that the
effects of environmental degradation cannot be fully assessed by the loss of rare or common taxa
alone, but both must be considered in tandem. We also noted that a sizeable portion of variation
in the mean value of the CSCI per watershed could be attributed to our modeled index. As this
index is comprised solely of a linear combination of measures α, β, and ζ10 diversity it suggests
future assessments of the health of BMI communities across a watershed could be done using a
more general framework than the one used to devise the CSCI.
Relationships between measures of watershed health and the physical
environment
We found that both our modeled index and the CSCI have similar patterns of variation
with our measurements of the physical environment (Figure 4).Variation in the average αtaxa and
αfunctional diversity per watershed appear to underpin a significant portion of the variation of both
the CSCI, itself a composite measure of taxonomic and functional richness (Mazor et al., 2016),
and our modeled index of watershed health (Table 2). Prior observations support the general
trend of the decline in BMI sample diversity associated with a rise in environmental stress at the
sample site (Paul and Meyer, 2001; Stepenuck et al., 2002; Ourso and Frenzel, 2003; Wallace
and Biastoch, 2016). This has been observed with αtaxa and upstream land use (Sponseller et al.,
2001; Stepenuck et al., 2002; Allan, 2004), αtaxa and altitude (Ward, 1986; Jacobsen et al., 1997;
Jacobsen, 2008), as well as with αfunctional and both altitude and land use (Huryn et al., 2002; Finn
and Poff, 2005; Heino, 2005). Our results then support the utility of using both average αtaxa and
αfunctional diversity per watershed as indicators of the health of their resident BMI communities.
30
We also observed that communities of BMIs tended to have greater degrees of similarity
(lower β and higher ζ10) both taxonomically as well as functionally, within watersheds at a higher
AL with the opposite trends observed in comparison to land use and distance (Figure 4,
Appendix 2: Figure S1, Appendix 1: Table S1-2). Our observed relationship between β diversity
and distance, where a greater geographic distance between samples is associated with greater
inter-sample dissimilarity, is in line with previous observations of communities of stream BMIs
(Astorga et al., 2012; Rouquette et al., 2013). However, our observations appear to run counter to
previous research indicating a positive correlation between the mean altitude of a set of samples
and their degree of community dissimilarity (Finn and Poff, 2005; Frady et al., 2007), as well as
a negative correlation with upstream land use (McGoff et al., 2013).
One potential mechanism for the observed correlations with dissimilarity between
sampled communities within watersheds, negative with altitude and positive with land use, may
be driven by an underlying correlation with habitat heterogeneity. Prior observations support a
positive relationship between the β diversity of multi-site assemblages of BMIs and variations in
the physical environment; habitat heterogeneity (Shostell and Williams 2007; Astorga et al.,
2014; Zhang et al., 2014). Using the standard deviation of land use on a per watershed basis as a
proxy for habitat heterogeneity, we find strong correlations with our measures of dissimilarity
between communities: βtaxa (r = 0.6, p < 10
-4
), βfunctional (r = 0.6, p < 10
-4
), ζ10,taxa (r = -0.4, p < 10
-
4
), and ζ10,functional (r = -0.5, p < 10
-4
). This suggests a measure such as the standard deviation of
land use could be used as a proxy for habitat heterogeneity on a watershed scale (Sponseller et
al., 2001; Astorga et al., 2014), and may be a metric to consider in subsequent studies. As this
proxy of habitat heterogeneity is negatively correlated with altitude (r = -0.4, p < 10
-4
) this
suggests that a decline in habitat heterogeneity with altitude may then be driving our observed
31
decline in dissimilarity between sampled communities within watersheds with altitude. A similar
relationship may also underpin the positive correlation between land use and dissimilarity as our
measure of habitat heterogeneity is positively correlated with land use (r = 0.6, p < 10
-4
).
𝞯 Diversity decay and evidence of niche differentiation
Our observations, whereby ζ diversity decays via a power-law with respect to ζ order,
support a model of stream community composition being determined more by niche
differentiation than a stochastic assembly process (Scheiner et al., 2011; Hui and McGeoch,
2014; Hui et al., 2018). Both our measures of turnover, β and ζ10 diversity, are more strongly
associated with the standard deviation of land use than altitude within both low and high altitude
watersheds (Tables S3-4). Treating the standard deviation of land use as a measure of landscape
heterogeneity this suggests variability in upstream land use underlies a process of niche
differentiation in metacommunities of BMIs on a watershed scale (Winemiller et al., 2010;
Astorga et al., 2014; Hammill et al., 2018).
Such niche differentiation, whereby the probability of detecting a member of a particular
taxon or functional feeding group within a particular sample varies across a landscape, is also
reflected in the probabilities of detecting members of a particular taxon or functional feeding
group within our set of BMI communities. For example, of the 334 unique taxa identified in our
data only members of three genera were detected in more than 50% of our samples. These were
members of the genera Baetis and Simulium and the subclass Oligochaeta, which are known to
be geographically widespread (Fend et al., 2005). We observed a similar skew in functional
diversity with the most common functional group, collector-gatherers, present in only 37.8% of
samples while the average functional group was only found in 12.5% of samples. This may be
reflective of the relative tolerance of each functional group to the abiotic stressors frequently
32
associated with human dominated land use. Such a case has been observed with members of the
collector-gatherer feeding group (Feld and Hering, 2007). Although we should note that the
BMIs within the family Chironomidae, whose members belong to the collector-gather functional
feeding group, were aggregated to the level of subfamily rather than genus.
Conclusion
Our results indicated that, for multiple communities within watersheds, both the mean
values of αtaxa and α diversity were found to negatively correlate with mean upstream land use
and positively with altitude. This is in line with prior observations of the response of
communities of stream BMIs to natural and anthropogenic variation in the environment (Ward,
1986; Jacobsen et al., 1997; Sponseller et al., 2001; Huryn et al., 2002; Stepenuck et al., 2002;
Allan, 2004; Finn and Poff, 2005; Heino, 2005; Jacobsen, 2008). Such relationships are also
observed in the behavior of the CSCI, itself a measure of taxonomic and functional α diversity
(Rehn et al., 2015; Mazor et al., 2016). What we have also demonstrated is a relationship
between measures of the diversity of our groups of sampled communities, described both by
measures of β and ζ10 diversity, and variations in the environment, both natural and
anthropogenically driven.
As the monitoring of the health of ecological communities transitions from methods
utilizing taxonomic diversity of particular groups of indicator taxa, or other forms of diversity
such as the functional feeding group memberships utilized in this study, to broader community
profiles enabled by metagenomic approaches (Baird and Hajibabaei, 2012; Bohman et al., 2014;
Deiner et al., 2016; Jackson et al., 2016) the need arises for a more generalized framework with
which to construct measures of single and multi-community health which could incorporate
33
information on the composition of full communities (e.g. bacteria, algae, etc.) in addition to
groups of well-defined indicator taxa (Hooper et al., 2003; Takamura et al., 2009; Pignata et al.,
2013). What our results illustrate then is the potential for creating a more generalizable
framework for assessing the health of communities using both established measures of diversity,
such as α and β, as well as comparisons of larger ensembles of communities enabled by the ζ
diversity framework.
Acknowledgements
We would like to acknowledge the support of the Environmental Protection Agency in
funding this analysis. Additionally, we would like to acknowledge Dr. Jed Fuhrman and Dr.
Carly Kenkel for their insights in landscape ecology used in this project.

References
Abell, R., Allan, J. D., & Lehner, B. (2007). Unlocking the potential of protected areas for
freshwaters. Biological Conservation, 134(1), 48-63.

Allan, J. D. (2004). Landscapes and riverscapes: the influence of land use on stream
ecosystems. Annu. Rev. Ecol. Evol. Syst., 35, 257-284.

Anderson, D. M., Glibert, P. M., & Burkholder, J. M. (2002). Harmful algal blooms and
eutrophication: nutrient sources, composition, and consequences. Estuaries, 25(4), 704-726.

34
Astorga, A., Oksanen, J., Luoto, M., Soininen, J., Virtanen, R., & Muotka, T. (2012). Distance
decay of similarity in freshwater communities: do macro-and microorganisms follow the same
rules?. Global Ecology and Biogeography, 21(3), 365-375.

Astorga, A., Death, R., Death, F., Paavola, R., Chakraborty, M., & Muotka, T. (2014). Habitat
heterogeneity drives the geographical distribution of beta diversity: the case of N ew Z ealand
stream invertebrates. Ecology and Evolution, 4(13), 2693-2702.

Baird, D. J., & Hajibabaei, M. (2012). Biomonitoring 2.0: a new paradigm in ecosystem
assessment made possible by next-generation DNA sequencing. Molecular ecology, 21(8), 2039-
2044.

Bohmann, K., Evans, A., Gilbert, M. T. P., Carvalho, G. R., Creer, S., Knapp, M., ... & De
Bruyn, M. (2014). Environmental DNA for wildlife biology and biodiversity monitoring. Trends
in ecology & evolution, 29(6), 358-367.

Brendonck, L., Jocqué, M., Tuytens, K., Timms, B. V., & Vanschoenwinkel, B. (2015).
Hydrological stability drives both local and regional diversity patterns in rock pool
metacommunities. Oikos, 124(6), 741-749.

Carpenter, S. R., Caraco, N. F., Correll, D. L., Howarth, R. W., Sharpley, A. N., & Smith, V. H.
(1998). Nonpoint pollution of surface waters with phosphorus and nitrogen. Ecological
applications, 8(3), 559-568.
35

Chao, A., Jost, L., Chiang, S. C., Jiang, Y. H., & Chazdon, R. L. (2008). A two-stage
probabilistic approach to multiple-community similarity indices. Biometrics, 64(4), 1178-1186.

Chase, J. M., Biro, E. G., Ryberg, W. A., & Smith, K. G. (2009). Predators temper the relative
importance of stochastic processes in the assembly of prey metacommunities. Ecology
letters, 12(11), 1210-1218.

Deiner, K., Fronhofer, E. A., Mächler, E., Walser, J. C., & Altermatt, F. (2016). Environmental
DNA reveals that rivers are conveyer belts of biodiversity information. Nature
communications, 7(1), 1-9.

Dudgeon, D., Arthington, A. H., Gessner, M. O., Kawabata, Z. I., Knowler, D. J., Lévêque, C., ...
& Sullivan, C. A. (2006). Freshwater biodiversity: importance, threats, status and conservation
challenges. Biological reviews, 81(2), 163-182.

Feld, C. K., & Hering, D. (2007). Community structure or function: effects of environmental
stress on benthic macroinvertebrates at different spatial scales. Freshwater Biology, 52(7), 1380-
1399.

Fend, S. V., Carter, J. L., & Kearns, F. R. (2005). Relationships of field habitat measurements,
visual habitat indices, and land cover to benthic macroinvertebrates in urbanized streams of the
Santa Clara Valley, California. In American Fisheries Society Symposium (Vol. 47, pp. 193-212).
36

Finn, D. S., & Leroy Poff, N. (2005). Variability and convergence in benthic communities along
the longitudinal gradients of four physically similar Rocky Mountain streams. Freshwater
Biology, 50(2), 243-261.

Frady, C., Johnson, S., & Li, J. (2007). Stream macroinvertebrate community responses as
legacies of forest harvest at the HJ Andrews Experimental Forest, Oregon. Forest Science, 53(2),
281-293.

Gaston, K. J., & Fuller, R. A. (2008). Commonness, population depletion and conservation
biology. Trends in ecology & evolution, 23(1), 14-19.

Gerth, W. J., & Herlihy, A. T. (2006). Effect of sampling different habitat types in regional
macroinvertebrate bioassessment surveys. Journal of the North American Benthological
Society, 25(2), 501-512.

Grömping, U. (2006). Relative importance for linear regression in R: the package
relaimpo. Journal of statistical software, 17(1), 1-27.

Grönroos, M., Heino, J., Siqueira, T., Landeiro, V. L., Kotanen, J., & Bini, L. M. (2013).
Metacommunity structuring in stream networks: roles of dispersal mode, distance type, and
regional environmental context. Ecology and evolution, 3(13), 4473-4487.

37
Hammill, E., Hawkins, C. P., Greig, H. S., Kratina, P., Shurin, J. B., & Atwood, T. B. (2018).
Landscape heterogeneity strengthens the relationship between β-diversity and ecosystem
function. Ecology, 99(11), 2467-2475.

Hawkins, C. P., Olson, J. R., & Hill, R. A. (2010). The reference condition: predicting
benchmarks for ecological and water-quality assessments. Journal of the North American
Benthological Society, 29(1), 312-343.

Heino, J. (2005). Functional biodiversity of macroinvertebrate assemblages along major
ecological gradients of boreal headwater streams. Freshwater Biology, 50(9), 1578-1587.

Heino, J. (2009). Biodiversity of aquatic insects: spatial gradients and environmental correlates
of assemblage-level measures at large scales. Freshwater reviews, 2(1), 1-29.

Hendrickx, F., MAELFAIT, J. P., Van Wingerden, W., Schweiger, O., Speelmans, M., Aviron,
S., ... & Burel, F. (2007). How landscape structure, land-use intensity and habitat diversity affect
components of total arthropod diversity in agricultural landscapes. Journal of Applied
Ecology, 44(2), 340-351.

Herbst, D. B., & Silldorff, E. L. (2006). Comparison of the performance of different
bioassessment methods: similar evaluations of biotic integrity from separate programs and
procedures. Journal of the north american Benthological society, 25(2), 513-530.

38
Hijmans, R. (2016). geosphere: spherical trigonometryR package version 1.5-5 https://CRAN. R-
project. org/package= geosphere.

Homer, C., Dewitz, J., Fry, J., Coan, M., Hossain, N., Larson, C., ... & Wickham, J. (2007).
Completion of the 2001 national land cover database for the counterminous United
States. Photogrammetric engineering and remote sensing, 73(4), 337.

Hooper, H. L., Sibly, R. M., Hutchinson, T. H., & Maund, S. J. (2003). The influence of larval
density, food availability and habitat longevity on the life history and population growth rate of
the midge Chironomus riparius. Oikos, 102(3), 515-524.

Hui, C., & McGeoch, M. A. (2014). Zeta diversity as a concept and metric that unifies incidence-
based biodiversity patterns. The American Naturalist, 184(5), 684-694.

Hui, C., Vermeulen, W., & Durrheim, G. (2018). Quantifying multiple-site compositional
turnover in an Afrotemperate forest, using zeta diversity. Forest Ecosystems, 5(1), 1-9.

Huryn, A. D., Butz Huryn, V. M., Arbuckle, C. J., & Tsomides, L. (2002). Catchment land-use,
macroinvertebrates and detritus processing in headwater streams: taxonomic richness versus
function. Freshwater Biology, 47(3), 401-415.

Jaccard, P. (1900). Contribution au problème de l'immigration post-glaciare de la flore
alpine. Bull Soc Vaudoise Sci Nat, 36, 87-130.
39

Jackson, M. C., Weyl, O. L. F., Altermatt, F., Durance, I., Friberg, N., Dumbrell, A. J., ... &
Leadley, P. W. (2016). Recommendations for the next generation of global freshwater biological
monitoring tools. In Advances in ecological research (Vol. 55, pp. 615-636). Academic Press.

Jacobsen, D., Schultz, R., & Encalada, A. (1997). Structure and diversity of stream invertebrate
assemblages: the influence of temperature with altitude and latitude. Freshwater Biology, 38(2),
247-261.

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

Latombe, G., McGeoch, M. A., Nipperess, D. A., & Hui, C. (2017). Zetadiv: functions to
compute compositional turnover using zeta diversity. R package version 1.1.1

Latombe, G., Hui, C., & McGeoch, M. A. (2017). Multi-site generalised dissimilarity modelling:
using zeta diversity to differentiate drivers of turnover in rare and widespread species. Methods
in Ecology and Evolution, 8(4), 431-442.

Mazor, R. D., Rehn, A. C., Ode, P. R., Engeln, M., Schiff, K. C., Stein, E. D., ... & Hawkins, C.
P. (2016). Bioassessment in complex environments: designing an index for consistent meaning in
different settings. Freshwater Science, 35(1), 249-271.

40
McGoff, E., Solimini, A. G., Pusch, M. T., Jurca, T., & Sandin, L. (2013). Does lake habitat
alteration and land-use pressure homogenize E uropean littoral macroinvertebrate
communities?. Journal of Applied Ecology, 50(4), 1010-1018.

Munoz, F., Couteron, P., & Ramesh, B. R. (2008). Beta diversity in spatially implicit neutral
models: a new way to assess species migration. The American Naturalist, 172(1), 116-127.

Ode, P. (2003). List of Californian macroinvertebrate taxa and standard taxonomic
effort. Rancho Cordova, CA: California Department of Fish and Game.

Oksanen, J., Blanchet, F. G., Kindt, R., Legendre, P., Minchin, P. R., O’hara, R. B., ... &
Wagner, H. (2013). Community ecology package. R package version, 2-0.

Ourso, R. T., & Frenzel, S. A. (2003). Identification of linear and threshold responses in streams
along a gradient of urbanization in Anchorage, Alaska. Hydrobiologia, 501(1-3), 117-131.

Paerl, H. W., Gardner, W. S., Havens, K. E., Joyner, A. R., McCarthy, M. J., Newell, S. E., ... &
Scott, J. T. (2016). Mitigating cyanobacterial harmful algal blooms in aquatic ecosystems
impacted by climate change and anthropogenic nutrients. Harmful Algae, 54, 213-222.

Paul, M. J., & Meyer, J. L. (2001). Streams in the urban landscape. Annual review of Ecology
and Systematics, 32(1), 333-365.

41
Pearman, P. B. (1997). Correlates of Amphibian Diversity in an Altered Landscape of
Amazonian Ecuador: Correlaciones de la Diversidad de Anfibios en un Paisaje Alterado de la
Amazonía Ecuatoriana. Conservation Biology, 11(5), 1211-1225.

Peck, D. V., Herlihy, A. T., Hill, B. H., Hughes, R. M., Kaufmann, P. R., Klemm, D. J., ... &
Magee, T. (2006). Environmental monitoring and assessment program—surface waters western
pilot study: field operations manual for wadeable streams. EPA 600/R-06/003.

Petrin, Z., Englund, G., & Malmqvist, B. (2008). Contrasting effects of anthropogenic and
natural acidity in streams: a meta-analysis. Proceedings of the Royal Society B: Biological
Sciences, 275(1639), 1143-1148.

Pignata, C., Morin, S., Scharl, A., Traversi, D., Schilirò, T., Degan, R., ... & Coste, M. (2013).
Application of European biomonitoring techniques in China: Are they a useful tool?. Ecological
indicators, 29, 489-500.

Pond, G. J. (2012). Biodiversity loss in Appalachian headwater streams (Kentucky, USA):
Plecoptera and Trichoptera communities. Hydrobiologia, 679(1), 97-117.

R Core Team, R. (2018). R: A language and environment for statistical computing.

42
Rawi, C. S. M., Al-Shami, S. A., Madrus, M. R., & Ahmad, A. H. (2013). Local effects of forest
fragmentation on diversity of aquatic insects in tropical forest streams: implications for
biological conservation. Aquatic Ecology, 47(1), 75-85.

Rehn, A. C., Ode, P. R., & Hawkins, C. P. (2007). Comparisons of targeted-riffle and reach-wide
benthic macroinvertebrate samples: implications for data sharing in stream-condition
assessments. Journal of the North American Benthological Society, 26(2), 332-348.

Rehn, A. C., Mazor, R. D., & Ode, P. R. (2015). The California stream condition index (CSCI): a
new statewide biological scoring tool for assessing the health of freshwater streams. Sacramento,
CA: Surface Water Ambient Monitoring Program (SWAMP).

Reynoldson, T. B., Norris, R. H., Resh, V. H., Day, K. E., & Rosenberg, D. M. (1997). The
reference condition: a comparison of multimetric and multivariate approaches to assess water-
quality impairment using benthic macroinvertebrates. Journal of the North American
Benthological Society, 16(4), 833-852.

Richards, A. B., & Rogers, D. C. (2006). List of freshwater macroinvertebrate taxa from
California and adjacent states including standard taxonomic effort levels. Southwest Association
of Freshwater Invertebrate Taxonomists (SAFIT).

43
Rouquette, J. R., Dallimer, M., Armsworth, P. R., Gaston, K. J., Maltby, L., & Warren, P. H.
(2013). Species turnover and geographic distance in an urban river network. Diversity and
Distributions, 19(11), 1429-1439.

Sala, O. E., Chapin, F. S., Armesto, J. J., Berlow, E., Bloomfield, J., Dirzo, R., ... & Leemans, R.
(2000). Global biodiversity scenarios for the year 2100. science, 287(5459), 1770-1774.

Scheiner, S. M., Chiarucci, A., Fox, G. A., Helmus, M. R., McGlinn, D. J., & Willig, M. R.
(2011). The underpinnings of the relationship of species richness with space and
time. Ecological Monographs, 81(2), 195-213.

Seaber, P. R., Kapinos, F. P., & Knapp, G. L. (1987). Hydrologic unit maps: US Geological
Survey water supply paper 2294. US Geological Survey.

Shade, A., Jones, S. E., & McMahon, K. D. (2008). The influence of habitat heterogeneity on
freshwater bacterial community composition and dynamics. Environmental microbiology, 10(4),
1057-1067.

Shostell, J. M., & Williams, B. S. (2007). Habitat complexity as a determinate of benthic
macroinvertebrate community structure in cypress tree reservoirs. Hydrobiologia, 575(1), 389-
399.

44
Siqueira, T., Bini, L. M., Roque, F. O., Marques Couceiro, S. R., Trivinho-Strixino, S., &
Cottenie, K. (2012). Common and rare species respond to similar niche processes in
macroinvertebrate metacommunities. Ecography, 35(2), 183-192.

Snodgrass, J. W. (1997). Temporal and spatial dynamics of beaver-created patches as influenced
by management practices in a south-eastern North American landscape. Journal of Applied
Ecology, 1043-1056.

Socolar, J. B., Gilroy, J. J., Kunin, W. E., & Edwards, D. P. (2016). How should beta-diversity
inform biodiversity conservation?. Trends in ecology & evolution, 31(1), 67-80.

Sponseller, R. A., Benfield, E. F., & Valett, H. M. (2001). Relationships between land use,
spatial scale and stream macroinvertebrate communities. Freshwater biology, 46(10), 1409-
1424.

Stepenuck, K. F., Crunkilton, R. L., & Wang, L. (2002). Impacts of urban landuse on
macroinvertebrate communities in southeastern Wisconsin streams 1. JAWRA Journal of the
American Water Resources Association, 38(4), 1041-1051.

Strayer, D. L., & Dudgeon, D. (2010). Freshwater biodiversity conservation: recent progress and
future challenges. Journal of the North American Benthological Society, 29(1), 344-358.

45
Takamura, N., Ito, T., Ueno, R., Ohtaka, A., Wakana, I., Nakagawa, M., ... & Nakajima, H.
(2009). Environmental gradients determining the distribution of benthic macroinvertebrates in
Lake Takkobu, Kushiro wetland, northern Japan. Ecological research, 24(2), 371-381.

Tonkin, J. D., Stoll, S., Jähnig, S. C., & Haase, P. (2016). Anthropogenic land-use stress alters
community concordance at the river-riparian interface. Ecological Indicators, 65, 133-141.

Urban, M. C., Skelly, D. K., Burchsted, D., Price, W., & Lowry, S. (2006). Stream communities
across a rural–urban landscape gradient. Diversity and distributions, 12(4), 337-350.

Vander Vorste, R., McElmurray, P., Bell, S., Eliason, K. M., & Brown, B. L. (2017). Does
stream size really explain biodiversity patterns in lotic systems? A call for mechanistic
explanations. Diversity, 9(3), 26.

Wallace, A. M., & Biastoch, R. G. (2016). Detecting changes in the benthic invertebrate
community in response to increasing chloride in streams in Toronto, Canada. Freshwater
Science, 35(1), 353-363.

Ward, J. V. (1986). Altitudinal zonation in a Rocky Moutain stream. Archiv für Hydrobiologie.
Supplementband. Monographische Beiträge, 74(2), 133-199.

Whittaker, R. H. (1972). Evolution and measurement of species diversity. Taxon, 21(2-3), 213-
251.
46

Winemiller, K. O., Flecker, A. S., & Hoeinghaus, D. J. (2010). Patch dynamics and
environmental heterogeneity in lotic ecosystems. Journal of the North American Benthological
Society, 29(1), 84-99.

Zhang, Y., Liu, L., Cheng, L., Cai, Y., Yin, H., Gao, J., & Gao, Y. (2014). Macroinvertebrate
assemblages in streams and rivers of a highly developed region (Lake Taihu Basin,
China). Aquatic biology, 23(1), 15-28.

47
Chapter 2
Using Co-occurrence Network Topology in Assessing Ecological Stress in Benthic
Macroinvertebrate Communities
Abstract
Ecological monitoring of streams has often focused on assessing the biotic integrity of
individual benthic macroinvertebrates (BMIs) communities through local measures of diversity,
such as taxonomic or functional richness. However, as individual BMI communities are
frequently linked by a variety of ecological processes at a regional scale, there is a need to assess
biotic integrity of groups of communities at the scale of watersheds. Using 4619 sampled
communities of streambed BMIs, we investigate this question using co-occurrence networks
generated from groups of communities selected within California watersheds under different
levels of stress due to upstream land use. Building on a number of arguments in theoretical
ecology and network theory, we propose a framework for the assessment of the biotic integrity of
watershed-scale groupings of BMI communities using measures of their co-occurrence network
topology. We found significant correlations between stress, as described by a mean measure of
upstream land use within a watershed, and topological measures of co-occurrence networks such
as network size (r = -0.81, p <10
-4
), connectance (r = 0.31, p <10
-4
), mean co-occurrence strength
(r = 0.25, p <10
-4
), degree heterogeneity (r = -0.10, p <10
-4
), and modularity (r = 0.11, p <10
-4
).
Using these five topological measures we constructed a linear model of biotic integrity, here a
composite of taxonomic and functional diversity known as the California Stream Condition
Index (CSCI), of groups of BMI communities within a watershed. This model can account for
48
66% of among-watershed variation in the mean biotic integrity of communities. These
observations imply a role for co-occurrence networks in assessing the current status of biotic
integrity for BMI communities, as well as their potential use in assessing other ecological
communities.
Introduction
Humanity can be considered a global scale force for ecosystem engineering (Vörösmarty
et al., 2010; Laurance et al., 2014; Guerry et al., 2015). Subsequent to the rise of anthropogenic
stressors on the environment there has been the recognition of the need for ecological monitoring
which can match the scale of human activity (Corona et al., 2011; Bergseth et al., 2015; Foley et
al., 2013; Schmeller et al., 2015; Steenweg et al., 2017).
Among the ecosystems being monitored, streams have been of long-term and ongoing
interest. Human activities are both dependent upon their ecological services (Anderson et al.,
2002; Dudgeon et al., 2006), and frequently a source of their environmental stress (Carpenter et
al., 1998; Paerl et al., 2016). Human-dominated environments, such as farms, tend to cover large
areas. For this reason, there have been efforts to monitor the state of streams across entire
watersheds rather than individual streams (Grönroos et al., 2013; Socolar et al., 2016), especially
in light of the importance of regional versus local measures of habitat quality with stream
communities (Stoll et al., 2016). These biomonitoring efforts have typically focused either on the
presence of certain indicator taxa (Fausch et al., 1990; Vieira et al., 2012), or on comparing the
composition of communities to an ‘undisturbed’ reference (Kerans and Karr, 1994; Masese et al.,
2009; Lakew and Moog, 2015; Mazor et al., 2015; Silva et al., 2017; Vile and Henning, 2018).
49
Historically, bioassessments of stream have tended to be based on data sets composed on
particular communities, such as BMIs (Maxted et al., 2000; Cuffney et al., 2010), organized by
morphological classifications. With the advent of high throughput metagenomic sequencing
there now exists the potential for rapidly constructing a pictures of community composition with
greater breadth (Stein et al., 2014; Elbrecht et al., 2017), taxonomic resolution, and reliability
(Sweeney et al., 2011; Baird and Hajibabaei, 2012). There is an opportunity then to create a
bioassessment framework for BMI communities in which one could readily incorporate
community composition data constructed from metagenomic methods (Goodwin et al., 2017;
Hering et al., 2018).
Here we propose the use of co-occurrence networks to the task of ecological monitoring.
These networks represent the likelihoods, represented by edges, of various unique categories of
taxa, represented by nodes, co-occurring in a landscape defined by the spatial extent of the
communities sampled and studied. These networks can be constructed from basic ecological
data, such as the presence or absences of a set of taxonomic groups across a set of sites (Gotelli,
2000; Arita, 2016; Morueta-Holme et al., 2016). Co-occurrence networks have been investigated
as a means of inferring ecological patterns, particularly when direct measurement of ecological
interactions proves infeasible. For example, the clustering of microbial species into distinct
modules within co-occurrence networks have been used to infer physiochemical niches for
various prokaryotic groups (Ruan et al., 2006; Fuhrman and Steele, 2008; Steele et al., 2011;
Widder et al., 2014; Larsen and Omerod, 2014; Mandakovic et al., 2018). In studies of larger
organisms, topological measures of these networks have also been used to illustrate a loss in both
diversity and the number of significant co-occurrences between reptiles in response to habitat
degradation (Kay et al., 2018). The previous diversity of scenarios where co-occurrence network
50
topology has been used in ecological analysis then implies that it could also be used to develop a
framework for the assessment of the biotic integrity of streams across an entire catchment area
(Moyle and Randall, 1998; Smith and Lamp, 2008; Ahn et al., 2017).
We hypothesized a number of relationships between ecological stress and five measures
of co-occurrence network topology (Table 1). We chose these measures based on prior analyses
of co-occurrence networks and their relationships with environmental stress. For example, spatial
aggregation of species across an environmentally heterogeneous landscape due to variations in
their ecological attributes (Bellisario et al., 2010; Borthagaray et al., 2014). These hypotheses
were then tested using co-occurrence networks generated from presence/absence data for benthic
macroinvertebrates gathered in streams across the state of California. We used upstream land use
as our measure of stress as it has been found to be a broad measure of anthropogenic stress in
stream communities (Novotny et al., 2005; Vaander Laan et al., 2013).

51
Table 1: Topological measures of co-occurrence networks, their ecological relevance, and predicted relationship with
an increase in stress due to upstream land use: network size, connectance, mean co-occurrence strength,
modularity, and degree heterogeneity.
Topological measure Ecological relevance Hypothesized
relationship with
stress
Network size The number of unique types of
taxa across a set of communities.
(-)
Connectance The fraction of significant co-
occurrences realized compared to
theoretical maximum for a
network.
(+)
Mean co-occurrence strength Correlation strength between
unique types of taxa.
(+)
Modularity How strongly patterns of co-
occurrence are partitioned into
sub-communities.
(+)
Degree heterogeneity How skewed the distribution of the
number of co-occurrences per
unique type of taxa are in a
community.
(-)

52
Ecological stress and co-occurrence network topology
Prior observations of BMI communities under stress have shown two trends. First, a
decline in taxonomic richness (Lenat and Crawford, 1994; Stepenuck et al., 2004; Voß and
Schäfer, 2017). Second, a predominance of members of generalist groups with broad ecological
niches (Büchi and Vuilleumier, 2014; Ducatez et al., 2014; Mykrä and Heino, 2017). Starting
with these trends we then hypothesized relationships between five topological measures and
ecological stress.
Network size
The sizes of our networks were determined by the number of unique BMI genera present
within a given set of sampling sites. Ecological stress in BMI communities has been found to be
associated with a decline in local taxonomic richness (Paul and Meyer, 2001; Stepenuck et al.,
2002; Ourso and Frenzel, 2003; Wallace and Biastoch, 2016). This is especially the case where
the stress is due to an increase in upstream land use (Sponseller et al., 2001; Stepenuck et al.,
2002; Allan, 2004). Given then the correspondence between the taxonomic richness present in a
group of sites, and the number of nodes in any resulting co-occurrence network, we expect
network size to be negatively correlated with stress (H1).
Connectance
As the number of edges in a variety of ecological networks may be sensitive to the
number of unique taxonomic groups (Goldwasser and Roughgarden, 1997; Nielsen and
Bascompte, 2007; Dormann et al, 2009) we then also calculated the connectance (Bell et al.,
2010). Generalists are expected to have a greater likelihood of co-occurring with a wider variety
of organisms (Fridley et al., 2007), and BMIs communities in degraded environments tend to
53
contain relatively more groups classified as generalist (Suga and Tanaka, 2013; von der Ohe and
Goedkoop, 2013). We then expect stress to be positively correlated with the fraction of realized
versus the potential number of edges (significant co-occurrences), that is the connectance, of a
resulting co-occurrence network (H2). With an expected increase in connectance associated with
stress, as well as a decline in number of nodes (number of unique BMIs) (Shaver et al., 1994;
Blann et al., 2009), we also expect a decline in the number of edges.
Mean co-occurrence strength
To make additional inferences on shifts in community co-occurrence patterns in relation
to environmental stress we then determined the mean strength of the co-occurrences found
within each network (Araújo and Rozenfeld, 2014). For this value we used the mean value of all
of the significant correlations, as described by standardized effect-size scores (Morueta-Holme et
al., 2016), within a network. For any two unique categories of organisms found in a group of
communities the standardized effect size score represents the conditional probability, as
compared to a null model, of observing one organism given the presence of the other. The mean
strength of correlations defining significant co-occurrences in a network has been observed to
decline with the number of edges (Cazalles et al., 2016). We then expect the number of edges in
a co-occurrence network to decline with stress, as described by a mean measure of upstream land
use within a watershed, along with a positive correlation between stress and the mean co-
occurrence strength (H3).
Modularity
Prior evidence suggests shifts in communities in response to environmental changes can
be better illustrated not just from the number or strength of co-occurrences, but from their
54
structural arrangement (Fortuna et al., 2010; Thébault and Fontaine, 2010; Tylianakis et al.,
2010). To measure these structural changes in our co-occurrence networks we used the
topological measures of modularity, defined here as the proportion of edges which occur within
sub-networks less the expected proportion of such edges (Clauset et al., 2004). With highly
modular networks This would be expected to lead to a co-occurrence network composed of
sparsely interconnected sub-networks.
Prior observations of stressed watersheds have shown both a decline in local diversity, as
well as a rise landscape diversity as a result of declining taxonomic similarities between
individual stream communities (Hawkins et al., 2015; Simons et al., 2019). Given these prior
observations, with regards to changes in patterns of diversity across watersheds in relation to
stress, we expect the taxonomic ‘space’ for co-occurrences to shrink with a rise in stress due to
land use (Figure 1), and with it a trend towards the fracturing of assembled co-occurrence
networks into weakly connected subnetworks. Similar relationships, between the modularity of
co-occurrence networks and ecological stress, have also been observed in various ecological
communities (Hu et al., 2017; Kay et al., 2018). To assess these trends we use modularity, the
degree to which networks are organized into clusters of weakly interconnected subnetworks
(Clauset et al., 2004; Barberán et al., 2011). With stress expected to drive greater dissimilarity
between communities we then expect a positive relationship between the modularity of co-
55
occurrence networks and the levels of stress experienced by the communities from which they
are constructed (H4).

Degree heterogeneity
To further investigate changes in the arrangement of co-occurrences we determined the
degree heterogeneity of each network, a measure of how skewed the distribution of edges per
node in a network is towards the most connected nodes (Yan et al., 2017). The distribution of
edges per node in ecological networks can be indicative of the structure of ecological
communities, such the likelihood of co-occurrence between generalist and specialist species
(Dormann et al., 2009; Williams, 2011).
Prior observations of co-occurrence networks assembled from communities at increasing
levels of anthropogenic disturbance have shown a trend towards the preferential loss of taxa of
low degree (Fournier et al., 2016; Tulloch et al., 2016). Additional trends regarding ecological
Figure 1: An example of a stressor reducing both the taxonomic richness of 3 communities, from an initial
state (𝞪 1, 𝞪 2, 𝞪 3) to a degraded state (𝞪 1’, 𝞪 2’, 𝞪 3’), as well as the number of unique categories of taxa held
in common between communities.
56
networks have also displayed trends towards the loss of highly keystone taxa due to
environmental stresses (Araújo and Rozenfeld, 2014; Morriën et al., 2017). Given both of these
trends, the loss of taxa of both high and low degree, we then expect co-occurrence networks
assembled by communities under stress have a narrow degree distribution and thus a low degree
heterogeneity.
Using prior arguments regarding connectance (H2) we also make an additional argument
regarding our expected trends in the degree with respect to land use. With ecological networks
connectance has been found to be negatively correlated with the skewness of their degree
distributions (Poisot and Gravel, 2014). With degree heterogeneity being a measure of skew for
the degree distribution of a network this then implies that the stress experienced by the
communities used to construct co-occurrence networks will be negatively correlated with their
degree heterogeneity (H5).
Materials and methods
Sample scope
The initial scope of data covered in this analysis consists of 4984 stream samples from
2997 unique geographic locations across the state of California, constituting a 23-year period
(1994-2016) (Mazor et al., 2015). Every sample contains the following data: benthic
macroinvertebrates (BMIs) enumerated and sorted to a standardized level (generally a genus-
level identification except chironomids which were identified to subfamily, Richards and Rogers,
2011), sample site altitude in meters, U.S. Geological Survey Hydrologic Unit Code 8 level
watershed (Seaber et al., 1987), and the percent developed land use (agricultural, urban, and
managed landscapes) within a 5 km clipped buffer of the watershed upstream of the sampling
57
site, and a bioassessment index score (CSCI) based on a composite of taxonomic and functional
diversity within BMI assemblages (Mazor et al., 2015).
Sample acquisitions and classifications
Approximately 55% of the BMI communities were sampled through a reach-wide
protocol of (Peck et al., 2006), with the remainder collected using a targeted riffle protocol, both
of which produce comparable data (Gerth and Herlihy, 2006; Herbst and Silldorff, 2006; Rehn et
al., 2007). Taxa were identified to one of 334 genera, with each genus then assigned to one of 8
functional feeding groups using CAMLnet (Ode, 2003). Of these 8 functional feeding groups we
could conclusively assign 5 of them to either generalist or specialist categories (Mihuc, 1997;
Barbour et al., 2006; Feld and Hering, 2007; Rawer-Jost et al., 2011; De Castro et al., 2016).
Using this information, we produced a measure of the number of generalist and specialist genera
per sample site.
Calculating the CSCI
Our measure of community biotic integrity at a given stream sample site was done using
the CSCI. This index compares observed taxa and metrics to values expected under undisturbed
reference conditions based on site-specific landscape-scale environmental variables, such as
watershed area, geology, and climate (Mazor et al., 2015). This index comprises two sets of
measurements using a standardized taxonomy for BMI communities (Richards and Rogers,
2011): the first being a ratio of observed-to-expected taxa (O/E), and the second a predictive
multi-metric index (pMMI) made of 6 metrics related to ecological structure and function of the
BMI assemblage describing the composition of community within a site. The CSCI, and its
components, were designed to have minimal influence from major natural gradients. This in turn
58
has allowed for it to be used as a measure of biological conditions with a consistent meaning in
different environmental settings (Reynoldson et al. 1997; Hawkins et al. 2010).
Land use
The type and geographic extent of land use in the upstream vicinity of each sampling site
data is derived from the National Land Cover Data set (NLCD) (Homer et al., 2007), with
developed land cover measured by the total percent of land cover in a designated area dedicated
to agriculture, urbanization, or otherwise managed vegetative landscapes such as golf courses.
The designated area for calculating percent developed land cover at each site is defined using a 5
km watershed-clipped buffer upstream of a stream sampling site using ArcGIS tools (version
10.3, Environmental Systems Research Institute, Redlands, California) (Mazor et al., 2015). The
values for land use were calculated from NLCD measurements acquired in the year 2000, though
it should be noted that the sample sites in our study were located in areas where the percent
developed land use was not significantly correlated with time over the duration of this study (r =
-0.02, p = 0.27).
Sample group selection
We first filtered our initial data by selecting watersheds with 15 or more unique samples.
This filtering reduced our overall data set from 4984 to 4619 unique samples in 2694 unique
geographic locations across 67 watersheds, while containing sample groups with sufficient data
density for co-occurrence network construction. From these remaining samples we then divided
both upstream land use and sample site altitude into quintiles. Samples for network generation
were then selected by randomly subsampling 10 samples within each watershed within quintiles
of upstream land use and sample site altitude. For each group of 10 samples we calculated the
59
mean sample site altitude (altitude), the mean percent developed upstream land use (land use),
and for mean geographic separation distance in meters between samples (distance) we used the
distm function within the R package geosphere (Hijmans et al., 2012). To obtain a measure of
environmental heterogeneity within each sample group we also calculated the standard
deviations on altitude and land use.
Network construction
Co-occurrence networks were then constructed using the R package netassoc (Morueta-
Holme et al., 2016), with presence/absence site by BMI genera as input. We chose to convert our
abundance data to presence/absence as the most conservative approach with representing our
assembled database of BMI communities. Observed co-occurrences were compared against 100
randomized null communities with the same taxonomic richness as the observed community.
The resulting edges were filtered so only correlations representing co-occurrences, as calculated
by standardized effect-size scores, with a significance and false discovery rate less that 10
'(

were kept. This process was repeated 100 times, with a set of 8208 co-occurrence networks kept
for analysis.
Topological measures
Topological measures of our networks, such as size and connectance, were calculated
using the packages igraph (v.1.2.2) (Csardi and Nepusz, 2014) and network (v.1.13.0.1) (Butts,
2015) in R (v.3.5.1). The mean co-occurrence strength values were calculated, using the R
package netassoc (Morueta-Holme et al., 2016), as the network mean of the significant
standardized effect-size scores. Modularity, defined as the proportion of edges which occur
within sub-networks less the expected proportion of such edges, was calculated using the
60
modularity function within the igraph package (Clauset et al., 2004). Degree heterogeneity was
calculated as ζ=
+,
-
.
+,.
-
, where k represents the mean number of edges per node in a network
(Yan et al., 2017).
Modeled biotic integrity index
Using the lm function in the stats R package (v3.5.1 R Core Team, 2018) we constructed
a best-fit linear model to predict the mean CSCI score of a set of samples, our measure of biotic
integrity, given the topological measures of their co-occurrence networks. We then applied a
backwards elimination method in order to select topological measures which make a significant
contribution to our model (Pearman, 1997; Snodgrass, 1997). In comparing the AIC scores of
linear models after the removal of each topological measure we found all five were significant.
We calculated coefficients for our linear models using a 10-fold cross-validation, with 100
repeats, within the ‘train’ function within the R package ‘caret’ (Kuhn, 2008). To determine the
relative importance of each topological measure in our linear models, and to adjust for any
collinearity between measures as a result, the function calc.relimp was used within the relaimpo
R package (Grömping, 2006). The relative importance of land use, altitude, and distance in
describing variations in both the mean CSCI score, as well as our modeled CSCI scores, were
also done using the calc.relimp function.
Results
In analyzing 8208 co-occurrence networks, generated from communities collected from
within-watershed groups with similar values for sample site land use and altitude, we found
general support (Table 2) for our hypotheses (Table 1).
61

Table 2: The relative importance of the topological measures used in our modeled stream health indices (p < 10
-4
).
Topological
measure
F(1,8208)
(Model 1)
Relative
importance (%)
(Model 1)
F(1,8208)
(Model 2)
Relative
importance (%)
(Model 2)
Network size 1.4 x 10
4
44.8 NA NA
Connectance 704.3 10.0 2086 19.3
Modularity 151.4 1.3 39.9 1.8
Mean co-
occurrence
strength
868.6 7.5 2462 13.4
Degree
heterogeneity
91.6 2.7 387.1 3.3

Trends in co-occurrence network topology
The size of our co-occurrence networks declined significantly with a rise in land use (r =
-0.81, p < 10
-4
). This reflects a general decline in both the number of genera found in an
individual sampling site (r = -0.52, p < 10
-4
), as well as the number of functional feeding groups
(r = -0.44, p < 10
-4
), in relation to land use.
While network size was found to have a strong negative correlation with land use, along
with the mean number of edges per node (r = -0.56, p < 10
-4
), we still found that connectance
tended to be larger in co-occurrence networks constructed from groups of stressed communities
with a rise in land use (r = 0.31, p < 10
-4
). This positive association between stress and
62
connectance appears to reflect a greater relative decline in the number of nodes relative to land
use (r = -0.81, p < 10
-4
) than with the number of edges (r = -0.70, p < 10
-4
).
These trends, a rise in connectance despite a decline in network size, may also reflect our
observations regarding the relative abundance of unique genera classified by membership of a
generalist or specialist functional feeding groups to land use. We found the proportion of genera
from specialist functional feeding groups (e.g. shredders and scrapers) tended decline with land
use, while those of generalist functional feeding groups (e.g. gatherers, filterers, and omnivores)
tended to increase with land use (Table 3).

Table 3: Coefficents of sample site altitude and land use in linear models describing linear models of the percent of
genera of BMIs per sample site per functional feeding group (All p < 10
-4
unless otherwise noted).
Generalist functional feeding groups Specialist functional feeding
groups
Gatherers Filterers Omnivores Scrapers Shredders
Coefficient
(land use)
1.6 x 10
-3
1.5 x 10
-4
3.5 x 10
-4
-9.0 x 10
-5
(p < 10
-2
)
-1.5 x 10
-4

Coefficient
(altitude)
8.0 x 10
-6
(p < 10
-2
)
-1.3 x 10
-5
-5.8 x 10
-6
-9.0 x 10
-5
1.1 x 10
-5

In addition to a rise in connectance, networks assembled from communities with higher
land use were, on average, found to have stronger co-occurrences (r = 0.25, p < 10
-4
). We also
found evidence of a negative relationship between both mean co-occurrence strength and the
number of co-occurrences (r = -0.32, p < 10
-4
), as well as connectance (r = -0.24, p < 10
-4
). This
63
potentially indicates a preferential loss of weak co-occurrences in networks assembled from
communities under high levels of land use.
Weaker trends were observed with regards to variables, modularity and degree
heterogeneity, which describe structural arrangements of co-occurrences. The mean modularity
of our networks (0.35) was found to be both greater than the common modularity threshold of
0.3 (Newman and Girvan, 2004), as well as that of our randomized null co-occurrence networks
(0.22). Using a Wilcoxon signed-rank test we found additional evidence for significant non-
random structuring in our networks as their modularity values were significantly larger than their
randomized null counterparts (p < 10
-4
). However, despite evidence of significant network
modularity there was only a relatively weak positive correlation with its value and land use (r =
0.11, p < 10
-4
).
Across our watersheds we observe a trend where land use is associated with a slight
decline in degree heterogeneity (r = -0.10, p < 10
-4
). Although we did find support for our
hypothesis that a decline in degree heterogeneity would be driven, at least in part, by a rise in
connectance (r = -0.60, p < 10
-4
). Similar to our results with modularity we found evidence,
using a Wilcoxon signed-rank test, for significantly higher values for degree heterogeneity in our
co-occurrence networks than their randomized null counterparts (p < 10
-4
). The higher mean
degree heterogeneity of our co-occurrence networks (1.82), as compared to that of the null
networks (1.14), indicates our networks are skewed more towards a relatively small number of
highly connected nodes than what would be expected by chance.
64
Linear models of watershed biotic integrity using co-occurrence network
topology
Using five measures of co-occurrence network topology, network size (N), connectance
(C), mean co-occurrence strength (S), modularity (M), and degree heterogeneity (𝞯), a linear
model was constructed to best predict the mean value of the CSCI score for a set of samples
(Table 2). The relationship between these topological measures and our first modeled mean
CSCI score per sample group are as follows:

MeanCSCI = 0.3 + 4.6 x 10
-3
x N - 1.2 x C – 1.8 x 10
-2
x S + 0.3 x M + 0.1 x ζ

This modeled index of watershed biotic integrity was found to be strongly correlated with the
observed variation in the mean value of the CSCI score for a set of samples (Figure 2). After
performing a 10-fold cross-
validation this model could still
account for approximately 66%
of the observed variation in the
mean CSCI score. This
modeled biotic integrity index
was also found to vary in
accordance with altitude, land
use, and distance for a set of samples in a similar fashion as the mean CSCI score, although this
first modeled index was less sensitive to altitude and the standard deviation on land use than the
mean CSCI score (Table 4). We also observed that most of the variations observed in both the
Figure 2: A comparison of the first modeled CSCI and mean CSCI
colored by land use (r = 0.81, p < 10
-4
).
65
mean and our first modeled CSCI scores were driven by land use and the standard deviation on
land use (Table 4).

Both network size and the
CSCI, our measure of biotic
integrity, represent measures based
on the taxonomic diversity of
sampled communities. To focus on
the potential role of the
characteristics and configuration of
our co-occurrences, rather than
measures of local diversity alone, we then generated a second model of the mean CSCI with
network size removed from our list of topological measures (Table 2). After performing a 10-
fold cross-validation we found this second linear model can account for 38% of the observed
variation in the mean CSCI score per group of samples, and it exhibits a similar trend compared
to the mean CSCI as with our first model (Figure 3).

Figure 3: A comparison of the second modeled CSCI and mean CSCI
colored by land use (r = 0.61, p < 10
-4
).
66
Table 4: The relative importance of altitude, standard deviation on altitude, land use, standard deviation of land use,
and distance in describing our linear models of the mean value of the CSCI and modeled index per HUC 8 watershed
(All p < 10
-4
unless otherwise noted).
CSCI Modeled index 1 Modeled index 2
Proportion of
variation due to
altitude F(1,8208)
690.5 9.0 (p < 10
-2
) 189.2
Proportion of
variation due to the
standard deviation on
altitude F(1,8208)
26.6 51.8 302.9
Proportion of
variation due to land
use F(1,8208)
2.2 x 10
4
1.6 x 10
4
3082
Proportion of
variation due to the
standard deviation on
land use F(1,8208)
776.2 212.5 217.7
Proportion of
variation due to
distance F(1,8208)
28.4 168.4 2.8
Relative importance
of altitude (%)
11.6 5.2 7.3
Relative importance
of the standard
7.4 4.1 2.7
67
deviation on altitude
(%)
Relative importance
of land use (%)
33.6 35.2 11.9
Relative importance
of the standard
deviation on land use
(%)
20.7 18.4 9.3
Relative importance
of distance (%)
1.0 3.3 0.4
Proportion of
variance explained
by model (%)
74.2 66.3 31.6
Discussion
We found changes in patterns of co-occurrence between genera of BMIs can play a role
in describing effects of land use on regional measures of biotic integrity. This is reflected in
evidence supporting our hypotheses regarding the relationships between land use and the
connectance of our co-occurrence networks (H2), as well as the mean strength of their co-
occurrences (H3). Evidence supporting our hypothesis regarding a negative correlation between
network size and land use (H1), reflecting a well-established link between environmental stress
and both the loss of biodiversity as well as measures of biotic integrity (Garie and McIntosh,
1986; Freeman and Schorr, 2004; Jun et al., 2011; Jun et al., 2016). We find the importance of
network size, along with network connectance and co-occurrence strength, reinforces prior
68
observations on the importance of both regional and local measures of environmental quality in
BMI communities (Stoll et al., 2016).
Trends relating the arrangement of co-occurrences within of our networks, as described
by our hypotheses regarding modularity (H4) and degree heterogeneity (H5), were less clear.
This may reflect limitations in our use of co-occurrence rather than co-abundance networks.
Although prior evidence from assessments of biotic integrity for stream communities of BMIs
have shown a strong correlation between results generated using community data sets described
through abundance or presence/absence (Beentjes et al., 2018). The more fundamental issue may
stem from differences between networks assembled from co-occurrences rather than interactions.
Analyses of co-occurrence networks have been used to identify candidate keystone taxa
(Berry and Widder, 2014), potential species interactions (Veech, 2013), and the simplification of
communities under ecological stress (Araújo et al., 2011). Though inferring co-occurrences,
rather than verifying interactions, is a far more tractable problem in complex ecological systems
we must acknowledge that co-occurrences do not necessarily imply interactions. An underlying
caveat with analyses involving co-occurrences is that various types of ecological interactions,
such as mutualism or similar environmental requirements, may produce similar patterns of co-
occurrence (Ovaskainen et al., 2010). In the context of our study, we observed trends between
land use and both network connectance and co-occurrence strength with our BMI communities
which may reflect changes in patterns of interaction between members of generalist genera.
Though our co-occurrence networks may also be incorporating information beyond potential
interactions between species, such as the tendency of organisms with similar ecological niches to
form co-occurrences, or for dispersal limitation to tend to limit them (Morueta-Holme et al.,
2016). While co-occurrence networks, such as the ones we’ve constructed, may only describe
69
potential interactions, they can still provide useful indications of changes in ecological systems
(Freilich et al., 2018).
Even with these limitations we found a simple linear model composed of topological
measures of co-occurrence networks could describe a significant portion of the observed
variation in the biotic integrity of our BMI communities (Table 2). Analysis of these models also
suggest the topology our networks reflect more than changes in local biodiversity. While
network size contributes a sizeable portion of the observed variation in biotic integrity, its
removal still leaves more than half of the remaining explanatory power of our linear model of the
mean CSCI score (Table 2). This suggests that we are not simply observing a decline in local
diversity in response to stress but a change in landscape diversity as well.
Variations in our models appeared to be driven more by both land use and the standard
deviation on land use than either altitude or geographic separation distance (Table 4). This
reflects our prior observations of this system (Simons et al., 2019), whereby both the mean and
standard deviation of land use, are strongly correlated with degree of taxonomic dissimilarity
between communities. These results are also in agreement with studies of other BMI
communities where measures of environmental heterogeneity, such as variations in upstream
land use between sample sites (Sponseller et al. 2001, Astorga et al. 2014), appear to drive
significant changes in patterns of co-occurrence (Shostell and Williams 2007; Heino, 2013a;
Larsen and Ormerod, 2014; Zhang et al. 2014).
These trends may reflect co-occurrence patterns unique to stream communities of BMIs.
However, the framework we’ve used to test these hypotheses is not dependent upon the
particular identities of taxa present in communities and may have the potential to be applied to
other systems. Assessments of the biotic integrity of freshwater ecosystems have been carried out
70
at a variety of spatial scales (Booth et al., 2004; King et al., 2011; Pratt and Chang, 2012),
regions (Waite et al., 2010; Weigel and Dimick, 2011; Jun et al., 2012), and biological
communities (Zalack et al., 2010; Ferreira et al., 2011; Fetscher et al., 2014). With co-occurrence
networks we then assert the potential for the development of a more flexible framework for the
monitoring of freshwater ecosystems, and find this direction warrants further research.
Synthesis and future directions
It is increasingly becoming feasible to characterize entire ecological communities, from
prokaryotes through metazoa, through metagenomic approaches (Baird and Hajibabaei, 2012;
Bohman et al., 2014; Deiner et al., 2016; Jackson et al., 2016). With the ability to generate such
broad and deep pictures of multiple communities there is a commensurate need to create a
framework which could evaluate in general patterns in ecological systems order to evaluate the
biotic integrity of ecosystems. Using these stream communities as an example our study suggests
significant relationships exist between ecological stress and the structure of co-occurrence
networks. We found our co-occurrence networks reflected changes in the structure of within-
watershed communities in accordance with both the mean and standard deviation on land use and
altitude, potentially reflecting segregation in BMI communities with a rise in measures of
environmental heterogeneity. Using only a small number of topological measures we were also
able to construct a simple linear model with a close correspondence to a well-accepted index of
biotic integrity. These networks are based on patterns derived from the co-occurrences of unique
organisms, rather than their identities. For such reasons we believe co-occurrence networks may
have the potential to describe the biotic integrity of a diverse array of ecological communities.
71
Acknowledgements
We would like to acknowledge Dr. Eric D. Stein and Dr. Sergey Nuzhdin for their
insights in ecology used in this project. This research was supported with an Environmental
Protection Agency (grant number 53-4814-5430). No conflict of interest exists in this
manuscript.
References
Ahn, S. R., & Kim, S. J. (2017). Assessment of integrated watershed health based on the natural
environment, hydrology, water quality, and aquatic ecology. Hydrology and Earth System
Sciences, 21(11), 5583.

Allan, J. D. (2004). Landscapes and riverscapes: the influence of land use on stream
ecosystems. Annu. Rev. Ecol. Evol. Syst., 35, 257-284.

Anderson, D. M., Glibert, P. M., & Burkholder, J. M. (2002). Harmful algal blooms and
eutrophication: nutrient sources, composition, and consequences. Estuaries, 25(4), 704-726.

Araújo, M. B., & Rozenfeld, A. (2014). The geographic scaling of biotic
interactions. Ecography, 37(5), 406-415.

Araújo, M. B., Rozenfeld, A., Rahbek, C., & Marquet, P. A. (2011). Using species co-occurrence
networks to assess the impacts of climate change. Ecography, 34(6), 897-908.

72
Arita, H. T. (2016). Species co-occurrence analysis: pairwise versus matrix-level
approaches. Global Ecology and Biogeography, 25(11), 1397-1400.

Astorga, A., Oksanen, J., Luoto, M., Soininen, J., Virtanen, R., & Muotka, T. (2012). Distance
decay of similarity in freshwater communities: do macro-and microorganisms follow the same
rules?. Global Ecology and Biogeography, 21(3), 365-375.

Baird, D. J., & Hajibabaei, M. (2012). Biomonitoring 2.0: a new paradigm in ecosystem
assessment made possible by next-generation DNA sequencing. Molecular ecology, 21(8), 2039-
2044.

Barberán, A., Bates, S. T., Casamayor, E. O., & Fierer, N. (2012). Using network analysis to
explore co-occurrence patterns in soil microbial communities. The ISME journal, 6(2), 343-351.

Barbour, M. T., Gerritsen, J., Griffith, 3., Frydenborg, R., McCarron, E., White, J. S., & Bastian,
M. L. (1996). A framework for biological criteria for Florida streams using benthic
macroinvertebrates. Journal of the North American Benthological Society, 15(2), 185-211.

Beentjes, K. K., Speksnijder, A. G., Schilthuizen, M., Schaub, B. E., & van der Hoorn, B. B.
(2018). The influence of macroinvertebrate abundance on the assessment of freshwater quality in
The Netherlands. Metabarcoding and Metagenomics, 2, e26744.

73
Bell, J. R., Andrew King, R., Bohan, D. A., & Symondson, W. O. (2010). Spatial co-occurrence
networks predict the feeding histories of polyphagous arthropod predators at field
scales. Ecography, 33(1), 64-72.

Bellisario, B., Cerfolli, F., & Nascetti, G. (2010). Spatial network structure and robustness of
detritus-based communities in a patchy environment. Ecological research, 25(4), 813-821.

Bergseth, B. J., Russ, G. R., & Cinner, J. E. (2015). Measuring and monitoring compliance in no-
take marine reserves. Fish and fisheries, 16(2), 240-258.

Berry, D., & Widder, S. (2014). Deciphering microbial interactions and detecting keystone
species with co-occurrence networks. Frontiers in microbiology, 5, 219.

Blann, K. L., Anderson, J. L., Sands, G. R., & Vondracek, B. (2009). Effects of agricultural
drainage on aquatic ecosystems: a review. Critical reviews in environmental science and
technology, 39(11), 909-1001.

Bohmann, K., Evans, A., Gilbert, M. T. P., Carvalho, G. R., Creer, S., Knapp, M., ... & De
Bruyn, M. (2014). Environmental DNA for wildlife biology and biodiversity monitoring. Trends
in ecology & evolution, 29(6), 358-367.

74
Booth, D. B., Karr, J. R., Schauman, S., Konrad, C. P., Morley, S. A., Larson, M. G., & Burges,
S. J. (2004). Reviving urban streams: Land use, hydrology, biology, and human behavior
1. JAWRA Journal of the American Water Resources Association, 40(5), 1351-1364.

Borthagaray, A. I., Arim, M., & Marquet, P. A. (2014). Inferring species roles in
metacommunity structure from species co-occurrence networks. Proceedings of the Royal
Society B: Biological Sciences, 281(1792), 20141425.

Büchi, L., & Vuilleumier, S. (2014). Coexistence of specialist and generalist species is shaped by
dispersal and environmental factors. The American Naturalist, 183(5), 612-624.

Butts, C. T. (2008). network: a Package for Managing Relational Data in R. Journal of statistical
software, 24(2), 1-36.

Carpenter, S. R., Caraco, N. F., Correll, D. L., Howarth, R. W., Sharpley, A. N., & Smith, V. H.
(1998). Nonpoint pollution of surface waters with phosphorus and nitrogen. Ecological
applications, 8(3), 559-568.

Cazelles, K., Araújo, M. B., Mouquet, N., & Gravel, D. (2016). A theory for species co-
occurrence in interaction networks. Theoretical Ecology, 9(1), 39-48.

Clauset, A., Newman, M. E., & Moore, C. (2004). Finding community structure in very large
networks. Physical review E, 70(6), 066111.
75

Corona, P., Chirici, G., McRoberts, R. E., Winter, S., & Barbati, A. (2011). Contribution of
large-scale forest inventories to biodiversity assessment and monitoring. Forest ecology and
management, 262(11), 2061-2069.

Csardi, G., & Nepusz, T. (2006). The igraph software package for complex network
research. InterJournal, complex systems, 1695(5), 1-9.

Cuffney, T. F., Brightbill, R. A., May, J. T., & Waite, I. R. (2010). Responses of benthic
macroinvertebrates to environmental changes associated with urbanization in nine metropolitan
areas. Ecological Applications, 20(5), 1384-1401.

Parreira de Castro, D. M., Reis de Carvalho, D., Pompeu, P. D. S., Moreira, M. Z., Nardoto, G.
B., & Callisto, M. (2016). Land use influences niche size and the assimilation of resources by
benthic macroinvertebrates in tropical headwater streams. PLoS One, 11(3), e0150527.

Deiner, K., Fronhofer, E. A., Mächler, E., Walser, J. C., & Altermatt, F. (2016). Environmental
DNA reveals that rivers are conveyer belts of biodiversity information. Nature
communications, 7(1), 1-9.

Dormann, C. F., Fründ, J., Blüthgen, N., & Gruber, B. (2009). Indices, graphs and null models:
analyzing bipartite ecological networks. The Open Ecology Journal, 2(1).

76
Ducatez, S., Tingley, R., & Shine, R. (2014). Using species co-occurrence patterns to quantify
relative habitat breadth in terrestrial vertebrates. Ecosphere, 5(12), 1-12.

Dudgeon, D., Arthington, A. H., Gessner, M. O., Kawabata, Z. I., Knowler, D. J., Lévêque, C., ...
& Sullivan, C. A. (2006). Freshwater biodiversity: importance, threats, status and conservation
challenges. Biological reviews, 81(2), 163-182.

Elbrecht, V., Vamos, E. E., Meissner, K., Aroviita, J., & Leese, F. (2017). Assessing strengths
and weaknesses of DNA metabarcoding-based macroinvertebrate identification for routine
stream monitoring. Methods in Ecology and Evolution, 8(10), 1265-1275.

Fausch, K. D., Lyons, J. O. H. N., Karr, J. R., & Angermeier, P. L. (1990, December). Fish
communities as indicators of environmental degradation. In American fisheries society
symposium (Vol. 8, No. 1, pp. 123-144).

Feld, C. K., & Hering, D. (2007). Community structure or function: effects of environmental
stress on benthic macroinvertebrates at different spatial scales. Freshwater Biology, 52(7), 1380-
1399.

Ferreira, W. R., Paiva, L. T., & Callisto, M. (2011). Development of a benthic multimetric index
for biomonitoring of a neotropical watershed. Brazilian Journal of Biology, 71(1), 15-25.

77
Fetscher, A. E., Stancheva, R., Kociolek, J. P., Sheath, R. G., Stein, E. D., Mazor, R. D., ... &
Busse, L. B. (2014). Development and comparison of stream indices of biotic integrity using
diatoms vs. non-diatom algae vs. a combination. Journal of applied phycology, 26(1), 433-450.

Foley, M. M., Armsby, M. H., Prahler, E. E., Caldwell, M. R., Erickson, A. L., Kittinger, J. N.,
... & Levin, P. S. (2013). Improving ocean management through the use of ecological principles
and integrated ecosystem assessments. BioScience, 63(8), 619-631.

Fortuna, M. A., Stouffer, D. B., Olesen, J. M., Jordano, P., Mouillot, D., Krasnov, B. R., ... &
Bascompte, J. (2010). Nestedness versus modularity in ecological networks: two sides of the
same coin?. Journal of Animal Ecology, 811-817.

Fournier, B., Mouly, A., & Gillet, F. (2016). Multiple assembly rules drive the co-occurrence of
orthopteran and plant species in grasslands: combining network, functional and phylogenetic
approaches. Frontiers in plant science, 7, 1224.

Freeman, P. L., & Schorr, M. S. (2004). Influence of watershed urbanization on fine sediment
and macroinvertebrate assemblage characteristics in Tennessee Ridge and Valley
Streams. Journal of Freshwater Ecology, 19(3), 353-362.

Freilich, M. A., Wieters, E., Broitman, B. R., Marquet, P. A., & Navarrete, S. A. (2018). Species
co-occurrence networks: Can they reveal trophic and non-trophic interactions in ecological
communities?. Ecology, 99(3), 690-699.
78

Fridley, J. D., Vandermast, D. B., Kuppinger, D. M., Manthey, M., & Peet, R. K. (2007). Co-
occurrence based assessment of habitat generalists and specialists: A new approach for the
measurement of niche width. Journal of ecology, 95(4), 707-722.

Fuhrman, J. A., & Steele, J. A. (2008). Community structure of marine bacterioplankton:
patterns, networks, and relationships to function. Aquatic Microbial Ecology, 53(1), 69-81.

Garie, H. L., & McIntosh, A. (1986). Distribution of benthic macroinvertebrates in a stream
exposed to urban runoff 1. JAWRA Journal of the American Water Resources Association, 22(3),
447-455.

Gerth, W. J., & Herlihy, A. T. (2006). Effect of sampling different habitat types in regional
macroinvertebrate bioassessment surveys. Journal of the North American Benthological
Society, 25(2), 501-512.

Goldwasser, L., & Roughgarden, J. (1997). Sampling effects and the estimation of food-web
properties. Ecology, 78(1), 41-54.

Goodwin, K. D., Thompson, L. R., Duarte, B., Kahlke, T., Thompson, A. R., Marques, J. C., &
Caçador, I. (2017). DNA sequencing as a tool to monitor marine ecological status. Frontiers in
Marine Science, 4, 107.

79
Gotelli, N. J. (2000). Null model analysis of species co-occurrence patterns. Ecology, 81(9),
2606-2621.

Grömping, U. (2006). Relative importance for linear regression in R: the package
relaimpo. Journal of statistical software, 17(1), 1-27.

Grönroos, M., Heino, J., Siqueira, T., Landeiro, V. L., Kotanen, J., & Bini, L. M. (2013).
Metacommunity structuring in stream networks: roles of dispersal mode, distance type, and
regional environmental context. Ecology and evolution, 3(13), 4473-4487.

Guerry, A. D., Polasky, S., Lubchenco, J., Chaplin-Kramer, R., Daily, G. C., Griffin, R., ... &
Feldman, M. W. (2015). Natural capital and ecosystem services informing decisions: From
promise to practice. Proceedings of the National Academy of Sciences, 112(24), 7348-7355.

Hawkins, C. P., Mykrä, H., Oksanen, J., & Vander Laan, J. J. (2015). Environmental disturbance
can increase beta diversity of stream macroinvertebrate assemblages. Global Ecology and
Biogeography, 24(4), 483-494.

Hawkins, C. P., Olson, J. R., & Hill, R. A. (2010). The reference condition: predicting
benchmarks for ecological and water-quality assessments. Journal of the North American
Benthological Society, 29(1), 312-343.

80
Heino, J. (2013). Environmental heterogeneity, dispersal mode, and co-occurrence in stream
macroinvertebrates. Ecology and Evolution, 3(2), 344-355.

Herbst, D. B., & Silldorff, E. L. (2006). Comparison of the performance of different
bioassessment methods: similar evaluations of biotic integrity from separate programs and
procedures. Journal of the north american Benthological society, 25(2), 513-530.

Hering, D., Borja, A., Jones, J. I., Pont, D., Boets, P., Bouchez, A., ... & Leese, F. (2018).
Implementation options for DNA-based identification into ecological status assessment under the
European Water Framework Directive. Water Research, 138, 192-205.

Hijmans, R. J., Williams, E., & Vennes, C. (2012). geosphere: Spherical Trigonometry. R
package version 1.2–28. CRAN. R-project. org/package= geosphere.

Homer, C., Dewitz, J., Fry, J., Coan, M., Hossain, N., Larson, C., ... & Wickham, J. (2007).
Completion of the 2001 national land cover database for the counterminous United
States. Photogrammetric engineering and remote sensing, 73(4), 337.

Hu, A., Ju, F., Hou, L., Li, J., Yang, X., Wang, H., ... & Yu, C. P. (2017). Strong impact of
anthropogenic contamination on the co-occurrence patterns of a riverine microbial
community. Environmental microbiology, 19(12), 4993-5009.

81
Jackson, M. C., Weyl, O. L. F., Altermatt, F., Durance, I., Friberg, N., Dumbrell, A. J., ... &
Leadley, P. W. (2016). Recommendations for the next generation of global freshwater biological
monitoring tools. In Advances in ecological research (Vol. 55, pp. 615-636). Academic Press.

Jun, Y. C., Kim, N. Y., Kim, S. H., Park, Y. S., Kong, D. S., & Hwang, S. J. (2016). Spatial
distribution of benthic macroinvertebrate assemblages in relation to environmental variables in
Korean nationwide streams. Water, 8(1), 27.

Jun, Y. C., Kim, N. Y., Kwon, S. J., Han, S. C., Hwang, I. C., Park, J. H., ... & Hwang, S. J.
(2011). Effects of land use on benthic macroinvertebrate communities: Comparison of two
mountain streams in Korea. In Annales de Limnologie-International Journal of Limnology (Vol.
47, No. S1, pp. S35-S49). EDP Sciences.

Jun, Y. C., Won, D. H., Lee, S. H., Kong, D. S., & Hwang, S. J. (2012). A multimetric benthic
macroinvertebrate index for the assessment of stream biotic integrity in Korea. International
journal of environmental research and public health, 9(10), 3599-3628.

Kay, G. M., Tulloch, A., Barton, P. S., Cunningham, S. A., Driscoll, D. A., & Lindenmayer, D.
B. (2018). Species co-occurrence networks show reptile community reorganization under
agricultural transformation. Ecography, 41(1), 113-125.

Kerans, B. L., & Karr, J. R. (1994). A benthic index of biotic integrity (B-IBI) for rivers of the
Tennessee Valley. Ecological applications, 4(4), 768-785.
82

King, R. S., Baker, M. E., Kazyak, P. F., & Weller, D. E. (2011). How novel is too novel?
Stream community thresholds at exceptionally low levels of catchment urbanization. Ecological
applications, 21(5), 1659-1678.

Kuhn, M. (2008). Building predictive models in R using the caret package. Journal of statistical
software, 28(5), 1-26.

Lakew, A., & Moog, O. (2015). A multimetric index based on benthic macroinvertebrates for
assessing the ecological status of streams and rivers in central and southeast highlands of
Ethiopia. Hydrobiologia, 751(1), 229-242.

Larsen, S., & Ormerod, S. J. (2014). Anthropogenic modification disrupts species co-occurrence
in stream invertebrates. Global change biology, 20(1), 51-60.

Laurance, W. F., Sayer, J., & Cassman, K. G. (2014). Agricultural expansion and its impacts on
tropical nature. Trends in ecology & evolution, 29(2), 107-116.

Lenat, D. R., & Crawford, J. K. (1994). Effects of land use on water quality and aquatic biota of
three North Carolina Piedmont streams. Hydrobiologia, 294(3), 185-199.

83
Mandakovic, D., Rojas, C., Maldonado, J., Latorre, M., Travisany, D., Delage, E., ... & Cabrera,
P. (2018). Structure and co-occurrence patterns in microbial communities under acute
environmental stress reveal ecological factors fostering resilience. Scientific reports, 8(1), 5875.

Masese, F. O., Raburu, P. O., & Muchiri, M. (2009). A preliminary benthic macroinvertebrate
index of biotic integrity (B-IBI) for monitoring the Moiben River, Lake Victoria Basin,
Kenya. African Journal of Aquatic Science, 34(1), 1-14.

Maxted, J. R., Barbour, M. T., Gerritsen, J., Poretti, V., Primrose, N., Silvia, A., ... & Renfrow,
R. (2000). Assessment framework for mid-Atlantic coastal plain streams using benthic
macroinvertebrates. Journal of the North American Benthological Society, 19(1), 128-144.

Mazor, R. D., Rehn, A. C., Ode, P. R., Engeln, M., Schiff, K. C., Stein, E. D., ... & Hawkins, C.
P. (2016). Bioassessment in complex environments: designing an index for consistent meaning in
different settings. Freshwater Science, 35(1), 249-271.

Mihuc, T. B. (1997). The functional trophic role of lotic primary consumers: generalist versus
specialist strategies. Freshwater biology, 37(2), 455-462.

Morriën, E., Hannula, S. E., Snoek, L. B., Helmsing, N. R., Zweers, H., De Hollander, M., ... &
Duyts, H. (2017). Soil networks become more connected and take up more carbon as nature
restoration progresses. Nature communications, 8(1), 1-10.

84
Morueta-Holme, N., Blonder, B., Sandel, B., McGill, B. J., Peet, R. K., Ott, J. E., ... & Svenning,
J. C. (2016). A network approach for inferring species associations from co-occurrence
data. Ecography, 39(12), 1139-1150.

Moyle, P. B., & Randall, P. J. (1998). Evaluating the biotic integrity of watersheds in the Sierra
Nevada, California. Conservation Biology, 12(6), 1318-1326.

Mykrä, H., & Heino, J. (2017). Decreased habitat specialization in macroinvertebrate
assemblages in anthropogenically disturbed streams. Ecological Complexity, 31, 181-188.

Newman, M. E., & Girvan, M. (2004). Finding and evaluating community structure in
networks. Physical review E, 69(2), 026113.

Nielsen, A., & Bascompte, J. (2007). Ecological networks, nestedness and sampling
effort. Journal of Ecology, 1134-1141.

Novotny, V., Bartošová, A., O’Reilly, N., & Ehlinger, T. (2005). Unlocking the relationship of
biotic integrity of impaired waters to anthropogenic stresses. Water Research, 39(1), 184-198.

Ode, P. R. (2003). CAMLnet: list of California macroinvertebrate taxa and standard taxonomic
effort. Aquatic Bioassessment Laboratory, Rancho Cordova. Retrieved September, 10, 2014.

85
Ourso, R. T., & Frenzel, S. A. (2003). Identification of linear and threshold responses in streams
along a gradient of urbanization in Anchorage, Alaska. Hydrobiologia, 501(1-3), 117-131.

Ovaskainen, O., Hottola, J., & Siitonen, J. (2010). Modeling species co-occurrence by
multivariate logistic regression generates new hypotheses on fungal interactions. Ecology, 91(9),
2514-2521.

Paerl, H. W., Gardner, W. S., Havens, K. E., Joyner, A. R., McCarthy, M. J., Newell, S. E., ... &
Scott, J. T. (2016). Mitigating cyanobacterial harmful algal blooms in aquatic ecosystems
impacted by climate change and anthropogenic nutrients. Harmful Algae, 54, 213-222.

Paul, M. J., & Meyer, J. L. (2001). Streams in the urban landscape. Annual review of Ecology
and Systematics, 32(1), 333-365.

Pearman, P. B. (1997). Correlates of Amphibian Diversity in an Altered Landscape of
Amazonian Ecuador: Correlaciones de la Diversidad de Anfibios en un Paisaje Alterado de la
Amazonía Ecuatoriana. Conservation Biology, 11(5), 1211-1225.

Peck, D. V., Herlihy, A. T., Hill, B. H., Hughes, R. M., Kaufmann, P. R., Klemm, D. J., ... &
Magee, T. (2006). Environmental monitoring and assessment program—surface waters western
pilot study: field operations manual for wadeable streams. EPA 600/R-06/003.

86
Poisot, T., & Gravel, D. (2014). When is an ecological network complex? Connectance drives
degree distribution and emerging network properties. PeerJ, 2, e251.

Pratt, B., & Chang, H. (2012). Effects of land cover, topography, and built structure on seasonal
water quality at multiple spatial scales. Journal of hazardous materials, 209, 48-58.

Rawer-Jost, C., Böhmer, J., Blank, J., & Rahmann, H. (2000). Macroinvertebrate functional
feeding group methods in ecological assessment. Hydrobiologia, 422, 225-232.

Team, R Core (2018). R: A language and environment for statistical computing. R Foundation
for Statistical Computing, Vienna, Austria.

Rehn, A. C., Ode, P. R., & Hawkins, C. P. (2007). Comparisons of targeted-riffle and reach-wide
benthic macroinvertebrate samples: implications for data sharing in stream-condition
assessments. Journal of the North American Benthological Society, 26(2), 332-348.

Reynoldson, T. B., Norris, R. H., Resh, V. H., Day, K. E., & Rosenberg, D. M. (1997). The
reference condition: a comparison of multimetric and multivariate approaches to assess water-
quality impairment using benthic macroinvertebrates. Journal of the North American
Benthological Society, 16(4), 833-852.

87
Richards, A. B., & Rogers, D. C. (2006). List of freshwater macroinvertebrate taxa from
California and adjacent states including standard taxonomic effort levels. Southwest Association
of Freshwater Invertebrate Taxonomists (SAFIT).

Ruan, Q., Dutta, D., Schwalbach, M. S., Steele, J. A., Fuhrman, J. A., & Sun, F. (2006). Local
similarity analysis reveals unique associations among marine bacterioplankton species and
environmental factors. Bioinformatics, 22(20), 2532-2538.

Schmeller, D. S., Julliard, R., Bellingham, P. J., Böhm, M., Brummitt, N., Chiarucci, A., ... &
Gregory, R. D. (2015). Towards a global terrestrial species monitoring program. Journal for
Nature Conservation, 25, 51-57.

Seaber, P. R., Kapinos, F. P., & Knapp, G. L. (1987). Hydrologic unit maps.

Shaver, E., Maxted, J., Curtis, G., & Carter, D. (1995). Watershed protection using an integrated
approach. In Stormwater NPDES related monitoring needs (pp. 435-459). ASCE.

Shostell, J. M., & Williams, B. S. (2007). Habitat complexity as a determinate of benthic
macroinvertebrate community structure in cypress tree reservoirs. Hydrobiologia, 575(1), 389-
399.

88
Silva, D. R., Herlihy, A. T., Hughes, R. M., & Callisto, M. (2017). An improved
macroinvertebrate multimetric index for the assessment of wadeable streams in the neotropical
savanna. Ecological Indicators, 81, 514-525.

Simons, A. L., Mazor, R., Stein, E. D., & Nuzhdin, S. (2019). Using alpha, beta, and zeta
diversity in describing the health of stream-based benthic macroinvertebrate
communities. Ecological Applications, 29(4), e01896.

Smith, R. F., & Lamp, W. O. (2008). Comparison of insect communities between adjacent
headwater and main-stem streams in urban and rural watersheds. Journal of the North American
Benthological Society, 27(1), 161-175.

Snodgrass, J. W. (1997). Temporal and spatial dynamics of beaver-created patches as influenced
by management practices in a south-eastern North American landscape. Journal of Applied
Ecology, 1043-1056.

Socolar, J. B., Gilroy, J. J., Kunin, W. E., & Edwards, D. P. (2016). How should beta-diversity
inform biodiversity conservation?. Trends in ecology & evolution, 31(1), 67-80.

Sponseller, R. A., Benfield, E. F., & Valett, H. M. (2001). Relationships between land use,
spatial scale and stream macroinvertebrate communities. Freshwater biology, 46(10), 1409-
1424.

89
Steele, J. A., Countway, P. D., Xia, L., Vigil, P. D., Beman, J. M., Kim, D. Y., ... & Rose, J. M.
(2011). Marine bacterial, archaeal and protistan association networks reveal ecological
linkages. The ISME journal, 5(9), 1414-1425.

Steenweg, R., Hebblewhite, M., Kays, R., Ahumada, J., Fisher, J. T., Burton, C., ... & Brodie, J.
(2017). Scaling-up camera traps: Monitoring the planet's biodiversity with networks of remote
sensors. Frontiers in Ecology and the Environment, 15(1), 26-34.

Stein, E. D., White, B. P., Mazor, R. D., Jackson, J. K., Battle, J. M., Miller, P. E., ... &
Sweeney, B. W. (2014). Does DNA barcoding improve performance of traditional stream
bioassessment metrics?. Freshwater Science, 33(1), 302-311.

Stepenuck, K. F., Crunkilton, R. L., & Wang, L. (2002). Impacts of urban landuse on
macroinvertebrate communities in southeastern Wisconsin streams 1. JAWRA Journal of the
American Water Resources Association, 38(4), 1041-1051.

Stoll, S., Breyer, P., Tonkin, J. D., Früh, D., & Haase, P. (2016). Scale-dependent effects of river
habitat quality on benthic invertebrate communities—implications for stream restoration
practice. Science of the Total Environment, 553, 495-503.

Suga, C. M., & Tanaka, M. O. (2013). Influence of a forest remnant on macroinvertebrate
communities in a degraded tropical stream. Hydrobiologia, 703(1), 203-213.

90
Sweeney, B. W., Battle, J. M., Jackson, J. K., & Dapkey, T. (2011). Can DNA barcodes of
stream macroinvertebrates improve descriptions of community structure and water
quality?. Journal of the North American Benthological Society, 30(1), 195-216.

Thébault, E., & Fontaine, C. (2010). Stability of ecological communities and the architecture of
mutualistic and trophic networks. Science, 329(5993), 853-856.

Tulloch, A. I., Chadès, I., Dujardin, Y., Westgate, M. J., Lane, P. W., & Lindenmayer, D. (2016).
Dynamic species co-occurrence networks require dynamic biodiversity
surrogates. Ecography, 39(12), 1185-1196.

Tylianakis, J. M., Laliberté, E., Nielsen, A., & Bascompte, J. (2010). Conservation of species
interaction networks. Biological conservation, 143(10), 2270-2279.

Vander Laan, J. J., Hawkins, C. P., Olson, J. R., & Hill, R. A. (2013). Linking land use, in-
stream stressors, and biological condition to infer causes of regional ecological impairment in
streams. Freshwater Science, 32(3), 801-820.

Veech, J. A. (2013). A probabilistic model for analysing species co-occurrence. Global Ecology
and Biogeography, 22(2), 252-260.

91
Vieira, C., Séneca, A., Sérgio, C., & Ferreira, M. T. (2012). Bryophyte taxonomic and functional
groups as indicators of fine scale ecological gradients in mountain streams. Ecological
Indicators, 18, 98-107.

Vile, J. S., & Henning, B. F. (2018). Development of indices of biotic integrity for high-gradient
wadeable rivers and headwater streams in New Jersey. Ecological Indicators, 90, 469-484.

von der Ohe, P. C., & Goedkoop, W. (2013). Distinguishing the effects of habitat degradation
and pesticide stress on benthic invertebrates using stressor-specific metrics. Science of the Total
Environment, 444, 480-490.

Vörösmarty, C. J., McIntyre, P. B., Gessner, M. O., Dudgeon, D., Prusevich, A., Green, P., ... &
Davies, P. M. (2010). Global threats to human water security and river
biodiversity. nature, 467(7315), 555-561.

Voß, K., & Schäfer, R. B. (2017). Taxonomic and functional diversity of stream invertebrates
along an environmental stress gradient. Ecological Indicators, 81, 235-242.

Waite, I. R., Brown, L. R., Kennen, J. G., May, J. T., Cuffney, T. F., Orlando, J. L., & Jones, K.
A. (2010). Comparison of watershed disturbance predictive models for stream benthic
macroinvertebrates for three distinct ecoregions in western US. Ecological Indicators, 10(6),
1125-1136.

92
Wallace, A. M., & Biastoch, R. G. (2016). Detecting changes in the benthic invertebrate
community in response to increasing chloride in streams in Toronto, Canada. Freshwater
Science, 35(1), 353-363.

Weigel, B. M., & Dimick, J. J. (2011). Development, validation, and application of a
macroinvertebrate-based Index of Biotic Integrity for nonwadeable rivers of Wisconsin. Journal
of the North American Benthological Society, 30(3), 665-679.

Widder, S., Besemer, K., Singer, G. A., Ceola, S., Bertuzzo, E., Quince, C., ... & Battin, T. J.
(2014). Fluvial network organization imprints on microbial co-occurrence
networks. Proceedings of the National Academy of Sciences, 111(35), 12799-12804.

Williams, R. J. (2011). Biology, methodology or chance? The degree distributions of bipartite
ecological networks. PLoS One, 6(3), e17645.

Yan, G., Martinez, N. D., & Liu, Y. Y. (2017). Degree heterogeneity and stability of ecological
networks. Journal of the Royal Society Interface, 14(131), 20170189.

Zalack, J. T., Smucker, N. J., & Vis, M. L. (2010). Development of a diatom index of biotic
integrity for acid mine drainage impacted streams. Ecological Indicators, 10(2), 287-295.

93
Zhang, Y., Liu, L., Cheng, L., Cai, Y., Yin, H., Gao, J., & Gao, Y. (2014). Macroinvertebrate
assemblages in streams and rivers of a highly developed region (Lake Taihu Basin,
China). Aquatic biology, 23(1), 15-28.

94
Chapter 3
Zeta diversity patterns in metabarcoded lotic communities as a tool for bioassessment
Abstract
Assessments of the ecological health of biological communities in streams have
traditionally focused on measures of the local diversity of morphologically classified benthic
macroinvertebrates (BMIs). Such communities are often connected through various ecological
processes, such as dispersal, and may be more accurately assessed as components of regional
scale metacommunities. With recent declines in the costs of sequencing and computation, it has
become increasingly feasible to investigate changes in patterns of landscape diversity across
metacommunities of BMIs and other groups as a means of regional scale bioassessment. Here we
investigate the use of zeta diversity as a novel means of quantifying landscape diversity for a
variety of stream communities as described using a variety of metabarcoding primers. Using
DNA extracted from 148 stream samples in California we used various orders of zeta diversity to
construct models of biotic integrity, as described by the California Streams Condition Index
(CSCI), for communities of prokaryotes, soft algae, diatoms, and BMIs. We found our models
could account for most of the observed variation in the CSCI using zeta diversity patterns of
BMIs, soft algae, and diatoms, but not prokaryotes. We also observed that aggregating
communities at higher taxonomic levels to be significantly associated with both improved model
performances and a greater likelihood of community assembly being described via a process of
niche assemble versus a stochastic one.
95
Introduction
The biological activity of stream communities provides for a variety of ecosystem
services ranging from the breakdown of organic matter in runoff (Schäfer et al., 2012),
sustaining the health of fisheries (Colvin et al., 2019), as well as regulating the formation of algal
blooms (Castilla et al., 2015). Assessments of the biological condition of streams have typically
involved the determination of the composition of specific communities, such as benthic
macroinvertebrates (BMIs) (Jun et al., 2012; Paller et al., 2017) or fish (Li et al., 2018; Vile et
al., 2018), through morphological identification. However, the process of sorting species by
morphology frequently requires a large investment of time and money and may yet still be
limited by mismatches between unique categories of biological taxa and morphotaxa (Herbert et
al., 2004; Pauls et al., 2010; Ristau et al., 2013; Brown et al., 2015; Lanzén et al., 2016).
Morphologically based methods are further limited in describing the composition of microbial
communities, which are significant factors in driving biogeochemical cycles in stream
environments (Wakelin et al., 2008; Battin et al., 2016; Ibekwe et al., 2016). These limitations
point to significant uncertainties in the stream communities described by morphologically based
methods, and the confidence of bioassessment methods based upon them.
With recent declines in the cost of genetic sequencing DNA metabarcoding has
demonstrated its potential as the basis for new bioassessment tools which are unencumbered by
the limitations of morphologically based methods (Ratnasingham & Hebert 2007; Hajibabaei et
al., 2011; Taberlet et al., 2012; Ji et al., 2013; Pawlowski et al., 2018). This technique involves
the use of high throughput next generation sequencing of DNA fragments, which are then
classified using databases of taxonomically identified reference sequences (Elbrecht & Leese,
2017; Keck et al., 2017; Banerji et al., 2018). Using metabarcoding the composition of a variety
96
of communities typically used in the bioassessment of freshwater streams have been evaluated
including BMIs (Hajibabaei et al., 2011; Carew et al., 2013; Stein et al., 2014a; Elbrecht & Leese
2017), soft algae (Banerji et al., 2018; Minerovic et al., 2020; Wolf & Vis, 2020), diatoms
(Vasselon et al., 2017; Bailet et al., 2019; Valentin et al., 2019), and prokaryotes (Chen et al.,
2018; Cordier et al., 2019). Use of metabarcoding based methods have demonstrated its potential
for both improved detection efficiency of rare organisms, as well as finer taxonomic resolution
(Baird & Sweeney 2011; Sweeney et al., 2011; Stein et al., 2013; Zhan et al., 2013; Gibson et al.,
2015), indicating its potential in forming the basis of far more sensitive and nuanced
bioassessments than those derived from traditional methods (Stein et al., 2014b; Elbrecht &
Leese, 2016; Mortágua et al., 2019).
For metabarcoding the use of a variety of gene fragments have been explored in
organizing genetic samples from the environment into the composition of various communities.
For studies of eukaryotic communities one of the more widely applied gene fragments of interest
is the V9 hypervariable region of the 18S rDNA gene, which can be amplified with a single set
of primers (Amaral-Zettler et al., 2009). The 18S-V9 region has been used in metabarcoding
studies of freshwater communities ranging from BMIs (Xie et al., 2017), soft algae (Abad et al.,
2016), and diatoms (Guo et al., 2016). For diatoms further improvements in both the
completeness and taxonomic resolution of communities have been found with the use of the
partial plastid rbcL gene (Kermarrec et al., 2014), especially for bioassessments of freshwater
environments (Hering et al., 2018). For studies of the composition of prokaryotic communities
the V4 hypervariable region of the 16S rDNA gene (Gilbert et al., 2014) has been frequently
used in a variety of environments, including streams (Paruch et al., 2019). Although not as
widely used as the 16S-V4 region, use of the 16S region spanning positions 515 to 926 has also
97
been investigated as a potentially less biased gene fragment for reproducing the composition of
known prokaryotic communities (Parada et al., 2016).
For communities quantified either through morphological or metabarcoding methods the
focus of bioassessments in streams has typically been based on various measures of alpha
diversity. For example, the degree of ecological degradation experienced in stream communities
within California have been assessed through a composite measure of taxonomic and functional
alpha diversity known as the California Stream Condition Index (CSCI; Rehn et al. 2015, Mazor
et al. 2016). However, individual stream communities do not exist in an ecological vacuum.
Stream communities, being connected through mechanisms such as dispersal, may be more
accurately assessed as part of regional metacommunities (Vanschoenwinkel et al., 2007;
Grönroos et al. 2013, Socolar et al. 2016). To provide a more complete form of bioassessment
for stream communities it is important to then consider both local and regional measurements of
their biodiversity (Heino, 2013; Socolar et al. 2016). To then assess potential relationships
between environmental gradients and metacommunity diversity various investigations have been
conducted using measures, such as various types of beta diversity, of landscape diversity (Ruhí et
al., 2017; Gál et al., 2019). Though beta diversity is a well-established ecological concept it is
somewhat limited in describing changes to metacommunities. First, as beta diversity involves
only a pairwise comparison of intercommunity diversity it cannot fully describe diversity
patterns in groups involving more than two communities (Chao et al. 2008). Second, although
various measures of average beta diversity have been developed they are potentially biased
towards the presence of categories of organisms which found only in a few sampled communities
(Beck et al., 2013; Latombe et al. 2017a).
98
To address limitations involved in relying solely on alpha or beta diversity, a new
measure, zeta diversity, has been recently developed (Figure 1). Zeta diversity is a generalized
extension of established measures of biodiversity, which was developed to describe the number
of unique categories of organisms held in common between an arbitrarily large set of
communities (Hui & McGeoch, 2014). In general zeta diversity, denoted as 𝞯N, quantifies the
average number of unique categories of organisms held in common between N communities
within a group of N or more communities. In being able to directly compare the composition of
large groups of communities use of zeta diversity can avoid limitations involved with beta
diversity. First, comparisons of intercommunity diversity more complex than pairwise measures
can be calculated. Second, in enabling calculations of the number of unique categories of
organisms held in common between three or more communities, zeta diversity can more
accurately assess changes in the presence of those unique organisms shared between a large
number of communities. Given prior evidence of the loss of commonly distributed organisms
acting as a reliable indicator of widespread environmental degradation (Gaston & Fuller 2008;
Pond 2012; Baker et al., 2019) this indicates the potential use for zeta diversity to form the basis
of bioassessments on regional scales. Further support for zeta diversity’s potential in forming the
basis for bioassessments of metacommunities has been demonstrated in recent work involving its
use in quantifying changes in communities in environments as varied as islands (Leihy et al.
99
2018), streams (Simons et al., 2019),
protected areas in aquatic
environments (Britton et al. 2017), and
forests (Lazarina et al., 2019).

Analysis of zeta diversity
patterns in stream metacommunities
can also be used to address a number
of ecological questions with relevance
in bioassessments, such as determining
the likeliest model of community
assembly. Prior searches for
communities of indicator taxa have
shown the importance of a diverse array of ecological niches, and thus environmental
sensitivities (Hilty and Merenlender, 2000). This suggests biological communities which tend to
follow a niche differentiation model of assembly, as opposed to a stochastic process, would be
more reliable as indicators in any future bioassessment framework. In the development of zeta
diversity, it was discovered that the value of 𝞯N in a metacommunity tends to decline, with
sample order N, following either as an exponential or power-law function. The likelihood of
these models, exponential or power-law, are respectively associated with either a stochastic or
niche differentiation process of community assembly (Hui & McGeoch, 2014; McGeoch et al.,
2019).
Figure 1: An illustration of the first four orders of zeta diversity.
The mean number of unique categories of organisms per
sample, a measure of alpha diversity, is represented by the
value 𝞯 1. In comparing two or more samples the average
number of unique categories of organism held in common
between any two samples is represented by the value 𝞯 2. This
process can be extended to N samples, allowing for a
determination of the values 𝞯 1 through 𝞯 N.
100
Additional use of the zeta diversity framework has also been used to quantify the roles
geographic and environment distances underlying intercommunity differences (Latombe et al.,
2017a; McGeoch et al., 2019). Both geographic distance and environmental conditions are
known to influence patterns of dispersal which connect stream metacommunities (Patrick and
Swan, 2011; Heino, 2013). This suggests additional potential in integrating zeta diversity into the
development of new bioassessments of streams on a regional scale.
Evidence of relationships between environmental gradients and measures of landscape
diversity has been found across a diverse array of biological communities found in streams
ranging from soft algae (Passy and Blanchet, 2007; Wu et al., 2018), diatoms (Soininen et al.,
2004; Jyrkänkallio‐Mikkola et al., 2016), BMIs (Costa and Melo, 2008; Maloney et al., 2011),
and prokaryotes (Besemer, 2015; Bandh et al., 2019). Using the zeta diversity framework, we
were also able to further demonstrate the use of a novel measure of landscape diversity in
modeling ecological degradation of streams on a regional scale (Simons et al., 2019). The goal of
this study is as follows:
(1) Classify stream metacommunities, using a variety of metabarcoding primer sets, into
groups of prokaryotes, diatoms, soft algae, and BMIs.
(2) Test the hypothesis that patterns in zeta diversity across stream metacommunities, can be
used as reliable indicators of ecological degradation at a regional scale as described by
the CSCI.
(3) Assess the likeliest models of community assembly for these stream metacommunities.
101
Materials and Methods
Sample scope
The samples used in this study were collected at wadeable stream reaches located
throughout California as part of four programs: the California State Water Board Surface Water
Ambient Monitoring Program (SWAMP) Perennial Stream Assessment survey and Reference
Condition Management Program, the Stormwater Monitoring Coalition, the San Diego Regional
Water Board, and the Bay Area Regional Monitoring Coalition. The data used in this study were
taken from 148 samples, covering 122 unique geographic locations, during the period 2016-
2017. Each sample contains the following data: BMIs enumerated and sorted by morphology to a
standardized level (generally to genus with the exception of chironomids, which were identified
to subfamily; Richards and Rogers, 2006), BMIs, soft algae, and diatoms classified using two
eukaryotic primer sets, prokaryotes classified using two primer sets, site altitude in meters, the
percent of total developed land use (agricultural, urban, and managed landscapes) within a 5 km
clipped buffer of the watershed upstream of the sampling site, and a CSCI score of biotic
integrity.
Acquiring and calculating data for the CSCI
For each sampled community biotic integrity was calculated using the CSCI, our current
standard score for the biotic integrity of stream communities. The CSCI compares
morphologically sorted BMI communities and metrics from observed samples to corresponding
values expected under undisturbed reference conditions as described by a site’s local
environment, such as geology, climate, and watershed area (Mazor et al., 2016). In obtaining
102
BMI communities approximately 55% were sampled using a reach-wide protocol (Peck et al.,
2006), with the remaining samples acquired using a targeted riffle protocol. Both sampling
methods have been found to produce comparable BMI community data (Gerth & Herlihy, 2006;
Herbst & Silldorff, 2006; Rehn et al., 2007). Once each BMI was sorted to a standardized level,
typically genus, it was then assigned to one of eight functional feeding groups using CAMLnet
(Ode, 2003).
To calculate the CSCI two sets of measurements are used which are based on BMI
communities as described by a standardized taxonomy (Richards and Rogers, 2006): the first is
a ratio of observed to expected BMIs (O/E), and the second is a predicted multimetric index
(pMMI) based on six measures of the sampled BMI community’s functional and ecological
structure. The CSCI was designed to minimize the impact of major natural gradients, which has
allowed for it to be a consistently useful measure of biological conditions in a variety of
environmental settings (Reynoldson et al., 2006; Hawkins et al., 2010).
Acquiring samples for metabarcoding
Samples collected for DNA analysis were done so using the SWAMP algal collection
method (Ode et al., 2016). This method entails the collection of 11 algal samples from various
substrates, such as sand or pebbles, across a 150 m stream reach transect. These samples are then
composited into a single representative sample in proportion of their substrate type to those
found within the transect.
For each algal sample between 10 and 50 ml of material was filtered onto 0.2 µm
Whatman Nucleopore polycarbonate filters (GE Healthcare Life Sciences, Buckinghamshire,
UK), with between one and three such filters collected from each composite sample. The sample
volume was determined on both its biomass density and the ability to filter it from the
103
environment. Each filter was first preserved in bead solution storage buffer (Mo Bio
Laboratories, Inc., Carlsbad, CA) and then stored at -80C. Following a thaw at 4C, DNA
extractions were then performed on the filters using PowerLyzer PowerSoil DNA kit and
protocols (Mo Bio Laboratories, Inc., Carlsbad, CA). Following extractions DNA yield was
quantified using a Nanodrop 8000 (Thermo Scientific, Wilmington, DE, USA).
DNA metabarcoding and bioinformatics
We targeted a variety of eukaryotic communities (BMIs, soft algae, and diatoms) using a
18S-V9 primer (Amaral-Zettler et al., 2009), diatom communities as described using a rbcL
primer (Vasselon et al., 2017), and prokaryotic communities using both a 16S-V4 primer
(Apprill et al., 2015) and one which targets the 515F-926R region of the 16S rRNA gene (Parada
et al., 2016) (Table 1). Paired-end sequencing of all of the resulting barcoded amplicons was
performed on an Illumina MiSeq platform by Laragen, Inc. (Culver City, CA, USA) using either
2x250 reads (rbcL and 16S 515F-926R) or 2x150 reads (16S-V4 and 18S-V9). The resulting
DNA amplicons were processed and demultiplexed using the Illumina MiSeq Recorder, which
enabled removal of barcodes, primers, and adapter sequences. The processed amplicons are
stored as BioProject number PRJNA545290 on the Sequence Read Archive
(https://www.ncbi.nlm.nih.gov/sra).
Table 1: Number of sequences used from each primer for each community.
18S-V9
soft algae
18S-V9
Diatoms
18S-V9
BMIs
16S-V4
Prokaryotes
16S
515F-926R
Prokaryotes
rbcL
Diatoms
104
Sequences 1,390,964 649,314 188,111 6,455,554 3,207,578 869,471

DNA sequences were then processed using QIIME2 2019.07 (Bolyen et al., 2019,
http://qiime2.org). In our workflow forward and reverse reads were first merged into contigs.
These contigs were then filtered by first removing chimeric sequences, then removing reads
which occurred less than two times in any individual sample. We then generated Amplicon
Sequence Variants (ASVs) within QIIME2 using DADA2 (Callahan et al., 2016). Taxonomic
assignments were then performed for 16S and 18S ASVs by using BLAST within QIIME2
against the SILVA v132 database (Quast et al., 2013). For ASVs classified using the rbcL primer
BLAST was used against the R-Syst database (Rimet et al., 2016) within QIIME2 in order to
perform diatom specific taxonomic assignments. ASV counts were then converted into a
percentage of total reads per sample as a means of normalizing differential read counts.
Rarefaction was not used on read counts as it has been found to increase the likelihood of false
positives (McMurdie and Holmes, 2014). Code documenting the workflow to produce ASV
tables are available online (https://github.com/stheroux/AlgaeDNA).
Additional processing of 18S ASVs was then run in order to assign amplicons to
communities of diatoms, soft algae, and BMIs. ASVs were assigned to BMI communities based
on shared taxonomies with BMIs used in calculating the CSCI (Richards and Rogers, 2006). For
communities described using the 18S-V9 primer we classified reads within the phylum
Ochrophyta as diatoms and those within the phyla Charophyta or Chlorophyta as soft algae.
105
Land use
To determine the geographic distribution of land use upstream of each sampling site we
used National Land Cover Data (NLCD) acquired in the year 2011 (Homer et al., 2015). We
define land use as the total percent of land cover in a designated area dedicated to urbanization,
agriculture, and managed landscapes such as parks. The area used in this calculation is defined
by a watershed-clipped buffer running 5 km upstream from the sampling site as determined with
ArcGIS tools (version 10.3; Environmental Systems Research Institute) (Mazor et al., 2016).
Sites were then classified into one of three categories of disturbance based on land use:
undisturbed (0-3%), intermediate (3-15%), and disturbed (15-100%). It should be noted that land
use at our sample site locations did not vary significantly over the duration of this study (r = -
0.13, p = 0.13).
Sample group selection
To analyze patterns in our measures of zeta diversity we first constructed equally sized
groups of 25 samples by random subsampling sites within each disturbance category. To
generate data to analyze trends in zeta diversity this process was repeated 100 times for sample
sites within each of the three disturbance categories. For each sample group the following values
were calculated: the mean percent developed upstream land use (land use), mean sample site
altitude (altitude), and mean geographic separation distance between samples (distance).
Distance between sample sites was calculated, in meters, used the function distm within the R
package geosphere (Hijmans, 2017).
106
Calculations of 𝞯 diversity
For this study we focused on three orders of zeta diversity to describe changes in both
local and landscape diversity for our sampled communities (Simons et al., 2019). To describe
changes to diversity within communities we used 𝞯1, which measures the mean alpha diversity
per group of samples. To describe changes to landscape diversity we used both 𝞯2 and 𝞯10, which
represent the mean number of unique taxonomic groups held in common across 2 and 10
sampled communities respectively. All measures of zeta diversity were calculated using the R
package ZETADIV (Latombe et al., 2017b).
To investigate how the similarities between sampled communities within groups decline
with zeta diversity order (N) we tested two decay models, a power-law of the form 𝞯N = 𝞯1N
-b

and an exponential of the form 𝞯N = 𝞯1e
b(N-1)
. The value of the exponent for our zeta diversity
decline models (b) is a measure of intercommunity dissimilarity within our sample groups, with
a greater value indicating a more rapid decline in similarities between communities. Both models
were fitted using the first 10 orders of zeta diversity, 𝞯1 through 𝞯10, for each group of 25
samples. Only the first 10 orders of zeta diversity were chosen as it was found to be sufficient to
demonstrate an asymptotic level of similarity in the average number of unique taxonomic groups
held in common across our sample groups. The likelihood of each model was assessed using an
Akaike Information Criterion (AIC) score within ZETADIV.
Modeled indices of biotic integrity
Linear models of the biotic integrity of a set of samples, as described using the mean
CSCI score per sample group, were constructed using the values of 𝞯1, 𝞯2, 𝞯10, and the ‘lm’
function within the R package stats (v3.5.1 Team, 2018). To determine the relative importance of
107
each order of zeta diversity in our modeled CSCI scores, and to adjust for any collinearity
between measures as a result, the function ‘calc.relimp’ was used within the relaimpo R package
(Grömping, 2006). The relative importance of land use, altitude, and distance in describing
variations in linear models of both the mean CSCI score and our modeled CSCI scores were also
done using the ‘calc.relimp’ function. Further testing of our modeled CSCI scores was done
using a 10‐fold cross‐validation, with 100 repeats, within the ‘train’ function of the R package
‘caret’ (Kuhn, 2008).
Results
For this project we obtained DNA from 148 California stream samples, classifying the
composition of their communities using a variety of metabarcoding primers. Using measures of
zeta diversity of those sampled communities, we constructed linear models of their biotic
integrity at a regional scale. We found partial support for our hypothesis, with strong correlations
between our modeled indices derived from zeta diversity patterns of communities of BMIs, soft
algae, and diatoms and the regional scale averages of biotic integrity as described by the CSCI
score (Table 2).
Table 2: 10-fold cross-validated Pearson correlations between mean and modeled CSCI using communities
aggregated to various taxonomic levels (p < 10
-4
).
18S-V9
Soft algae
18S-V9
Diatoms
18S-V9
BMIs
16S-V4
Prokaryotes
16S
515F-926R
Prokaryotes
rbcL
Diatoms
Species 0.76 0.52 0.50 0.63 0.48 0.84
108
Genus 0.83 0.25 0.74 0.10 0.12 0.96
Family 0.84 0.56 0.66 0.36 0.39 0.96
Order 0.93 0.63 0.61 0.65 0.72 0.97

The strength of these correlations varied by metabarcoding primer, community type, and the
taxonomic level at which community members were aggregated with diatoms, classified using a
rbcL primer, providing the best results at all levels of taxonomic aggregation. We also were able
to use models of how zeta diversity decays with zeta order to determine the likeliest process of
assembly for all of our community types. With this analysis we observed community assembly
models following a process of niche differentiation to be more probable than a stochastic one at
higher levels of taxonomic aggregation (Table 3).
Table 3: AIC scores of Exponential (Power-law) models of zeta diversity decline with sampling order (p < 10
-4
).
18S-V9
Soft algae
18S-V9
Diatoms
18S-V9
BMIs
16S-V4
Prokaryotes
16S
515F-926R
Prokaryotes
rbcL
Diatoms
Species -8.1 (17.1) -15.4 (-12.2) -19.0 (2.3) -2.1 (-24.7) -0.7 (-18.4) 1.7 (-24.3)
Genus -9.1 (-7.9) -18.4 (-24.9) -19.7 (2.4) -13.0 (-63.3) -11.8 (-58.6) -9.0 (-47.4)
Family -12.9 (-5.8) -16.1 (-36.1) -11.4 (-21.3) -23.7 (-73.9) -19.3 (-59.3) -16.1 (-55.2)
Order -14.1 (-24.3) -22.7 (-43.2) -16.9 (-25.4) -31.1 (-77.6) -25.2 (-62.6) -19.8 (-48.7)

109
Sequencing statistics for metabarcoded amplicons
We found a similar number of sequences for diatoms classified through both the rbcL and
18S-V9 primers (Table 1), although the resolution of taxonomic assignments was far higher
using a rbcL primer (Figure 2A). Using
the 18S-V9 primer we identified
approximately an order of magnitude
more soft algal sequences and ASVs
than those for BMIs, although we could
taxonomically identify a far higher
portion of ASVs for BMIs than for soft
algae. Of the two prokaryotic primers
used the 16S-V4 provided approximately
twice as many sequences as with the 16S
515F-926R primer. For prokaryotes,
using either primer, we found taxonomic
identifications could only be provided for at least 50% of ASVs at the level of family or higher.
Using the rbcL primer we were able to identify a far larger number of unique species of
diatoms than with the 18S-V9 primer (Figure 2B). With the 18S-V9 primer we were also able to
identify a three to five times the number of unique categories of organisms for soft algal rather
than communities across various taxonomic levels. For both prokaryotic primers we could
identify a similar number of unique categories of organism from the level of species through
order.
Figure 2: (A) The fraction of ASVs from each primer and
community, and the number of ASVs resolved to various
taxonomic levels. (B) The number of unique groups resolved to
various taxonomic levels.
110
Linear models of diversity and biotic integrity
Using three parameters, 𝞯1, 𝞯2, and 𝞯10, linear models were constructed to predict the
mean value of the CSCI score for various communities aggregated at taxonomic levels ranging
from species to order. The proportion of the observed variation in the mean value of the CSCI
score for communities within a set of samples ranged from 10%, for prokaryotic communities
classified using
the 16S-V4
primer and
aggregated to
genus, to 97%,
for diatom
communities
classified using
the rbcL primer
and aggregated to
order (Table 2). We observed that linear models which could consistently account for the
greatest proportion of variation in the mean CSCI score were constructed using diatom
communities classified using the rbcL primer (Table 2). We also observed that the modeled
indices of the health of stream communities constructed using these diatom communities were
found to vary most similarly in accordance with land use, altitude, and distance as the mean
CSCI score across all of the taxonomic levels considered (Figures 3A-D).
Figure 3: Pearson correlations between environmental variables, measures of zeta
diversity, mean CSCI scores, and modeled indices using communities of diatoms
classified using the rbcL primer and aggregated to various taxonomic levels (p < 10
-4
).
111
𝞯 diversity patterns and the environment
In our analysis of 100 permutations of regional groups of samples we observed the
following trends between our measures of local (𝞯1) and landscape (𝞯2 and 𝞯10) diversity and the
environment (land use, altitude, and distance).
𝞯 diversity and land use
We
observed four
significant trends
between land use
and our measures
of local and
landscape
diversity. Land use
was found to be
associated with a decline in local diversity in communities of BMIs (Figures 4A-D), prokaryotes
classified using the 16S 515F-926R primer and aggregated above the level of species (Figures
5B-D), and diatoms as classified by the rbcL primer (Figures 3A-D). For diatoms classified using
Figure 4: Pearson correlations between environmental variables, measures of zeta
diversity, mean CSCI scores, and modeled indices using BMI communities classified
using the 18S-V9 primer and aggregated to various taxonomic levels (p < 10
-4
).
112
a 18S-V9 primer local diversity was positively correlated with land use (Figures 6B-D), although
this trend was only significant with communities aggregated above the level of species. No
significant trends were observed between land use and the local diversity for communities of soft
algae (Figures 7A-D) and prokaryotes classified using the 16S-V4 primer (Figures 8A-D).
Landscape diversity was found to be negatively correlated with land use in communities of
prokaryotes aggregated to the level of species or order (Figures 5A, 5D, 8A, 8D), BMIs
aggregated to the level of family or order (Figures 4C-D), and diatoms as classified by a rbcL
primer (Figures 3A-D). For soft algae (Figures 7A-D) and diatoms (Figures 6B-D) aggregated
above the level of species, classified using a 18S-V9 primer, landscape diversity was found to be
positively correlated
with land use.
Figure 5: Pearson correlations between environmental variables, measures of zeta
diversity, mean CSCI scores, and modeled indices using prokaryotic communities
classified using the 16S 515F-926R primer and aggregated to various taxonomic levels (p
< 10
-4
).
113
𝞯 diversity and altitude
With
altitude we found
four significant
trends in our
measures of local
and landscape
diversity, although
not across all of our
community types.
We found altitude to be positively correlated with the local diversity of diatoms classified using
the rbcL primer (Figures 3A-D), soft algal communities (Figures 7A, B, D), BMIs (Figures 4A-
D), prokaryotes
classified using the
16S 515F-926R
primer aggregated
below the level of
order (Figures 5A-
C), and
prokaryotes
classified using the
16S-V4 primer aggregated at the level of species or family (Figures 8A-C). For communities of
diatoms classified using the 18S-V9 primer we found a negative correlation between altitude and
Figure 6: Pearson correlations between environmental variables, measures of zeta
diversity, mean CSCI scores, and modeled indices using diatom communities classified
using the 18S-V9 primer and aggregated to various taxonomic levels (p < 10
-4
).
Figure 7: Pearson correlations between environmental variables, measures of zeta
diversity, mean CSCI scores, and modeled indices using soft algal communities classified
using the 18S-V9 primer and aggregated to various taxonomic levels (p < 10
-4
).
114
local diversity (Figures 6A-D). Landscape diversity was found to increase with altitude in
communities of prokaryotes classified using the 16S-V4 primer and aggregated to the level of
species or order (Figures 8A, D), prokaryotes classified using the 16S 515F-926R primer and
aggregated to order (Figure 5D), BMIs aggregated to order (Figure 4D), and with communities
of diatoms classified using the rbcL primer (Figures 3A-D). Measures of landscape diversity
were found to be negatively correlated with altitude in communities of soft algae (Figures 7A-
D), and diatoms classified using the 18S-V9 primer and aggregated to the level of genus or
higher (Figures 6B-D).

𝞯 diversity and distance
For our
measures of
diversity, we
observed the
following four
trends with
relation to
distance. An
increase in
distance was
found to be associated with an increase in local diversity for communities of BMIs (Figures 4A-
D), prokaryotes classified using the 16S 515F-926R primer and aggregated above the level of
genus (Figures 5C-D), and diatoms classified using the rbcL primer (Figures 3A-D). The
Figure 6: Pearson correlations between environmental variables, measures of zeta
diversity, mean CSCI scores, and modeled indices using prokaryotic communities
classified using the 16S-V4 primer and aggregated to various taxonomic levels (p < 10
-4
).
115
negative trend between local diversity and distance was observed for diatoms classified by the
18S-V9 primer, but only for communities aggregated above the level of genus (Figures 6C-D).
No consistent trends between our measure of local diversity and distance were observed for
prokaryotes classified using the 16S-V4 primer (Figures 8A-D) or soft algae (Figures 7A-D). For
measures of landscape diversity we found positive correlations with measures of landscape
diversity for diatoms classified using the rbcL primer (Figures 3A-D), prokaryotes classified
using the 16S-V4 primer and aggregated at the level of species or order (Figures 8A, D),
prokaryotes classified using the 16S 515F-926R primer and aggregated to order (Figure 5D) and
BMIs (Figures 4A-D). A negative trend was observed between distance and measures of
landscape diversity for communities of soft algae (Figures 7A-D), and diatoms classified using
the 18S-V9 primer and aggregated above the level of genus (Figures 6C-D).
𝞯 diversity and models of community assembly
For prokaryotes, as well as diatoms class using the rbcL primer, we found a power-law
model was likelier than an exponential one in describing the decay of zeta diversity with sample
number across all of the taxonomic levels considered (Table 3). These decay patterns in zeta
diversity indicate communities composed of organisms classified using 16S-V4, 16S 515F-926R,
and rbcL primers are likely following a niche differentiation model of community assembly
across taxonomic levels of aggregation ranging from species to order. With organisms classified
using the 18S-V9 primer we found a greater likelihood for a niche differentiated model of
community assembly above the taxonomic level of genus for BMIs, above species for diatoms,
and above family for soft algae. At lower taxonomic levels of aggregation, the communities
described using the 18S-V9 primer appeared to follow a stochastic model of community
assembly.
116
The relative importance of diversity measures in our modeled indices
Across all taxonomic levels considered we found 𝞯1, our mean measure of local diversity,
for both BMIs and diatoms classified using the rbcL primer to be a relatively important
component in our models of the mean CSCI score (Table 4). For prokaryotes, soft algae, and
diatoms our measures of landscape diversity, 𝞯2 and 𝞯10, tended to be of greater relative
importance in modeling the mean CSCI score.
Table 4: The relative importance of zeta diversity parameters in describing linear models of the mean value of the
CSCI with communities aggregated to various taxonomic levels.
18S-V9
Soft algae
18S-V9
Diatoms
18S-V9
BMIs
16S-V4
Prokaryotes
16S
515F-926R
Prokaryotes
rbcL
Diatoms
Species 𝞯 1: 4.7
𝞯 2: 67.2
𝞯 10: 3.5
𝞯 1: 19.1
𝞯 2: 23.7
𝞯 10: 8.3
𝞯 1: 44.7
𝞯 2: 3.8
𝞯1 0: 1.3
𝞯 1: 2.0
𝞯 2: 34.8
𝞯1 0: 25.0
𝞯 1: 0.1
𝞯 2: 45.5
𝞯 10: 2.1
𝞯 1: 32.5
𝞯 2: 31.2
𝞯 10: 20.0
Genus 𝞯 1: 2.4
𝞯 2: 66.3
𝞯 10: 13.4
𝞯 1: 5.7
𝞯 2: 4.3
𝞯1 0: 12.6
𝞯 1: 62.7
𝞯 2: 9.7
𝞯 10: 1.6
𝞯 1: 0.2
𝞯 2: 1.9
𝞯 10: 6.7
𝞯 1: 6.8
𝞯 2: 2.1
𝞯 10: 2.7
𝞯 1: 27.5
𝞯 2: 37.7
𝞯 10: 30.8
Family 𝞯 1: 5.5
𝞯 2: 64.8
𝞯 10: 12.7
𝞯 1: 7.6
𝞯 2: 17.0
𝞯1 0: 30.0
𝞯 1: 33.5
𝞯 2: 7.3
𝞯1 0: 24.1
𝞯 1: 12.7
𝞯 2: 10.4
𝞯 10: 12.5
𝞯 1: 30.9
𝞯 2: 3.1
𝞯 10: 4.2
𝞯 1: 31.1
𝞯 2: 33.9
𝞯 10: 30.4
order 𝞯 1: 1.4
𝞯 2: 86.9
𝞯 1: 17.0
𝞯 2: 34.8
𝞯 1: 12.9
𝞯 2: 32.3
𝞯 1: 7.0
𝞯 2: 17.5
𝞯 1: 8.3
𝞯 2: 28.6
𝞯 1: 28.3
𝞯 2: 33.5
117
𝞯 10: 4.8 𝞯 10: 10.0 𝞯 10: 15.5 𝞯 10: 40.3 𝞯 10: 35.2 𝞯 10: 34.7

The relative importance of environmental measures in our modeled indices
We observed similar patterns in relative importance (land use > altitude > distance) of
our measures of the physical environment in both the mean CSCI score and our modeled indices
for communities of soft algae, BMIs, and diatoms classified using the rbcL primer (Table 5). For
communities of prokaryotes, and diatoms classified using the 18S-V9 primer, the relative
importance of environmental variables were both diminished and inconsistent across taxonomic
levels between the modeled and mean CSCI values. However, in all communities we tended to
see land use as being the most important of the three environmental variables considered.

118

Table 5: The relative importance of land use, altitude, and distance in describing linear models for modeled and
(mean) CSCI values using communities aggregated to various taxonomic levels.
18S-V9
Soft algae
18S-V9
Diatoms
18S-V9
BMIs
16S-V4
Prokaryotes
16S
515F-926R
Prokaryotes
rbcL
Diatoms
Species Land use:
17.5 (26.4)
Altitude:
37.1 (45.3)
Distance:
23.1 (27.1)
Land use:
31.3 (47.9)
Altitude:
10.6 (26.4)
Distance:
12.9 (24.4)
Land use:
34.3 (48.2)
Altitude:
12.5 (27.0)
Distance:
19.4 (23.2)
Land use:
36.5 (49.4)
Altitude:
16.8 (26.9)
Distance:
11.2 (22.4)
Land use:
17.0 (42.3)
Altitude:
22.4 (25.7)
Distance:
17.8 (29.3)
Land use:
70.9 (62.6)
Altitude:
16.5 (25.3)
Distance:
5.7 (11.4)
Genus Land use:
40.2 (47.0)
Altitude:
22.5 (26.3)
Distance:
20.2 (25.4)
Land use:
10.1 (48.2)
Altitude:
15.9 (27.0)
Distance:
3.3 (23.2)
Land use:
40.5 (47.3)
Altitude:
15.9 (26.4)
Distance:
21.9 (25.1)
Land use:
10.4 (42.4)
Altitude:
10.7 (26.3)
Distance:
11.8 (28.6)
Land use:
5.2 (42.3)
Altitude:
13.8 (25.7)
Distance:
3.0 (29.3)
Land use:
61.6 (62.6)
Altitude:
23.7 (25.3)
Distance:
11.0 (11.4)
family Land use:
40.0 (46.9)
Altitude:
22.8 (25.4)
Land use:
22.2 (46.2)
Altitude:
20.5 (26.1)
Land use:
31.0 (46.2)
Altitude:
15.5 (25.6)
Land use:
14.2 (42.4)
Altitude:
12.8 (26.3)
Land use:
14.5 (42.3)
Altitude:
14.4 (25.7)
Land use:
58.7 (62.6)
Altitude:
26.1 (25.3)
119
Distance:
20.3 (26.4)
Distance:
13.0 (26.5)
Distance:
20.6 (26.9)
Distance:
8.7 (28.6)
Distance:
12.1 (29.3)
Distance:
11.3 (11.4)
Order Land use:
27.1 (27.5)
Altitude:
43.4 (46.3)
Distance:
23.4 (24.9)
Land use:
24.4 (46.9)
Altitude:
31.1 (26.2)
Distance:
14.0 (25.6)
Land use:
35.0 (46.6)
Altitude:
13.2 (26.6)
Distance:
21.8 (25.5)
Land use:
35.0 (49.4)
Altitude:
12.7 (26.9)
Distance:
21.4 (22.4)
Land use:
41.6 (48.0)
Altitude:
14.5 (26.4)
Distance:
22.2 (24.4)
Land use:
64.9 (62.6)
Altitude:
23.0 (25.3)
Distance:
9.8 (11.4)

Discussion
We found zeta diversity to play a potential role in the investigation of likely models of
stream community assembly and, at least for some types of biological communities, a viable tool
for the development of new bioassessments. In applying zeta diversity analysis to stream
communities we found it to have potential in both assessing the health of stream communities as
well as investigating the likelihood of models for community assembly for prokaryotes, soft
algae, diatoms, and BMIs. In comparing the likelihoods of models for how zeta diversity decays
with order we were able to determine the most probable model of assembly for a variety of
communities in stream. In this analysis we found communities of prokaryotes and diatoms
classified using the rbcL primer were most likely assembled following a niche differentiation
process. For soft algae, BMIs, and diatoms classified using the 18S-V9 primer the likelihood of
either model of community assembly depended on the level of taxonomic aggregation, with the
120
general trend of a niche differentiated model becoming likelier at higher levels of taxonomic
aggregation. While lower orders of zeta diversity, 𝞯1 and 𝞯2, were often found to explain sizable
fractions of the observed variation on the mean CSCI score, our measure of regional stream
community health, we generally found 𝞯10 to be of significant importance as well in our models.
We used zeta diversity, which is a generalized framework for analyzing presence/absence
patterns for unique categories of organisms across communities, to construct modeled indices of
the health of stream communities. As these modeled indices based on zeta diversity are not
reliant on the use of undisturbed reference communities, or the presence of particular indicator
species, we believe they have potential use in constructing tools to assess the health of stream
communities beyond the original geographic scope of the data used in constructing the CSCI.
Metabarcoding and limitations with taxonomic assignments
We found sizable differences between community types in terms of the number of unique
groups which could be resolved per taxonomic level (Figure 2B). This is likely due to differences
in the ability of different metabarcoding primers to resolve members of different groups of
organisms to particular taxonomic levels (Clarke et al., 2014; Debroas et al., 2017; Elbrecht and
Leese, 2017), as well as gaps in various reference databases used to assign taxonomies
(Somervuo et al., 2017; von Ammon et al., 2018). For example, the larger proportion of resolved
taxonomic groups in diatom communities classified using the rbcL primer versus ones targeted
ribosomal sequences most likely reflects both mechanisms (Kermarrec et al., 2013).
Completeness of reference databases may also reflect biases towards well relatively studied
groups, such as freshwater BMIs (Elbrecht et al., 2017). While particular groups of prokaryotes
may be well studied we found significant limitations in our taxonomic resolution for a large
portion of our 16S amplicons. This may simply reflect limitations in being able to distinguish
121
prokaryotes based solely on the identification of individual hypervariable regions (Johnson et al.,
2019).
The importance of local and landscape diversity in modeling regional
measures of stream community health
While we found our measures of landscape diversity, 𝞯2 and 𝞯10, frequently played
significant roles in our models of biotic integrity, they are potentially capturing different patterns
across the communities being sampled. As zeta diversity order increases, for example from 𝞯2 to
𝞯10, we find an increased sensitivity to the presence of common categories of organisms
(Latombe et al. 2017a). To varying degrees, across all of the communities analyzed, we observed
a significant role for both our measures of landscape diversity in our regional models of stream
health. As a result, we propose that assessments of the impacts of environmental degradation
cannot be described solely from the loss of relatively rare or common categories of organisms,
but that such changes must be accounted for together. For a number of communities, in particular
those where most of the ASVs have taxonomic assignments, we found a sizable portion of the
observed variation in mean CSCI values could be attributed to our modeled indices (Figures 3-
8). These indices are simple linear models of only three orders of zeta diversity: 𝞯1, 𝞯2, and 𝞯10.
This simplicity suggests that future assessments of the health of streams and a regional scale
could be accomplished with a more generalized framework than the one used to construct the
CSCI.
122
Measures of regional stream health and environmental variables
Of the environmental measures considered, we tended to observe land use as having the
greatest relative importance for our modeled indices (Table 5). We also tended to observe
opposite trends comparing modeled indices to either altitude or distance as those compared to
land use (Figures 3-8). While we used calculations of relative importance to help account for
collinearity in our environmental variables, we do note broadly negative trends between land use
and altitude (r = -0.86, p < 10
-4
) as well as land use and distance (r = -0.85, p < 10
-4
) across our
sample groups.
For BMI communities we found, at least below the taxonomic level of order, a high
relative importance for 𝞯1 and our modeled indices (Table 4). The value of 𝞯1 is equivalent to the
average local diversity for a set of samples. This then corresponds to expected behavior with the
CSCI, which is a composite measure of both taxonomic and functional feeding group local
diversity for BMIs (Mazor et al., 2016). We observed a negative correlation between modeled
indices for BMI communities and land use (Figures 4A-D), reflecting prior observations of a
decline in local BMI diversity associated with an increase in environmental stress in general
(Paul and Meyer 2001, Stepenuck et al. 2002, Ourso and Frenzel 2003, Wallace and Biastoch
2016) and land use in particular (Sponseller et al. 2001, Stepenuck et al. 2002, Allan 2004;
Fugère et al., 2016). Reflecting prior observations between the health of BMI communities and
altitude (Ward, 1986; Jacobsen et al. 1997, Jacobsen 2008), which tends to decline with land use,
we found a positive correlation between altitude and our modeled indices for BMI communities.
For BMI communities, aggregated above the level of genus, we observed significant negative
relationships between our measures of landscape diversity and land use (Figures 4C-D). These
findings are in accordance with prior observations which suggest a greater dissimilarity between
123
the composition of BMI communities and levels of environmental disturbance across a landscape
(Gutiérrez‐Cánovas et al., 2013; Hawkins et al., 2015; Simons et al., 2019). These results lend
support to the utility of building assessments of the health of regional groupings of BMI
communities using zeta diversity and metabarcoding.
We found that the modeled indices which performed best, in accounting for the observed
variation in the mean CSCI, were those based on measures of zeta diversity for communities of
diatoms classified using the rbcL primer (Table 2, Figures 3A-D). This may reflect the
approximately equal relative importance of measures of both local and landscape diversity to
modeled indices constructed using these diatom community data sets (Table 4), which implies a
sensitivity to changes in both widespread and localized taxonomic groups. The high dispersal
ability of diatoms, compared to other groups such as soft algae and BMIs (Padial et al., 2014),
may in part underpin the sensitivity of changes to their metacommunity structure in response to
environmental factors. We also find prior evidence of highly similar changes in patterns of
diversity for diatoms and BMIs, at both a local and regional scale, associated with environmental
disturbances (Passy et al., 2004; Johnson et al., 2007; Johnson and Hering, 2010; Vázquez et al.,
2011). Although these similarities may simply reflect that our accepted measure of stream health,
based on changes to the composition of BMI communities, may be expected to vary across
environmental gradients in a very similar way to the composition of diatom communities.
For soft algal communities we tended to observe a convergence in community
composition associated with a rise in land use (Figures 7A-D). This potential relationship
between diversity and disturbance is in agreement prior observations of streambed soft algal
communities (Passy and Blanchet, 2007), which in turn may simply reflect strong environmental
filtering in soft algal communities due to anthropogenic disturbances (Dunck et al., 2019).
124
Across the taxonomic levels considered we observed no significant trend between
environmental disturbance, as described by land use, and local diversity of prokaryotes classified
by the 16S-V4 primer (Figures 8A-D). A similar lack of correlation has been observed in prior
studies of prokaryotic communities in streams along disturbance gradients (Lear et al., 2011;
Simonin et al., 2019). While local diversity appears to be unaffected by an increase in land use,
for measures of landscape diversity we found some evidence that similarities in the composition
of prokaryotic communities decline with it (Figures 5A, 5D, 6A, 6D). This again reflects prior
observations of prokaryotic communities in streams subjected to anthropogenic stress (Lear et
al., 2011; Hosen et al., 2017; Laperriere et al., 2020).
For prokaryotes classified using the 16S 515F-926R primer we tended to observe a
negative trend between environmental disturbance, as described by land use, and local
prokaryotic diversity, as described by 𝞯1 (Figures 5B-D). However, no such consistent trend
could be found with prokaryotic communities classified using the 16S-V4 primer (Figures 8A-
D). This may primarily be due to the greater length of the region amplified using the 16S 515F-
926R (381 bp) rather than the 16S-V4 primer (291 bp), which has been found to better capture
the diversity of prokaryotic communities (Parada et al., 2016; Yang et al., 2016). While a lack of
correlation has been observed in prior studies of prokaryotic communities in streams along
disturbance gradients (Lear et al., 2011; Simonin et al., 2019), more recent analysis indicates a
negative relationship between land use and local prokaryotic diversity (Laperriere et al., 2020).
For measures of landscape diversity, we tended to see similarities in the composition of
prokaryotic communities decline with it, although this pattern was far clearer when prokaryotes
were classified using the 16S 515F-926R rather than 16S-V4 primer. This again reflects prior
125
observations of prokaryotic communities in streams subjected to anthropogenic stress (Lear et
al., 2011; Hosen et al., 2017; Laperriere et al., 2020).
Of the environmental measures considered, we tended to observe land use as having the
greatest relative importance for our modeled indices (Table 5). We also tended to observe
opposite trends comparing modeled indices to either altitude or distance as those compared to
land use (Figures 3-8). While we used calculations of relative importance to help account for
collinearity in our environmental variables, we do note broadly negative trends between land use
and altitude (r = -0.86, p < 10
-4
) as well as land use and distance (r = -0.85, p < 10
-4
) across our
sample groups.
For BMI communities we found, at least below the taxonomic level of order, a high
relative importance for 𝞯1 and our modeled indices (Table 4). The value of 𝞯1 is equivalent to the
average local diversity for a set of samples. This then corresponds to expected behavior with the
CSCI, which is a composite measure of both taxonomic and functional feeding group local
diversity for BMIs (Mazor et al., 2016). We observed a negative correlation between modeled
indices for BMI communities and land use (Figures 4A-D), reflecting prior observations of a
decline in local BMI diversity associated with an increase in environmental stress in general
(Paul and Meyer 2001, Stepenuck et al. 2002, Ourso and Frenzel 2003, Wallace and Biastoch
2016) and land use in particular (Sponseller et al. 2001, Stepenuck et al. 2002, Allan 2004;
Fugère et al., 2016). Reflecting prior observations between the health of BMI communities and
altitude (Ward, 1986; Jacobsen et al. 1997, Jacobsen 2008), which tends to decline with land use,
we found a positive correlation between altitude and our modeled indices for BMI communities.
For BMI communities, aggregated above the level of genus, we observed significant negative
relationships between our measures of landscape diversity and land use (Figures 4C-D). These
126
findings are in accordance with prior observations which suggest a greater dissimilarity between
the composition of BMI communities and levels of environmental disturbance across a landscape
(Gutiérrez‐Cánovas et al., 2013; Hawkins et al., 2015; Simons et al., 2019). These results lend
support to the utility of building assessments of the health of regional groupings of BMI
communities using zeta diversity and metabarcoding.
We found that the modeled indices which performed best, in accounting for the observed
variation in the mean CSCI, were those based on measures of zeta diversity for communities of
diatoms classified using the rbcL primer (Table 2, Figures 3A-D). This may reflect the
approximately equal relative importance of measures of both local and landscape diversity to
modeled indices constructed using these diatom community data sets (Table 4), which implies a
sensitivity to changes in both widespread and localized taxonomic groups. The high dispersal
ability of diatoms, compared to other groups such as soft algae and BMIs (Padial et al., 2014),
may in part underpin the sensitivity of changes to their metacommunity structure in response to
environmental factors. We also find prior evidence of highly similar changes in patterns of
diversity for diatoms and BMIs, at both a local and regional scale, associated with environmental
disturbances (Passy et al., 2004; Johnson et al., 2007; Johnson and Hering, 2010; Vázquez et al.,
2011). Although these similarities may simply reflect that our accepted measure of stream health,
based on changes to the composition of BMI communities, may be expected to vary across
environmental gradients in a very similar way to the composition of diatom communities.
For soft algal communities we tended to observe a convergence in community
composition associated with a rise in land use (Figures 7A-D). This potential relationship
between diversity and disturbance is in agreement prior observations of streambed soft algal
127
communities (Passy and Blanchet, 2007), which in turn may simply reflect strong environmental
filtering in soft algal communities due to anthropogenic disturbances (Dunck et al., 2019).
Across the taxonomic levels considered we observed no significant trend between
environmental disturbance, as described by land use, and local diversity of prokaryotes classified
by the 16S-V4 primer (Figures 8A-D). A similar lack of correlation has been observed in prior
studies of prokaryotic communities in streams along disturbance gradients (Lear et al., 2011;
Simonin et al., 2019). While local diversity appears to be unaffected by an increase in land use,
for measures of landscape diversity we found some evidence that similarities in the composition
of prokaryotic communities decline with it (Figures 5A, 5D, 6A, 6D). This again reflects prior
observations of prokaryotic communities in streams subjected to anthropogenic stress (Lear et
al., 2011; Hosen et al., 2017; Laperriere et al., 2020).
For prokaryotes classified using the 16S 515F-926R primer we tended to observe a
negative trend between environmental disturbance, as described by land use, and local
prokaryotic diversity, as described by 𝞯1 (Figures 5B-D). However, no such consistent trend
could be found with prokaryotic communities classified using the 16S-V4 primer (Figures 8A-
D). This may primarily be due to the greater length of the region amplified using the 16S 515F-
926R (381 bp) rather than the 16S-V4 primer (291 bp), which has been found to better capture
the diversity of prokaryotic communities (Parada et al., 2016; Yang et al., 2016). While a lack of
correlation has been observed in prior studies of prokaryotic communities in streams along
disturbance gradients (Lear et al., 2011; Simonin et al., 2019), more recent analysis indicates a
negative relationship between land use and local prokaryotic diversity (Laperriere et al., 2020).
For measures of landscape diversity, we tended to see similarities in the composition of
prokaryotic communities decline with it, although this pattern was far clearer when prokaryotes
128
were classified using the 16S 515F-926R rather than 16S-V4 primer. This again reflects prior
observations of prokaryotic communities in streams subjected to anthropogenic stress (Lear et
al., 2011; Hosen et al., 2017; Laperriere et al., 2020).
𝞯 diversity decay and models of community assembly
In observing how zeta diversity decays with respect to zeta order we found a power-law
decay model was more likely than an exponential one for communities of prokaryotes and
diatoms classified using the rbcL primer (Table 3). This evidence supports a model of
community assembly for these prokaryotes and diatoms, aggregated from the levels of species
through order, determined more by niche differentiation than a stochastic process (Scheiner et al.
2011, Hui and McGeoch 2014, Hui et al. 2018). These results support prior observations of niche
differentiation being the driving force for community assembly in stream communities of both
prokaryotes (Van der Gucht et al., 2007) and diatoms (Pan et al., 1999; Soininen 2007). For all
three community types classified using the 18S-V9 primer we also found evidence of a
community assembly process determined by niche differentiation, but only at taxonomic levels at
or above order, genus, and family for soft algae, diatoms, and BMIs respectively. At lower
taxonomic levels a stochastic assembly process appears more likely for the communities
classified using a 18S-V9 primer.
In general, we found the power-law decay model for how zeta diversity decays with zeta
order becomes increasingly likely compared to an exponential one as we move from
communities aggregated at the level of species to those aggregated at the level of order (Table 3).
Prior analyses support the hypothesis of an increasingly strong, albeit broad, relationship
between functional and taxonomic diversity at higher taxonomic levels (Anacker and Strauss,
2016; Keck et al., 2016; Lu et al., 2016). Our data then suggests a greater likelihood for
129
communities to assemble via a process of niche differentiation process when studied at higher
taxonomic levels corresponding more strongly to functional groupings. This is in line with
observations of stronger relationships between functional diversity and the biotic filters driving
niche differentiation as a process of community assembly than with taxonomic diversity (Siefert
et al., 2013; Khalil et al., 2018).
Conclusion
We found diversity patterns of diatom communities to be the most responsive to
environmental gradients amongst the stream communities we studied, forming the basis for our
most accurate indices of stream health. Although it should be noted that we could only obtain
previously observed trends relating the diversity of diatom communities and environmental
gradients (Passy et al., 2004; Johnson et al., 2007; Johnson and Hering, 2010; Vázquez et al.,
2011), as well as significantly higher taxonomic resolution, when communities were classified
using a rbcL rather than 18S-V9 primer.
We have observed a number of significant relationships between measures of landscape
diversity, as described by 𝞯2 and 𝞯10, and geospatial variations in environmental conditions. Our
results tend to indicate negative correlations between the average local diversity of regional
communities of BMIs, diatoms classified using the rbcL primer, and prokaryotes and land use,
while the opposite trend was found for communities of soft algae and diatoms classified using
the 18S-V9 primer. With relation to altitude we observed a positive correlation between the
average local diversity of all community types, except for diatoms classified using the 18S-V9
primer. These trends are in general agreement with prior observations regarding environmental
gradients and the local diversity of communities of BMIs (Ward 1986, Jacobsen et al. 1997,
130
Sponseller et al. 2001, Huryn et al. 2002, Stepenuck et al. 2002, Allan 2004, Finn and Leroy Poff
2005, Jacobsen 2008; Fugère et al., 2016), soft algae (Passy and Blanchet, 2007; Vázquez et al.,
2011; Dunck et al., 2019), and prokaryotes (Laperriere et al., 2020).
We have demonstrated the potential utility of applying analysis of zeta diversity patterns
of various types of stream communities as a means of assessing their ecological health. Although
the accuracy of this approach appears to be strongly dependent on the completeness of the
reconstructed structure of communities. Which are in turn dependent on the community type,
primer, and the completeness of the reference library used to assign taxonomies (Marcy et al.,
2007; Baker et al., 2010; Lan et al., 2012; Mande et al., 2012; Rinke et al., 2013; Christensen and
Olsen, 2018). To overcome these limitations, future work on the development of biological
assessments of streams using zeta diversity and sequence data could move towards reference-free
based approaches of classifying sequences and organizing them into communities (Dubinika et
al., 2016; Tang et al., 2019, Ren et al. 2018, Zielezinski et al. 2019, Linard et al., 2019). Such
approaches have shown a significant improvement over standard taxonomic approaches in the
completeness with which sequences can be incorporated into profiles of community composition
(Apothéloz-Perret-Gentil et al., 2017), and have begun to demonstrate promise in building tools
for biological assessment (Pawlowski et al., 2018). Future work on the development of biological
assessments, based on patterns of zeta diversity, are likely to then improve with the incorporation
of community data constructed using reference-free approaches to the classification of
sequences.

131
References
Allan, J. D. (2004). Landscapes and riverscapes: the influence of land use on stream ecosystems.
Annu. Rev. Ecol. Evol. Syst., 35, 257-284.

Amaral-Zettler, L. A., McCliment, E. A., Ducklow, H. W., & Huse, S. M. (2009). A method for
studying protistan diversity using massively parallel sequencing of V9 hypervariable regions of
small-subunit ribosomal RNA genes. PloS one, 4(7).

Anacker, B. L., & Strauss, S. Y. (2016). Ecological similarity is related to phylogenetic distance
between species in a cross‐niche field transplant experiment. Ecology, 97(7), 1807-1818.

Apothéloz‐Perret‐Gentil, L., Cordonier, A., Straub, F., Iseli, J., Esling, P., & Pawlowski, J.
(2017). Taxonomy‐free molecular diatom index for high‐throughput eDNA biomonitoring.
Molecular ecology resources, 17(6), 1231-1242.

Apprill, A., McNally, S., Parsons, R., & Weber, L. (2015). Minor revision to V4 region SSU
rRNA 806R gene primer greatly increases detection of SAR11 bacterioplankton. Aquatic
Microbial Ecology, 75(2), 129-137.

Bailet, B., Bouchez, A., Franc, A., Frigerio, J. M., Keck, F., Karjalainen, S. M., ... & Kahlert, M.
(2019). Molecular versus morphological data for benthic diatoms biomonitoring in Northern
Europe freshwater and consequences for ecological status.

132
Baird, D. J., & Sweeney, B. W. (2011). Applying DNA barcoding in benthology: the state of the
science. Journal of the North American Benthological Society, 30(1), 122-124.

Baker, B. J., Comolli, L. R., Dick, G. J., Hauser, L. J., Hyatt, D., Dill, B. D., ... & Banfield, J. F.
(2010). Enigmatic, ultrasmall, uncultivated Archaea. Proceedings of the National Academy of
Sciences, 107(19), 8806-8811.

Baker, David J., Stephen T. Garnett, James O'Connor, Glenn Ehmke, Rohan H. Clarke, John CZ
Woinarski, and Melodie A. McGeoch. "Conserving the abundance of nonthreatened species."
Conservation Biology 33, no. 2 (2019): 319-328.

Bandh, S. A., Shafi, S., Shameem, N., Dar, R., Kamili, A. N., & Ganai, B. A. (2019). Spatio-
temporal patterns of bacterial diversity along environmental gradients and bacterial attachment to
organic aggregates. In Freshwater Microbiology (pp. 137-174). Academic Press.

Banerji, A., Bagley, M., Elk, M., Pilgrim, E., Martinson, J., & Santo Domingo, J. (2018). Spatial
and temporal dynamics of a freshwater eukaryotic plankton community revealed via 18S rRNA
gene metabarcoding. Hydrobiologia, 818(1), 71-86.

Battin, T. J., Besemer, K., Bengtsson, M. M., Romani, A. M., & Packmann, A. I. (2016). The
ecology and biogeochemistry of stream biofilms. Nature Reviews Microbiology, 14(4), 251.

133
Beck, Jan, Jeremy D. Holloway, and Wolfgang Schwanghart. "Undersampling and the
measurement of beta diversity." Methods in Ecology and Evolution 4, no. 4 (2013): 370-382.

Besemer, K. (2015). Biodiversity, community structure and function of biofilms in stream
ecosystems. Research in microbiology, 166(10), 774-781.

Bolyen, E., Rideout, J. R., Dillon, M. R., Bokulich, N. A., Abnet, C. C., Al-Ghalith, G. A., ... &
Bai, Y. (2019). Reproducible, interactive, scalable and extensible microbiome data science using
QIIME 2. Nature biotechnology, 37(8), 852-857.

Britton, A. W., J. J. Day, C. J. Doble, B. P. Ngatunga, K. M. Kemp, C. Carbone, and D. J.
Murrell. "Terrestrial-focused protected areas are effective for conservation of freshwater fish
diversity in Lake Tanganyika." Biological Conservation 212 (2017): 120-129.

Brown, E. A., Chain, F. J., Crease, T. J., MacIsaac, H. J., & Cristescu, M. E. (2015). Divergence
thresholds and divergent biodiversity estimates: can metabarcoding reliably describe
zooplankton communities?. Ecology and Evolution, 5(11), 2234-2251.

Callahan, B. J., McMurdie, P. J., Rosen, M. J., Han, A. W., Johnson, A. J. A., & Holmes, S. P.
(2016). DADA2: high-resolution sample inference from Illumina amplicon data. Nature
methods, 13(7), 581.

134
Carew, M. E., Pettigrove, V. J., Metzeling, L., & Hoffmann, A. A. (2013). Environmental
monitoring using next generation sequencing: rapid identification of macroinvertebrate
bioindicator species. Frontiers in zoology, 10(1), 45.

Castilla, E. P., Cunha, D. G. F., Lee, F. W. F., Loiselle, S., Ho, K. C., & Hall, C. (2015).
Quantification of phytoplankton bloom dynamics by citizen scientists in urban and peri-urban
environments. Environmental monitoring and assessment, 187(11), 690.

Chao, Anne, Lou Jost, S. C. Chiang, Y‐H. Jiang, and Robin L. Chazdon. "A two‐stage
probabilistic approach to multiple‐community similarity indices." Biometrics 64, no. 4 (2008):
1178-1186.

Chen, W., Wilkes, G., Khan, I. U., Pintar, K. D., Thomas, J. L., Lévesque, C. A., ... & Lapen, D.
R. (2018). Aquatic bacterial communities associated with land use and environmental factors in
agricultural landscapes using a metabarcoding approach. Frontiers in microbiology, 9, 2301.

Christensen, H., & Olsen, J. E. (2018). Full Shotgun DNA Metagenomics. In Introduction to
Bioinformatics in Microbiology (pp. 163-175). Springer, Cham.

Clarke, L. J., Soubrier, J., Weyrich, L. S., & Cooper, A. (2014). Environmental metabarcodes for
insects: in silico PCR reveals potential for taxonomic bias. Molecular ecology resources, 14(6),
1160-1170.

135
Colvin, S. A., Sullivan, S. M. P., Shirey, P. D., Colvin, R. W., Winemiller, K. O., Hughes, R. M.,
... & Danehy, R. J. (2019). Headwater streams and wetlands are critical for sustaining fish,
fisheries, and ecosystem services. Fisheries, 44(2), 73-91.

Cordier, T., Lanzén, A., Apothéloz-Perret-Gentil, L., Stoeck, T., & Pawlowski, J. (2019).
Embracing environmental genomics and machine learning for routine biomonitoring. Trends in
microbiology, 27(5), 387-397.

Costa, S. S., & Melo, A. S. (2008). Beta diversity in stream macroinvertebrate assemblages:
among-site and among-microhabitat components. Hydrobiologia, 598(1), 131-138.

Debroas, D., Domaizon, I., Humbert, J. F., Jardillier, L., Lepère, C., Oudart, A., & Taïb, N.
(2017). Overview of freshwater microbial eukaryotes diversity: a first analysis of publicly
available metabarcoding data. FEMS Microbiology Ecology, 93(4), fix023.

Dubinkina, V. B., Ischenko, D. S., Ulyantsev, V. I., Tyakht, A. V., & Alexeev, D. G. (2016).
Assessment of k-mer spectrum applicability for metagenomic dissimilarity analysis. BMC
bioinformatics, 17(1), 38.

Dunck, B., Felisberto, S. A., & de Souza Nogueira, I. (2019). Effects of freshwater
eutrophication on species and functional beta diversity of periphytic algae. Hydrobiologia,
837(1), 195-204.

136
Elbrecht, V., & Leese, F. (2016). Validation and development of freshwater invertebrate
metabarcoding COI primers for Environmental Impact Assessment (No. e2044v3). PeerJ
Preprints.

Elbrecht, V., & Leese, F. (2017). Validation and development of COI metabarcoding primers for
freshwater macroinvertebrate bioassessment. Frontiers in Environmental Science, 5, 11.

Finn, D. S., & Leroy Poff, N. (2005). Variability and convergence in benthic communities along
the longitudinal gradients of four physically similar Rocky Mountain streams. Freshwater
Biology, 50(2), 243-261.

Fugère, V., Kasangaki, A., & Chapman, L. J. (2016). Land use changes in an afrotropical
biodiversity hotspot affect stream alpha and beta diversity. Ecosphere, 7(6), e01355.

Gál, B., Szivák, I., Heino, J., & Schmera, D. (2019). The effect of urbanization on freshwater
macroinvertebrates–knowledge gaps and future research directions. Ecological indicators, 104,
357-364.

Gaston, Kevin J., and Richard A. Fuller. "Commonness, population depletion and conservation
biology." Trends in ecology & evolution 23, no. 1 (2008): 14-19.

137
Gerth, W. J., & Herlihy, A. T. (2006). Effect of sampling different habitat types in regional
macroinvertebrate bioassessment surveys. Journal of the North American Benthological Society,
25(2), 501-512.

Gibson, J. F., Shokralla, S., Curry, C., Baird, D. J., Monk, W. A., King, I., & Hajibabaei, M.
(2015). Large-scale biomonitoring of remote and threatened ecosystems via high-throughput
sequencing. PloS one, 10(10).

Gilbert, J. A., Jansson, J. K., & Knight, R. (2014). The Earth Microbiome project: successes and
aspirations. BMC biology, 12(1), 69.

Grömping, U. (2006). Relative importance for linear regression in R: the package relaimpo.
Journal of statistical software, 17(1), 1-27.

Grönroos, Mira, Jani Heino, Tadeu Siqueira, Victor L. Landeiro, Juho Kotanen, and Luis M.
Bini. "Metacommunity structuring in stream networks: roles of dispersal mode, distance type,
and regional environmental context." Ecology and evolution 3, no. 13 (2013): 4473-4487.

Guo, L., Sui, Z., & Liu, Y. (2016). Quantitative analysis of dinoflagellates and diatoms
community via Miseq sequencing of actin gene and v9 region of 18S rDNA. Scientific reports, 6,
34709.

138
Gutiérrez‐Cánovas, C., Millán, A., Velasco, J., Vaughan, I. P., & Ormerod, S. J. (2013).
Contrasting effects of natural and anthropogenic stressors on beta diversity in river organisms.
Global Ecology and Biogeography, 22(7), 796-805.

Hajibabaei, M., Shokralla, S., Zhou, X., Singer, G. A., & Baird, D. J. (2011). Environmental
barcoding: a next-generation sequencing approach for biomonitoring applications using river
benthos. PloS one, 6(4).

Hawkins, C. P., Olson, J. R., & Hill, R. A. (2010). The reference condition: predicting
benchmarks for ecological and water-quality assessments. Journal of the North American
Benthological Society, 29(1), 312-343.

Hawkins, C. P., Mykrä, H., Oksanen, J., & Vander Laan, J. J. (2015). Environmental disturbance
can increase beta diversity of stream macroinvertebrate assemblages. Global Ecology and
Biogeography, 24(4), 483-494.

Hebert, P. D., Penton, E. H., Burns, J. M., Janzen, D. H., & Hallwachs, W. (2004). Ten species
in one: DNA barcoding reveals cryptic species in the neotropical skipper butterfly Astraptes
fulgerator. Proceedings of the National Academy of Sciences, 101(41), 14812-14817.

Heino, Jani. "The importance of metacommunity ecology for environmental assessment research
in the freshwater realm." Biological Reviews 88, no. 1 (2013): 166-178.

139
Herbst, D. B., & Silldorff, E. L. (2006). Comparison of the performance of different
bioassessment methods: similar evaluations of biotic integrity from separate programs and
procedures. Journal of the North American Benthological Society, 25(2), 513-530.

Hering, D., Borja, A., Jones, J. I., Pont, D., Boets, P., Bouchez, A., ... & Leese, F. (2018).
Implementation options for DNA-based identification into ecological status assessment under the
European Water Framework Directive. Water Research, 138, 192-205.

Hijmans, R. J., Williams, E., Vennes, C., & Hijmans, M. R. J. (2017). Package ‘geosphere’.
Spherical trigonometry, 1, 7.

Hilty, J., & Merenlender, A. (2000). Faunal indicator taxa selection for monitoring ecosystem
health. Biological conservation, 92(2), 185-197.

Homer, C., Dewitz, J., Yang, L., Jin, S., Danielson, P., Xian, G., ... & Megown, K. (2015).
Completion of the 2011 National Land Cover Database for the conterminous United States–
representing a decade of land cover change information. Photogrammetric Engineering &
Remote Sensing, 81(5), 345-354.

Hosen, J. D., Febria, C. M., Crump, B. C., & Palmer, M. A. (2017). Watershed urbanization
linked to differences in stream bacterial community composition. Frontiers in microbiology, 8,
1452.

140
Hui, Cang, and Melodie A. McGeoch. "Zeta diversity as a concept and metric that unifies
incidence-based biodiversity patterns." The American Naturalist 184, no. 5 (2014): 684-694.

Hui, C., Vermeulen, W., & Durrheim, G. (2018). Quantifying multiple-site compositional
turnover in an Afrotemperate forest, using zeta diversity. Forest Ecosystems, 5(1), 1-9.

Huryn, A. D., Butz Huryn, V. M., Arbuckle, C. J., & Tsomides, L. (2002). Catchment land‐use,
macroinvertebrates and detritus processing in headwater streams: taxonomic richness versus
function. Freshwater Biology, 47(3), 401-415.

Ibekwe, A. M., Ma, J., & Murinda, S. E. (2016). Bacterial community composition and structure
in an Urban River impacted by different pollutant sources. Science of the Total Environment,
566, 1176-1185.

Jacobsen, D., Schultz, R., & Encalada, A. (1997). Structure and diversity of stream invertebrate
assemblages: the influence of temperature with altitude and latitude. Freshwater Biology, 38(2),
247-261.

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

141
Ji, Y., Ashton, L., Pedley, S. M., Edwards, D. P., Tang, Y., Nakamura, A., ... & Larsen, T. H.
(2013). Reliable, verifiable and efficient monitoring of biodiversity via metabarcoding. Ecology
letters, 16(10), 1245-1257.

Johnson, R. K., Furse, M. T., Hering, D., & Sandin, L. (2007). Ecological relationships between
stream communities and spatial scale: implications for designing catchment‐level monitoring
programmes. Freshwater Biology, 52(5), 939-958.

Johnson, R. K., & Hering, D. (2010). Spatial congruency of benthic diatom, invertebrate,
macrophyte, and fish assemblages in European streams. Ecological Applications, 20(4), 978-992.

Johnson, J. S., Spakowicz, D. J., Hong, B. Y., Petersen, L. M., Demkowicz, P., Chen, L., ... &
Sodergren, E. (2019). Evaluation of 16S rRNA gene sequencing for species and strain-level
microbiome analysis. Nature communications, 10(1), 1-11.

Jun, Y. C., Won, D. H., Lee, S. H., Kong, D. S., & Hwang, S. J. (2012). A multimetric benthic
macroinvertebrate index for the assessment of stream biotic integrity in Korea. International
journal of environmental research and public health, 9(10), 3599-3628.

Jyrkänkallio‐Mikkola, J., Heino, J., & Soininen, J. (2016). Beta diversity of stream diatoms at
two hierarchical spatial scales: implications for biomonitoring. Freshwater Biology, 61(2), 239-
250.

142
Keck, F., Bouchez, A., Franc, A., & Rimet, F. (2016). Linking phylogenetic similarity and
pollution sensitivity to develop ecological assessment methods: a test with river diatoms. Journal
of Applied Ecology, 53(3), 856-864.

Keck, F., Vasselon, V., Tapolczai, K., Rimet, F., & Bouchez, A. (2017). Freshwater
biomonitoring in the Information Age. Frontiers in Ecology and the Environment, 15(5), 266-
274.

Kermarrec, L., Franc, A., Rimet, F., Chaumeil, P., Humbert, J. F., & Bouchez, A. (2013). Next‐
generation sequencing to inventory taxonomic diversity in eukaryotic communities: a test for
freshwater diatoms. Molecular ecology resources, 13(4), 607-619.

Kermarrec, L., Franc, A., Rimet, F., Chaumeil, P., Frigerio, J. M., Humbert, J. F., & Bouchez, A.
(2014). A next-generation sequencing approach to river biomonitoring using benthic diatoms.
Freshwater Science, 33(1), 349-363.

Khalil, M. I., Gibson, D. J., Baer, S. G., & Willand, J. E. (2018). Functional diversity is more
sensitive to biotic filters than phylogenetic diversity during community assembly. Ecosphere,
9(3), e02164.

Kuhn, M. (2008). Building predictive models in R using the caret package. Journal of statistical
software, 28(5), 1-26.

143
Lan, Y., Wang, Q., Cole, J. R., & Rosen, G. L. (2012). Using the RDP classifier to predict
taxonomic novelty and reduce the search space for finding novel organisms. PLoS one, 7(3).

Lanzén, A., Lekang, K., Jonassen, I., Thompson, E. M., & Troedsson, C. (2016). High‐
throughput metabarcoding of eukaryotic diversity for environmental monitoring of offshore oil‐
drilling activities. Molecular ecology, 25(17), 4392-4406.

Laperriere, S. M., Hilderbrand, R. H., Keller, S. R., Trott, R., & Santoro, A. E. (2020).
Headwater stream microbial diversity and function across agricultural and urban land use
gradients. Applied and Environmental Microbiology.

Latombe, Guillaume, Cang Hui, and Melodie A. McGeoch. "Multi‐site generalised dissimilarity
modelling: using zeta diversity to differentiate drivers of turnover in rare and widespread
species." Methods in Ecology and Evolution 8, no. 4 (2017a): 431-442.

Latombe, G., McGeoch, M. A., Nipperess, D. A., & Hui, C. (2017b). Zetadiv: functions to
compute compositional turnover using zeta diversity. R package version 1.0.

Lazarina, Maria, Athanasios S. Kallimanis, Panayotis Dimopoulos, Maria Psaralexi, Danai-Eleni
Michailidou, and Stefanos P. Sgardelis. "Patterns and drivers of species richness and turnover of
neo-endemic and palaeo-endemic vascular plants in a Mediterranean hotspot: the case of Crete,
Greece." Journal of Biological Research-Thessaloniki 26, no. 1 (2019): 12.

144
Lear, G., Dopheide, A., Ancion, P., & Lewis, G. D. (2011). A comparison of bacterial, ciliate
and macroinvertebrate indicators of stream ecological health. Aquatic Ecology, 45(4), 517-527.

Leihy, Rachel I., Grant A. Duffy, and Steven L. Chown. "Species richness and turnover among
indigenous and introduced plants and insects of the Southern Ocean Islands." Ecosphere 9, no. 7
(2018): e02358.

Li, T., Huang, X., Jiang, X., & Wang, X. (2018). Assessment of ecosystem health of the Yellow
River with fish index of biotic integrity. Hydrobiologia, 814(1), 31-43.

Linard, B., Swenson, K., & Pardi, F. (2019). Rapid alignment-free phylogenetic identification of
metagenomic sequences. Bioinformatics, 35(18), 3303-3312.

Lu, H. P., Yeh, Y. C., Sastri, A. R., Shiah, F. K., Gong, G. C., & Hsieh, C. H. (2016). Evaluating
community–environment relationships along fine to broad taxonomic resolutions reveals
evolutionary forces underlying community assembly. The ISME journal, 10(12), 2867-2878.

Mande, S. S., Mohammed, M. H., & Ghosh, T. S. (2012). Classification of metagenomic
sequences: methods and challenges. Briefings in bioinformatics, 13(6), 669-681.

Marcy, Y., Ouverney, C., Bik, E. M., Lösekann, T., Ivanova, N., Martin, H. G., ... & Quake, S.
R. (2007). Dissecting biological “dark matter” with single-cell genetic analysis of rare and
145
uncultivated TM7 microbes from the human mouth. Proceedings of the National Academy of
Sciences, 104(29), 11889-11894.

Mazor, R. D., Rehn, A. C., Ode, P. R., Engeln, M., Schiff, K. C., Stein, E. D., ... & Hawkins, C.
P. (2016). Bioassessment in complex environments: designing an index for consistent meaning in
different settings. Freshwater Science, 35(1), 249-271.

McGeoch, Melodie A., Guillaume Latombe, Nigel R. Andrew, Shinichi Nakagawa, David A.
Nipperess, Mariona Roigé, Ezequiel M. Marzinelli et al. "Measuring continuous compositional
change using decline and decay in zeta diversity." Ecology 100, no. 11 (2019): e02832.

McMurdie, P. J., & Holmes, S. (2014). Waste not, want not: why rarefying microbiome data is
inadmissible. PLoS computational biology, 10(4).

Minerovic, A. D., Potapova, M. G., Sales, C. M., Price, J. R., & Enache, M. D. (2020). 18S-V9
DNA metabarcoding detects the effect of water-quality impairment on stream biofilm eukaryotic
assemblages. Ecological Indicators, 113, 106225.

Ode, P. R., Fetscher, A. E., & Busse, L. B. (2016). Standard Operating Procedures (SOP) for the
Collection of Field Data for Bioassessments of California Wadeable Streams: Benthic
Macroinvertebrates, Algae, and Physical Habitat. California State Water Resources Control
Board Surface Water Ambient Monitoring Program, Sacramento, Calif., USA.

146
Maloney, K. O., Munguia, P., & Mitchell, R. M. (2011). Anthropogenic disturbance and
landscape patterns affect diversity patterns of aquatic benthic macroinvertebrates. Journal of the
North American Benthological Society, 30(1), 284-295.

Mortágua, A., Vasselon, V., Oliveira, R., Elias, C., Chardon, C., Bouchez, A., ... & Almeida, S.
F. (2019). Applicability of DNA metabarcoding approach in the bioassessment of Portuguese
rivers using diatoms. Ecological Indicators, 106, 105470.

Ode, P. R. (2003). CAMLnet: list of California macroinvertebrate taxa and standard taxonomic
effort. Aquatic Bioassessment Laboratory, Rancho Cordova. Retrieved September, 10, 2014.

Ourso, R. T., & Frenzel, S. A. (2003). Identification of linear and threshold responses in streams
along a gradient of urbanization in Anchorage, Alaska. Hydrobiologia, 501(1-3), 117-131.

Padial, A. A., Ceschin, F., Declerck, S. A., De Meester, L., Bonecker, C. C., Lansac-Tôha, F. A.,
... & Bini, L. M. (2014). Dispersal ability determines the role of environmental, spatial and
temporal drivers of metacommunity structure. PloS one, 9(10).

Paller, M. H., Kosnicki, E., Prusha, B. A., Fletcher, D. E., Sefick, S. A., & Feminella, J. W.
(2017). Development of an index of biotic integrity for the Sand Hills ecoregion of the
southeastern United States. Transactions of the American Fisheries Society, 146(1), 112-127.

147
Pan, Y., Stevenson, R. J., Hill, B. H., Kaufmann, P. R., & Herlihy, A. T. (1999). Spatial patterns
and ecological determinants of benthic algal assemblages in Mid‐Atlantic streams, USA. Journal
of Phycology, 35(3), 460-468.

Parada, A. E., Needham, D. M., & Fuhrman, J. A. (2016). Every base matters: assessing small
subunit rRNA primers for marine microbiomes with mock communities, time series and global
field samples. Environmental microbiology, 18(5), 1403-1414.

Paruch, L., Paruch, A. M., Eiken, H. G., & Sørheim, R. (2019). Faecal pollution affects
abundance and diversity of aquatic microbial community in anthropo-zoogenically influenced
lotic ecosystems. Scientific Reports, 9(1), 1-13.

Passy, S. I., Bode, R. W., Carlson, D. M., & Novak, M. A. (2004). Comparative environmental
assessment in the studies of benthic diatom, macroinvertebrate, and fish communities.
International Review of Hydrobiology: A Journal Covering all Aspects of Limnology and
Marine Biology, 89(2), 121-138.

Passy, S. I., & Blanchet, F. G. (2007). Algal communities in human‐impacted stream ecosystems
suffer beta‐diversity decline. Diversity and Distributions, 13(6), 670-679.

Patrick, C. J., & Swan, C. M. (2011). Reconstructing the assembly of a stream-insect
metacommunity. Journal of the North American Benthological Society, 30(1), 259-272.

148
Paul, M. J., & Meyer, J. L. (2001). Streams in the urban landscape. Annual review of Ecology
and Systematics, 32(1), 333-365.

Pauls, S. U., Blahnik, R. J., Zhou, X., Wardwell, C. T., & Holzenthal, R. W. (2010). DNA
barcode data confirm new species and reveal cryptic diversity in Chilean Smicridea
(Smicridea)(Trichoptera: Hydropsychidae). Journal of the North American Benthological
Society, 29(3), 1058-1074.

Pawlowski, J., Kelly-Quinn, M., Altermatt, F., Apothéloz-Perret-Gentil, L., Beja, P., Boggero,
A., ... & Feio, M. J. (2018). The future of biotic indices in the ecogenomic era: Integrating (e)
DNA metabarcoding in biological assessment of aquatic ecosystems. Science of the Total
Environment, 637, 1295-1310.

Peck, D. V., Lazorchak, J. M., & Klemm, D. J. (Eds.). (2002). Environmental Monitoring and
Assessment Program--surface Waters: Western Pilot Study Field Operations Manual for
Wadeable Streams. National Health and Environmental Effects Research Laboratory [and]
National Exposure Research Laboratory, Office of Research and Development, US
Environmental Protection Agency.

Pond, Gregory J. "Biodiversity loss in Appalachian headwater streams (Kentucky, USA):
Plecoptera and Trichoptera communities." Hydrobiologia 679, no. 1 (2012): 97-117.

149
Ratnasingham, S., & Hebert, P. D. (2007). BOLD: The Barcode of Life Data System
(http://www. barcodinglife. org). Molecular ecology notes, 7(3), 355-364.

Rehn, A. C., Mazor, R. D., & Ode, P. R. (2015). The California stream condition index (CSCI): a
new statewide biological scoring tool for assessing the health of freshwater streams. Sacramento,
CA: Surface Water Ambient Monitoring Program (SWAMP).

Ren, S., Ahmed, N., Bertels, K., & Al-Ars, Z. (2018, October). An Efficient GPU-Based de
Bruijn Graph Construction Algorithm for Micro-Assembly. In 2018 IEEE 18th International
Conference on Bioinformatics and Bioengineering (BIBE) (pp. 67-72). IEEE.

Reynoldson, T. B., Norris, R. H., Resh, V. H., Day, K. E., & Rosenberg, D. M. (1997). The
reference condition: a comparison of multimetric and multivariate approaches to assess water-
quality impairment using benthic macroinvertebrates. Journal of the North American
Benthological Society, 16(4), 833-852.

Richards, A. B., & Rogers, D. C. (2006). List of the freshwater macroinvertebrate taxa from
California and adjacent states including standard taxonomic effort levels. Southwest Association
of Freshwater Invertebrate Taxonomists, Chico, California. Accessed online March, 1, 2008.

Rimet, F., Chaumeil, P., Keck, F., Kermarrec, L., Vasselon, V., Kahlert, M., ... & Bouchez, A.
(2016). R-Syst:: diatom: an open-access and curated barcode database for diatoms and
freshwater monitoring. Database, 2016.
150

Rinke, C., Schwientek, P., Sczyrba, A., Ivanova, N. N., Anderson, I. J., Cheng, J. F., ... &
Dodsworth, J. A. (2013). Insights into the phylogeny and coding potential of microbial dark
matter. Nature, 499(7459), 431-437.

Ristau, K., Steinfartz, S., & Traunspurger, W. (2013). First evidence of cryptic species diversity
and significant population structure in a widespread freshwater nematode morphospecies (T
obrilus gracilis). Molecular Ecology, 22(17), 4562-4575.

Ruhí, A., Datry, T., & Sabo, J. L. (2017). Interpreting beta‐diversity components over time to
conserve metacommunities in highly dynamic ecosystems. Conservation Biology, 31(6), 1459-
1468.

Schäfer, R. B., Bundschuh, M., Rouch, D. A., Szöcs, E., Peter, C., Pettigrove, V., ... & Kefford,
B. J. (2012). Effects of pesticide toxicity, salinity and other environmental variables on selected
ecosystem functions in streams and the relevance for ecosystem services. Science of the Total
Environment, 415, 69-78.

Scheiner, S. M., Chiarucci, A., Fox, G. A., Helmus, M. R., McGlinn, D. J., & Willig, M. R.
(2011). The underpinnings of the relationship of species richness with space and time. Ecological
Monographs, 81(2), 195-213.

151
Siefert, A., Ravenscroft, C., Weiser, M. D., & Swenson, N. G. (2013). Functional beta‐diversity
patterns reveal deterministic community assembly processes in eastern North American trees.
Global Ecology and Biogeography, 22(6), 682-691.

Simonin, M., Voss, K. A., Hassett, B. A., Rocca, J. D., Wang, S. Y., Bier, R. L., ... & Bernhardt,
E. S. (2019). In search of microbial indicator taxa: shifts in stream bacterial communities along
an urbanization gradient. Environmental microbiology, 21(10), 3653-3668.

Simons, Ariel Levi, Raphael Mazor, Eric D. Stein, and Sergey Nuzhdin. "Using alpha, beta, and
zeta diversity in describing the health of stream‐based benthic macroinvertebrate communities."
Ecological Applications 29, no. 4 (2019): e01896.

Socolar, Jacob B., James J. Gilroy, William E. Kunin, and David P. Edwards. "How should beta-
diversity inform biodiversity conservation?." Trends in ecology & evolution 31, no. 1 (2016):
67-80.

Soininen, J. (2007). Environmental and spatial control of freshwater diatoms—a review. Diatom
research, 22(2), 473-490.

Soininen, J., Paavola, R., & Muotka, T. (2004). Benthic diatom communities in boreal streams:
community structure in relation to environmental and spatial gradients. Ecography, 27(3), 330-
342.

152
Somervuo, P., Yu, D. W., Xu, C. C., Ji, Y., Hultman, J., Wirta, H., & Ovaskainen, O. (2017).
Quantifying uncertainty of taxonomic placement in DNA barcoding and metabarcoding.
Methods in Ecology and Evolution, 8(4), 398-407.

Sponseller, R. A., Benfield, E. F., & Valett, H. M. (2001). Relationships between land use,
spatial scale and stream macroinvertebrate communities. Freshwater biology, 46(10), 1409-1424.

Stein, E. D., Martinez, M. C., Stiles, S., Miller, P. E., & Zakharov, E. V. (2014a). Is DNA
barcoding actually cheaper and faster than traditional morphological methods: results from a
survey of freshwater bioassessment efforts in the United States?. PloS one, 9(4).

Stein, E. D., White, B. P., Mazor, R. D., Jackson, J. K., Battle, J. M., Miller, P. E., ... &
Sweeney, B. W. (2014b). Does DNA barcoding improve performance of traditional stream
bioassessment metrics?. Freshwater Science, 33(1), 302-311.

Stepenuck, K. F., Crunkilton, R. L., & Wang, L. (2002). Impacts of urban landuse on
macroinvertebrate communities in southeastern Wisconsin streams 1. JAWRA Journal of the
American Water Resources Association, 38(4), 1041-1051.

Sweeney, B. W., Battle, J. M., Jackson, J. K., & Dapkey, T. (2011). Can DNA barcodes of
stream macroinvertebrates improve descriptions of community structure and water quality?.
Journal of the North American Benthological Society, 30(1), 195-216.

153
Taberlet, P., Coissac, E., Pompanon, F., Brochmann, C., & Willerslev, E. (2012). Towards next‐
generation biodiversity assessment using DNA metabarcoding. Molecular ecology, 21(8), 2045-
2050.

Tang, K., Ren, J., Cronn, R., Erickson, D. L., Milligan, B. G., Parker-Forney, M., ... & Sun, F.
(2018). Alignment-free genome comparison enables accurate geographic sourcing of white oak
DNA. BMC genomics, 19(1), 896.

Team, R. C. (2017). R: A Language and environment for statistical computing, v. 3.5. 1 (The R
Foundation for Statistical Computing, Vienna). Available at www. r-project. org. Accessed
January, 3, 2018.

Valentin, V., Frédéric, R., Isabelle, D., Olivier, M., Yorick, R., & Agnès, B. (2019). Assessing
pollution of aquatic environments with diatoms’ DNA metabarcoding: experience and
developments from France Water Framework Directive networks. Metabarcoding and
Metagenomics, 3, e39646.

Van der Gucht, K., Cottenie, K., Muylaert, K., Vloemans, N., Cousin, S., Declerck, S., ... &
Degans, H. (2007). The power of species sorting: local factors drive bacterial community
composition over a wide range of spatial scales. Proceedings of the National Academy of
Sciences, 104(51), 20404-20409.

154
Vanschoenwinkel, Bram, Chris De Vries, Maitland Seaman, and Luc Brendonck. "The role of
metacommunity processes in shaping invertebrate rock pool communities along a dispersal
gradient." Oikos 116, no. 8 (2007): 1255-1266.

Vasselon, V., Domaizon, I., Rimet, F., Kahlert, M., & Bouchez, A. (2017). Application of high-
throughput sequencing (HTS) metabarcoding to diatom biomonitoring: Do DNA extraction
methods matter?. Freshwater Science, 36(1), 162-177.

Vázquez, G., Aké-Castillo, J. A., & Favila, M. E. (2011). Algal assemblages and their
relationship with water quality in tropical Mexican streams with different land uses.
Hydrobiologia, 667(1), 173-189.

Vile, J. S., & Henning, B. F. (2018). Development of indices of biotic integrity for high-gradient
wadeable rivers and headwater streams in New Jersey. Ecological Indicators, 90, 469-484.

von Ammon, U., Wood, S. A., Laroche, O., Zaiko, A., Tait, L., Lavery, S., ... & Pochon, X.
(2018). Combining morpho-taxonomy and metabarcoding enhances the detection of non-
indigenous marine pests in biofouling communities. Scientific reports, 8(1), 1-11.

Wakelin, S. A., Colloff, M. J., & Kookana, R. S. (2008). Effect of wastewater treatment plant
effluent on microbial function and community structure in the sediment of a freshwater stream
with variable seasonal flow. Appl. Environ. Microbiol., 74(9), 2659-2668.

155
Wallace, A. M., & Biastoch, R. G. (2016). Detecting changes in the benthic invertebrate
community in response to increasing chloride in streams in Toronto, Canada. Freshwater
Science, 35(1), 353-363.

Ward, J. V. (1986). Altitudinal zonation in a Rocky Moutain stream. Archiv für Hydrobiologie.
Supplementband. Monographische Beiträge, 74(2), 133-199.

Wolf, D. I., & Vis, M. L. (2020). Stream Algal Biofilm Community Diversity Along An Acid
Mine Drainage Recovery Gradient Using Multimarker Metabarcoding. Journal of Phycology,
56(1), 11-22.

Wu, N., Qu, Y., Guse, B., Makarevičiūtė, K., To, S., Riis, T., & Fohrer, N. (2018). Hydrological
and environmental variables outperform spatial factors in structuring species, trait composition,
and beta diversity of pelagic algae. Ecology and evolution, 8(5), 2947-2961.

Xie, Y., Hong, S., Kim, S., Zhang, X., Yang, J., Giesy, J. P., ... & Khim, J. S. (2017).
Ecogenomic responses of benthic communities under multiple stressors along the marine and
adjacent riverine areas of northern Bohai Sea, China. Chemosphere, 172, 166-174.

Yang, B., Wang, Y., & Qian, P. Y. (2016). Sensitivity and correlation of hypervariable regions in
16S rRNA genes in phylogenetic analysis. BMC bioinformatics, 17(1), 135.

156
Zhan, A., Hulák, M., Sylvester, F., Huang, X., Adebayo, A. A., Abbott, C. L., ... & MacIsaac, H.
J. (2013). High sensitivity of 454 pyrosequencing for detection of rare species in aquatic
communities. Methods in Ecology and Evolution, 4(6), 558-565.

Zielezinski, A., Girgis, H. Z., Bernard, G., Leimeister, C. A., Tang, K., Dencker, T., ... & Comin,
M. (2019). Benchmarking of alignment-free sequence comparison methods. Genome biology,
20(1), 144.

Conclusion
This dissertation covers a series of investigations on relationships between novel
measures of landscape diversity for lotic communities across California, environmental
gradients, and assessments of their biotic integrity. Chapter 1 focused on patterns of zeta
diversity for communities of morphologically sorted benthic macroinvertebrates (BMIs), and
found that they could explain most of the observed regional variation in their biotic integrity.
Chapter 2 approaches the same data as chapter 1, but with a focus on how changes in the
topology of spatial co-occurrence networks for lotic BMIs can be used as a basis for modeling
their regional biotic integrity. Chapter 3 revists the use of zeta diversity in modeling regional
biotic integrity, but with a focus on a diverse array of biological communities as classified with
the use of environmental DNA.
For all three chapters we used a variety of data sets, collected by the Southern California
Coastal Water Research Project (SCCWRP), describing both the compositions of lotic
communities as well as measures of their physical environment. The source data for the analyses
157
within the first two chapters were composed of morphologically sorted BMIs. In the last chapter
we used a number of metabarcoding primers to describe the composition of communities of
BMIs as well as diatoms, soft algae (non-diatomaceous algae), and prokaryotes.
With metabarcoding a much more diverse set of biological communities can be classified
than with established taxonomic methods. The potential breadth of ecological data which
metabarcoding can capture, along with its relative low cost, underlies our motivation in
developing methods of bioassessment which build upon existing methods of bioassessment to
describe the ecological health of lotic communities ranging from prokaryotes to BMIs. We
investigated both zeta diversity and topological measures of co-occurrence networks as a means
of developing tools which could investigate if predictable relationships exist between ecological
stress and patterns of landscape diversity. We found BMIs, sorted morphologically or classified
via metabarcoding, tended follow a pattern whereby increased ecological stress, due to upstream
land use, was associated with a decline both in local lotic community diversity as well as the
average amount of taxonomic overlap between impacted communities. For our prokaryotic
communities, classified using two different 16S primers, we tended to see their compositions
diverge under ecological stress, although the signal was weaker and less consistent across
taxonomic levels than with BMIs. With soft algae we tended to find community compositions
converging under ecological stress. We found our communities of diatoms behaved in a similar
fashion as BMIs, and that changes in their patterns of landscape diversity could even account for
more observed variation in regional biotic integrity, although only when using diatoms classified
with a diatom specific rbcL primer.
Our results in chapter 3 point to the potential for extending our methods in chapters 1 and
2 to assess the regional biotic integrity of lotic communities described via metabarcoding.
158
However, our findings in chapter 3 also point to limitations with the use of communities
described via metabarcoding. While relationships between our measures of landscape diversity
and ecological stress appear similar between communities of prokaryotes classified with
different primers, or between communities of BMIs classified via metabarcoding or
morphological sorting, for diatoms we found significant differences depending on the primers
used in metabarcoding based classifications. These differences may be due to the ability of
different primers to fully capture the diversity of different biological communities, such as BMIs
versus diatoms, given different lengths for metabarcoded regions. We also found large
differences in the completeness of the various databases we used to assign taxonomic identities
to our metabarcoded sequences. Limitations such as these may underlie the differences we found
in patterns of landscape diversity of diatoms classified via an 18S V9 versus a rbcL primer, given
the much greater taxonomic resolution of the later.
Although commonly used methods of classifying biological communities via
metabarcoding have a number of limitations, future work may be able to overcome them and
more fully capture measures of diversity through various techniques such as shotgun
metagenomics and reference-free approaches to classifying sequences. Ultimately, we find that
with continuing advances in sequencing and computational technology there is a great potential
in developing new techniques in bioassessment which synthesize measures of landscape ecology
with environmental DNA.
159

Appendix A: supplementary material

160
Chapter 1: supplementary material

Table S1: Linear models of taxonomic diversity measures with land use, altitude, and distance.
Diversity measure α taxa β taxa ζ 10,taxa b taxa
Intercept 45.5 7.3 x 10
-1
1.9 -1.6
Coefficient (altitude) -9.3 x 10
-4
-2.4 x 10
-5
1.5 x 10
-4
1.5 x 10
-4
Coefficient (land use) -3.6 x 10
-1
7.1 x 10
-4
-1.8 x 10
-2
-4.6 x 10
-3
Coefficient (distance) 4.5 x 10
-5
1.5 x 10
-6
-2.3 x 10
-5
-1.4 x 10
-5
Proportion of variation due
to altitude
F(1,5596)
1909

1166 308 348
P-value, variation due to
altitude
< 10
-4
< 10
-4
< 10
-4
< 10
-4

Proportion of variation due
to land use
F(1,5596)
6003 103 180 11.7
P-value, variation due to
land use
< 10
-4
< 10
-4
< 10
-4
< 10
-3
Proportion of variation due
to distance
F(1,5596)
29.5 578 207 387
P-value, variation due to
distance
< 10
-4
< 10
-4
< 10
-4
< 10
-4

161
Relative importance of
altitude (%)
7.1 10.5 2.7 3.8
Relative importance of land
use (%)
47.3 7.7 6.1 2.5
Relative importance of
distance
(%)
4.3 6.6 2.2 5.5
Proportion of variance
explained by model (%)
58.7 24.8 11.0 11.8

162
Table S2: Linear models of functional diversity measures with land use, altitude, and distance.
Diversity measure α functional β functional ζ 10,functional b functional
Intercept 6.6 2.0 x 10
-1
4.3 -1.9 x 10
-1
Coefficient (altitude) -3.1 x 10
-6
-6.0 x 10
-6
1.1 x 10
-4
1.1 x 10
-5
Coefficient (land use) -2.1 x 10
-2
1.6 x 10
-3
-2.7 x 10
-2
-2.3 x 10
-3

Coefficient (distance) 2.9 x 10
-6
4.7 x 10
-7
-2.8 x 10
-6
-6.0 x 10
-7
Proportion of variation due to
altitude
F(1,5596)
2252 774 1068 454
P-value, variation due to
altitude
< 10
-4
< 10
-4
< 10
-4
< 10
-4

Proportion of variation due to
land use
F(1,5596)
5497 1113 1610 599
P-value, variation due to land
use
< 10
-4
< 10
-4
< 10
-4
< 10
-4

Proportion of variation due to
distance
F(1,5596)
31.0 35.6 5.9 15.0
P-value, variation due to
distance
< 10
-4
< 10
-4
< 5 x 10
-2
< 10
-3

Relative importance of
altitude (%)
8.3 5.2 6.6 3.5
Relative importance of land
use (%)
45.8 19.8 24.9 12.2
163
Relative importance of
distance
(%)
4.0 0.6 1.0 0.4
Proportion of variance
explained by model (%)
58.2 25.6 32.4 16.0

164
Table S3: Linear models of taxonomic β and ζ 10 diversity measures with land use and altitude.
Diversity measure β taxa β functional ζ 10,taxa ζ 10,functional
Intercept 7.3 x 10
-1
2.1 x 10
-1
1.7 3.9
Coefficient (altitude) 1.0 x 10
-5
-7.1 x 10
-5
3.4 x 10
-4
2.4 x 10
-3

Coefficient (standard
deviation of land use)
2.4 x 10
-3
2.8 x 10
-3
-4.5 x 10
-2
-4.0 x 10
-2

Proportion of variation due to
altitude
F(1,2797)
24.4 4.4 2.0 76.0
P-value, variation due to
altitude
< 10
-4
< 5 x 10
-2
n.s. < 10
-4

Proportion of variation due to
the standard deviation of land
use
F(1,2797)
1175 646.9 707.4 647.3
P-value, variation due to the
standard deviation of land use
< 10
-4
< 10
-4
< 10
-4
< 10
-4

Relative importance of
altitude (%)
3.2 x 10
-1
4.1 x 10
-1
6.2 x 10
-2
2.9
Relative importance of the
standard deviation of land use
(%)
29.7 18.5 20.2 17.6
Proportion of variance
explained by model (%)
30.0 18.9 20.2 20.6

165

166
Table S4: Linear models of taxonomic β and ζ 10 diversity measures with land use and altitude.
Diversity measure β taxa β functional ζ 10,taxa ζ 10,functional
Intercept 7.6 x 10
-1
1.8 x 10
-1
1.5 4.6
Coefficient (altitude) -2.7 x 10
-5
-2.2 x 10
-6
2.4 x 10
-4
1.0 x 10
-4

Coefficient (standard
deviation of land use)
2.6 x 10
-3
3.5 x 10
-3
-4.3 x 10
-2
-5.0 x 10
-2

Proportion of variation due to
altitude
F(1,2797)
1143 522.5 202.5 478.6

P-value, variation due to
altitude
< 10
-4
< 10
-4
< 10
-4
< 10
-4

Proportion of variation due to
the standard deviation of land
use
F(1,2797)
984.5 1441 282.7 1059
P-value, variation due to the
standard deviation of land use
< 10
-4
< 10
-4
< 10
-4
< 10
-4

Relative importance of
altitude (%)
13.5 5.5 3.3 5.6
Relative importance of the
standard deviation of land use
(%)
29.7 35.7 11.5 29.9
Proportion of variance
explained by model (%)
43.2 41.2 14.8 35.5

167

168

Figure S1: Trends in the mean values per HUC 8 watershed of α functional, β functional, and ζ 10,functional versus changes in
altitude, land use, and distance. Note that a rise in β or a decline in ζ 10 diversity indicates a greater level of
dissimilarity between samples within a group.

Abstract (if available)

Abstract Freshwater streams are of great ecological importance, helping to regulate a variety of biogeochemical cycles. While the ecosystem services provided by streams help to support a number of human activities, they are frequently degraded by a number of anthropogenic stressors. Efforts to monitor the ecological health of stream ecosystems have frequently focused on measurements based on the local diversity of benthic macroinvertebrates (BMIs) sorted by morphology. However, individual BMI communities are effectively linked across a landscape through processes such as dispersal, and as such there is a need to measure the health of multiple communities on a regional basis. Here we investigate the potential role of three techniques in development of such regional assessments: zeta diversity, topological measures of co-occurrence networks, and environmental DNA (eDNA). Through our investigations we were able to demonstrate the following key findings: (1) patterns of landscape diversity, as described by either zeta diversity or topological measures of co-occurrence networks, can account for a significant proportion in the regional health of watersheds as described by the California Streams Condition Index (CSCI), (2) eDNA is a viable alternative to classifying and organizing biological communities found in freshwater streams. Our results indicate the potential for the development of regional scale assessments of the health of freshwater streams which can be derived novel patterns of landscape diversity in biological communities classified via eDNA.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Marine protistan diversity, spatiotemporal dynamics, and physiological responses to environmental cues

PDF

Spatial and temporal dynamics of marine microbial communities and their diazotrophs in the Southern California Bight

PDF

Annual pattern and response of the bacterial and microbial eukaryotic communities in an aquatic ecosystem restructured by disturbance

PDF

Future impacts of warming and other global change variables on phytoplankton communities of coastal Antarctica and California

Asset Metadata

Creator Simons, Ariel Levi (author)

Core Title The development of novel measures of landscape diversity in assessing the biotic integrity of lotic communities

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Biology (Marine Biology and Biological Oceanography)

Publication Date 09/21/2020

Defense Date 08/12/2020

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag benthic macroinvertebrates,biotic integrity,co-occurrence networks,environmental DNA,landscape diversity,landscape ecology,OAI-PMH Harvest,stream ecosystems,zeta diversity

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Nuzhdin, Sergey (committee chair), Fuhrman, Jed (committee member), Kenkel, Carly (committee member), Leach, William (committee member)

Creator Email alsimons@usc.edu,levisimons@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-373243

Unique identifier UC11666202

Identifier etd-SimonsArie-8998.pdf (filename),usctheses-c89-373243 (legacy record id)

Legacy Identifier etd-SimonsArie-8998.pdf

Dmrecord 373243

Document Type Dissertation

Rights Simons, Ariel Levi

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA