Dimensions of biodiversity in Chesapeake Bay demersal fishes: patterns and drivers through space and time

Biodiversity has typically been described in terms of species richness and composition, but theory and growing empirical evidence indicate that the diversity of functional traits, the breadth of evolutionary relationships, and the equitability with which individuals or biomass are distributed among species better characterize patterns and processes within ecosystems. Yet, the advantages of including such data come at the expense of measuring traits, sequencing genes, and counting or weighing individuals, and it remains unclear whether this greater resolution yields substantial benefits in describing diversity. We summarized a decade of high-resolution trawl data from a bimonthly trawl survey to investigate spatial and seasonal patterns of demersal fish diversity in the Chesapeake Bay, USA, with the goal of identifying areas and times of mismatch between different dimensions of diversity, and their response to environmental forcing. We found moderate to strong positive relationships among all metrics of diversity, and that functional and phylogenetic differences were well-reflected in an index derived from taxonomic (Linnaean) hierarchy. Compared with species richness and species diversity, functional, phylogenetic, and taxonomic indices peaked later in the year, which was a consequence of the distribution of biomass among functionally and evolutionarily divergent species. Generalized additive models revealed that spatial, temporal, and environmental variables explained roughly similar proportions of deviance across all aspects of diversity, suggesting that these three factors do not differentially affect the functional and phylogenetic aspects of community structure. We conclude that an index of diversity derived from taxonomic hierarchy served well as a practical surrogate for functional and phylogenetic diversity of the demersal fish community in this system. We also emphasize the importance of evenness in understanding diversity patterns, especially since most ecological communities in nature are dominated by one or few species.


INTRODUCTION
Biodiversity has long been described in terms of species richness and composition, but it is increasingly recognized that species' roles in ecosystems are intrinsically linked to their functional and phylogenetic differences. Unlike traditional measures of diversity, functional and phylogenetic approaches recognize that species vary by degrees, i.e., that some species are more functionally or phylogenetically similar than others (Webb et al. 2002, McGill et al. 2006. Consequently, these approaches allow for a finer understanding of ecological patterns, the mechanisms driving observed patterns and maintenance of diversity in nature (Webb et al. 2002, Petchey andGaston 2006), and the role of diversity in ecosystem processes (Hooper et al. 2005, Srivastava et al. 2012). Yet, these advantages come with a tradeoff in effort, as calculations of functional and phylogenetic diversity require considerable additional data beyond simple species counts. While functional and phylogenetic methods are being increasingly applied across systems and taxa (Cadotte et al. 2011), it remains unclear when and where investment in functional and phylogenetic information yields the greatest benefits when characterizing biodiversity.
Functional diversity (hereafter FD) quantifies the community-wide variation in organismal functional traits. For our study, we define a functional trait as any character that influences how an organism interacts with its environment and/or other organisms. Our definition is closest to that proposed by Violle et al. (2007) in that we consider functional traits to be morphological, physiological, and life historical indicators of fitness at the individual level. However, we also include traits related to performance that contribute to ecosystem processes as well, such as body mass, and thus more closely resemble effect traits as defined by Lavorel and Garnier (2002). Incorporating information on multiple kinds of traits (functional, effect, etc.) can provide deeper insight into differences among organisms at the community level (e.g., Petchey et al. 2007, Devictor et al. 2010. A trait-based approach to diversity is attractive because functional traits can, in principle, be directly linked to ecosystem processes (Díaz and Cabido 2001). Thus, FD may provide stronger predictions about those processes than richness alone (e.g., Cadotte et al. 2009). Calculation of FD requires identification of traits that are relevant to the pattern or process under investigation, although there is often a lack of consensus as to what constitute relevant traits for many taxa (Petchey and Gaston 2006). Since trait measurements can be a costly and/or time-intensive process, most traits used in calculations of FD to date are relatively easy to obtain. However, they may not provide the most direct link to fitness or performance (so-called soft traits, sensu Hodgson et al. 1999).
Phylogenetic diversity (hereafter PD) quantifies the amount of unshared evolutionary history within an assemblage, and is used with the assumption that underlying evolutionary variation among organisms is an indicator of functional dissimilarity (Webb et al. 2002). As such, PD is often invoked as a holistic measure of diversity, especially if relevant functional traits are not known or measured (Srivastava et al. 2012). For this reason, PD has been shown to be a better predictor of ecosystem functioning than FD in some cases (Cadotte et al. 2009, but see Flynn et al. 2011). This approach assumes that functional differences are evolutionarily conserved, but a comprehensive review of 26 different datasets showed that only 60% of 103 ecological traits exhibited significant phylogenetic signal (Freckleton et al. 2002). These synthetic results, coupled with those from more recent individual studies (e.g., Flynn et al. 2011), suggest that the utility of PD as a surrogate for FD may vary depending on the organisms, their traits, and the scale of the investigation (Srivastava et al. 2012).
A larger roadblock to the widespread use of phylogenetic approaches is that robust molecular phylogenies are not always available, and DNA sequencing can be cost prohibitive or difficult. Several studies have suggested substituting existing Linnaean taxonomy as a proxy for phylogeny (Clarke and Warwick 1998, Shimatani 2001, Crozier et al. 2005, Ricotta 2005, Maherali and Klironomos 2007, Ricotta et al. 2012), or at least using taxonomy to inform poorly resolved phylogenies (Davies et al. 2012). With taxonomy, researchers can easily draw on centuries worth of rigorous data. However, taxonomic distances may not wholly reflect evolutionary relationships suggested by molecular phylogeny, especially when species have originated recently and/or radiated rapidly (Crozier et al. 2005). Moreover, discrete taxonomic levels (e.g., genus, family, etc.) limit the resolution of evolutionary differences compared to those derived from continuous DNA sequence variation and thus may not provide much insight when species within an assemblage are closely related. In the case of microbes, evolutionary relationships may not be easily defined without the use of molecular tools (Wayne et al. 1987). There have been few direct tests of whether measures of diversity derived from a taxonomic hierarchy-which we refer to as taxonomic diversity (hereafter, TD)-provide similar inferences to PD (see Ricotta et al. 2012).
In addition to functional and phylogenetic differences among species, the relative abundance or biomass of species, defined as evenness, also influences ecosystem processes and may have strong implications for the interpretation of FD and PD (reviewed in Hillebrand et al. 2008). Evidence suggests that communities with a more even distribution of individuals or biomass among species alters biodiversity patterns (Stuart-Smith et al. 2013) and enhances ecosystem processes Potvin 2000, Kirwan et al. 2007), although there have been comparatively few examples of how evenness modifies the inferences made from FD or PD in particular. In terrestrial grasslands, indices of FD and PD that did not incorporate relative abundance were better predictors of ecosystem processes than those that did (Mouillot et al. 2011b, Cadotte et al. 2012. Similarly, the inclusion or exclusion of relative abundance information did not alter areas of spatial mismatch of FD and PD in a bird assemblage (Devictor et al. 2010). However, turnover in dominant species altered the functional structure of estuarine fish communities after disturbance, despite no change in species richness (Villéger et al. 2010). Additionally, incorporating information on both functional traits and relative abundances shifted global diversity hotspots of reef fish diversity towards temperate regions (Stuart-Smith et al. 2013).
Thus, it appears that in at least some circumstances, evenness may play an important role in altering the relationships between FD, and potentially PD, and more traditional indices of diversity.
Here, we investigated how different dimensions of diversity explain spatial and seasonal patterns in the demersal fish community in the Chesapeake Bay, the largest estuary in the USA. Chesapeake Bay is characterized by a strong spatial gradient in salinity driven by tidal exchange with the Atlantic Ocean at the mouth and riverine input, and experiences some of the most extreme intra-annual fluctuations in temperature in the world (Murdy and Musick 2013). The Bay provides valuable economic, ecological, and recreational resources, and supports numerous fisheries (Murdy and Musick 2013). Approximately 350 fish species have been reported in the Chesapeake Bay, although many of these species are vagrants from freshwater or marine habitats, and are not found regularly (Jung andHoude 2003, Murdy andMusick 2013). Some species are only present in the Bay for a portion of the year to feed and/or spawn, spending the rest of their time offshore (Murdy and Musick 2013). We focus here on 50 species of bottom-dwelling fishes that represent the average annual demersal fish community in the Chesapeake Bay.
A recent study characterized spatial and temporal patterns in species richness and Simpson diversity in Chesapeake Bay demersal fishes using high-resolution data from 10 years of fisheries-independent bottom trawl data (Buchheister et al. 2013). They found that species richness and composition were primarily structured along strong environmental gradients, principally salinity, temperature, and dissolved oxygen, with highest richness at the northern and southern ends of the Bay's mainstem channel (Buchheister et al. 2013). Given that an increasing number of studies have found strong disagreement in spatial patterns between species richness, functional, and/or phylogenetic diversity (Petchey et al. 2007, Devictor et al. 2010, Mouillot et al. 2011a, Safi et al. 2011, Stuart-Smith et al. 2013, we supplemented the trawl survey data presented in Buchheister et al. (2013) with newly collected functional trait and genetic information.
Our specific objectives were to compare and contrast patterns in species richness, evenness, Gini-Simpson, functional, phylogenetic, and taxonomic diversity through space and time, and to identify the unique contributions of space, time, and environmental variables to diversity using v www.esajournals.org generalized additive models and deviance partitioning. We expected to find areas of spatial and temporal mismatch among the different measures of diversity based on previous investigations of estuarine (Villéger et al. 2010(Villéger et al. , 2012 and marine fishes (Mouillot et al. 2011a), as well as other ecological communities (Petchey et al. 2007, Devictor et al. 2010, Safi et al. 2011, Stuart-Smith et al. 2013. Based on well-delineated evolutionary relationships among species within this system, and evidence from other systems (Ricotta et al. 2012), we expected PD and TD to be suitable proxies for one another. Given the strong dominance of a small number of species in this system (Jung andHoude 2003, Buchheister et al. 2013), we predicted evenness would play a large role in determining diversity patterns. Finally, we expected to find strong spatial, temporal, and environmental signal in patterns of diversity, as has been shown previously with richness in this system (Buchheister et al. 2013) and estuarine systems in general (Odum 1988).

ChesMMAP survey
Data for this study were obtained from the Chesapeake Bay Multispecies Monitoring and Assessment Program (ChesMMAP; Latour et al. 2003, Bonzek et al. 2010. This program conducts a bottom trawl survey that operates bimonthly from March to November, and monitors the relative abundance, biomass, age-and sizestructure, and trophic interactions of demersal fishes within the mainstem waters of the Chesapeake Bay. Here, we present analyses of data from the years 2002-2011. The ChesMMAP survey employs a stratified random design, with strata determined by water depth and latitude, the latter dividing the Bay roughly by salinity regime (ranging from euryhaline at the mouth to oligohaline at the head). The survey samples approximately 80 randomly selected stations per bimonthly cruise, with the number of stations proportional to the area of each stratum. At each station, a twenty-minute tow is performed with an otter trawl (7.6-cm mesh in the cod end), which targets late-juvenile and adult fishes. Vessel GPS and trawl mensuration gear are used to quantify the area swept for each tow. Finally, measurements of temperature, salinity, dissolved oxygen, and depth are also collected at each station.
Following Buchheister et al. (2013), we omitted any species that were not adequately sampled by the survey gear (i.e., excluded all pelagic species and any demersal species with , 10 individuals or , 1 kg captured). This exercise left a total of 50 species for this analysis that together accounted for . 99% of the total biomass sampled by the survey (see Table A1 for the full list of species used in the analyses). We further omitted any tows that contained 1 species (19% of all tows), as these samples did not inform us about diversity and yielded values of evenness that were not interpretable. Ultimately, we were left with a total of 2,530 independent tows across all months and years.

Calculating diversity
To calculate diversity, we used Rao's quadratic entropy (hereafter Q, Rao 1982). Rao's Q is calculated as: where S is the number of species (richness), p i is the relative abundance of species i, p j is the relative abundance of species j, and d ij is the distance (functional, phylogenetic, etc.) between species i and j. Q is therefore the average difference between individuals in each tow. Most relevantly, the distances d ij can be derived from any measure of dissimilarity, including functional and phylogenetic (Shimatani 2001, Botta-Dukát 2005. Thus, Rao's Q provides a unifying calculation with which to investigate and compare across the different dimensions of diversity. For the analysis presented in the main text, we chose to weight by relative biomass instead of relative abundance. Biomass is a better proxy for ecological processes such as production and trophic transfer (Brett and Groves 1979), and has recently been shown to better predict biodiversity patterns in marine demersal fauna than abundance (Certain et al. 2013). However, we also calculated Q using relative abundance and presence-absence data to facilitate comparisons among the different weighting schemes, and to relate back to previous studies.
For purposes of interpretability and comparav www.esajournals.org bility, we converted values of Q into 'effective numbers' (Jost 2006) using the following equation: This transformation yields the 'effective number' of species in the sample if all species were equally abundant and maximally distinct. Effective numbers have several ideal attributes, including the 'doubling property,' which states that if two communities have the same diversity X but no shared species, then their combined diversity is 2X (Jost 2006). Additionally, as long as distances are ultrametric (i.e., between 0 and 1), then effective numbers can be interpreted identically and directly compared across all dimensions of diversity, including richness (Pavoine et al. 2005, Jost 2006).
We calculated richness as the number of species captured at a given station. For our index of species diversity, we standardized the distances d ij to either 0 (individuals belong to the same species) or 1 (different species). In this case, Q collapses to the well-known Gini-Simpson index (Ricotta and Szeidl 2006). In the event that species are all present in equal proportions, then the Gini-Simpson index will be equal to species richness. For evenness, we calculated the relative evenness index from Jost (2010): The relative evenness index is analogous to Pielou's evenness and is bound between 0 (maximally uneven) and 1 (maximally even).
To derive functional distances, we first measured and scored 25 functional traits for all 50 species. Functional traits included continuous measures of morphology (e.g., Villéger et al. 2010), estimates of diet composition (e.g., Petchey et al. 2007, Devictor et al. 2010, and life history characteristics. As the 25 traits contained both continuous and categorical values, we calculated Gower distances (Gower 1971) using the correction from Podani (1999) to account for ordered traits. We then employed hierarchical agglomerative clustering to produce an ultrametric functional dendrogram (Petchey and Gaston 2002). To account for the sensitivity of the dendrogram to the clustering algorithm used, we employed the approach suggested by Mouchet et al. (2008) of using multiple algorithms and choosing the one that best preserves the original, non-ultrametric distances ). Finally, we extracted the ultrametric distance matrix from the dendrogram and scaled by the maximum so that all values were between 0-1 before calculating Q. (See Appendix A for a more detailed discussion of trait selection, measurement, and dendrogram reconstruction).
To derive phylogenetic distances, we constructed a gene tree as a proxy for a molecular phylogeny, using cytochrome oxidase I (COI ) sequences. Sequence data were available on GenBank for 34 species in the dataset. We extracted, amplified, and sequenced COI for a further 3 species, and substituted congeneric relatives for the remaining 13 species. We aligned all sequences, including five outgroups from GenBank, using the online alignment tool MUS-CLE (Edgar 2004). The aligned sequences were used to construct a maximum-likelihood gene tree using the program RAxML (v.2.2.3, Stamatakis 2006). To convert the tree to ultrametric, we estimated absolute evolutionary rates using penalized-likelihood and verified the rates using a data-driven cross-validation criterion (Sanderson 2002). Scaling by the absolute rate fixed the age of the root at 1, making the tree ultrametric. The pairwise distance matrix was extracted from the ultrametric tree and divided by the maximum so all values were between 0 and 1 before calculating Q. (See Appendix A for a detailed discussion of sequencing protocol, congener selection, and phylogeny construction.) To test whether molecular phylogenetic data provided more information than Linnaean taxonomy, we obtained classification information for all 50 species using the latest accepted taxonomy from the Integrated Taxonomic Information System (http://www.itis.gov). Distances were computed as the number of nodes separating each pair of species on the hierarchy (Clarke and Warwick 1998), weighted by the proportional decrease in taxon richness at each successive level on the hierarchy (Clarke and Warwick 1999). Since these distances were already ultrametric (Pavoine et al. 2005), they were not transformed before calculating Q. With this approach, we preserve the definition of 'taxonomic diversity' (TD) as put forth by Shimatani v www.esajournals.org (2001), and note that this definition of TD is distinct from that used by Devictor et al. (2010) and Villéger et al. (2010), for example, who instead use 'taxonomic diversity' to refer to an index of species diversity (analogous to our Gini-Simpson index).

Statistical analysis
To explore the relationships among different indices of diversity, we calculated Spearman's rank correlation (r s ). This index tests for a monotonic relationship between two variables but does not assume bivariate normality (Quinn and Keough 2002), a condition which was lacking between some of our variables (see Fig. 1). Values of r s were deemed significantly different (i.e., there is no monotonic relationship between two variables) based on the test in Best and Roberts (1975).
To investigate phylogenetic and taxonomic signal in functional traits, we calculated matrix correlations between the functional, phylogenetic, and taxonomic distance matrices using Mantel's test. Mantel's test of matrix correlations removes the potentially confounding effect of richness and evenness on the relationships between FD, PD, and TD, and because it is based on distances, allows for the evaluation of signal across multiple functional traits simultaneously, including those that are not continuous (Pavoine et al. 2013).
To examine spatial and intra-annual trends in the different dimensions of demersal fish diversity, we first interpolated diversity for each month of the survey using ordinary kriging. This approach accounts for the spatial structure of the response, and has been used to address similar questions in the past (e.g., Devictor et al. 2010, Safi et al. 2011. All responses were log-transformed to better meet the assumptions of normality of errors and constant variance, except for evenness, which was arcsine-square root transformed. Predicted values were then backtransformed for plotting. In addition to interpolating the raw values, we also extracted and interpolated the residuals of each index fitted against each of the others using generalized additive models (GAMS, to account for nonlinearities, see below). The plotting of interpolated residuals allowed us to better visualize potential areas of mismatch between the different indices of diversity (e.g., Safi et al. 2011). To further clarify these mismatches, we plotted line graphs of the mean local diversity for each stratum of the sampling design, i.e., the average of all tows within a given region of the Bay (see legend of Fig. 3).
To understand the drivers of diversity patterns, we fitted generalized additive models (GAMs) to each index separately. GAMs provided a general and flexible modeling framework in which the effects of continuous covariates were modeled using non-parametric smoothing functions, and categorical factors were modeled parametrically to determine their mean effect sizes (Wood 2006, Zuur et al. 2009). We evaluated multiple explanatory variables and classified them into three general components of influence: temporal factors (encompassing year and month), spatial variables (latitude and longitude), and environmental covariates (salinity, water temperature, dissolved oxygen, and water depth). The fully parameterized GAM was therefore defined as: where y i is the calculated diversity for station i, a's are the estimated mean effects for each level of year (YEAR) and month (MONTH), and g's are nonparametric smoothing functions for each covariate, including salinity (SAL), temperature (TEMP), dissolved oxygen (DO), and water depth (DEPTH). We incorporated a separate longitude-latitude (LONG, LAT) smoother (g 1 ) for each month to account for any potential space-time interaction in the data following Buchheister et al. (2013). Thin plate regression splines with shrinkage terms were used as the basis to smooth all covariates. The intercept, a, scales the model prediction to the appropriate level of the response because each smooth estimate (g) is constrained to average to 0 over the entire dataset (Ciannelli et al. 2008). The residual error, e i , was assumed to be independent and identically distributed with a mean of zero and constant variance. No predictors were sufficiently collinear to warrant dropping from the analysis (see also Buchheister et al. 2013). This model was fitted to richness using a Poisson distribution and a log-link, while all other indices were log-transformed (except evenness which was arcsine-square root transformed) and fitted to a normal distribution. Graphical and statistical analyses indicated that differences in sampling effort (area swept) did not have discernible effects on diversity across stations and that all tows were sufficiently long to capture a representative sample of biological diversity. Therefore no offset was necessary to account for sampling effort in the GAMs (Zuur et al. 2009).
In order to assess the relative explanatory power of each component, we calculated partial deviances by sequentially removing suites of predictors from the full model (Eq. 4) corresponding to indicators of space, time, or environment, repeated this procedure for all possible permutations, averaged the deviances for all models in which a predictor appeared, and finally calculated a standard error (e.g., Schmiing et al. 2013). The partial deviances are therefore the proportion of total explained deviance of the model uniquely explained by space, time, or the v www.esajournals.org environment, accounting for the other predictors in the model. This approach is analogous to variance partitioning approaches in a linear modeling framework (Legendre and Legendre 1998), only extended to non-linear models. We performed all analyses in R (v.3.0.1, R Development Core Team 2013) using the following packages: ape (Paradis et al. 2004), automap (Hiemstra et al. 2009), clue (Hornik 2013), ggplot2 (Wickham 2009), and mgcv (Wood 2011).

Comparing dimensions of diversity
Across all tows and dates, we found a moderate positive association between species richness and other components of diversity (0.57 Spearman's rank correlation r s 0.62; Fig. 1). Evenness (computed from relative biomass) was strongly related to the other dimensions (0.69 r s 0.80, Fig. 1). Integrating both richness and evenness in the Gini-Simpson index, we found still stronger positive relationships with FD, PD, and TD (0.81 r s 0.91; Fig. 1). Using presenceabsence data revealed a slightly diminished but still strong relationship between the Gini-Simpson index and FD, PD, and TD (0.69 r s 0.87; Appendix B: Fig. B2), suggesting that the relationships between these dimensions were only marginally affected by the inclusion of relative biomass. Finally, FD, PD, and TD were nearly redundant in describing demersal fish diversity in this system (0.95 r s 0.97; Fig. 1), and weighting by presence-absence did not alter this conclusion (0.89 r s 0.92; Appendix B: Fig. B2). All values of Spearman's r s were significantly different from zero (P , 0.001). Mantel correlations confirmed that functional distances between species were strongly conserved in phylogenetic and taxonomic distances (0.74 Mantel's R 0.87). Consistent with our initial hypothesis, the strongest matrix correlation occurred between taxonomic and phylogenetic distances (R ¼ 0.87).

Spatial and seasonal patterns in diversity
Despite the moderate to strong correlations among the different dimensions of fish diversity across all samples, we did identify areas and times of incongruence. Richness was highest at the mouth of the Bay during July and September, decreased in the middle Bay, and increased again at the head of the Bay (Figs. 2A, 3A). This pattern of higher richness at the two extremes of the salinity gradient adhered to a well-established pattern along estuarine gradients (Odum 1988).
Evenness peaked early in the year (Figs. 2B, 3B), although this trend can be explained by the fact that most communities during March were composed of few species ( Fig. 2A). Evenness was lowest in the middle Bay during May and July (Figs. 2B, 3B) due to the overwhelming dominance of Atlantic croaker (Micropogonius undulatus), and increased in this region during September and November with the migration of M. undulatus out of the Bay (Figs. 2B, 3B). In the upper Bay, evenness decreased slightly from July to November, tracking the increasing dominance of white perch (Morone americana) and the concomitant migration of euryhaline fishes, such as spot (Leiostomus xanthurus) and weakfish (Cynoscion regalis), out of the Bay at the end of the summer (Figs. 2B, 3B).
The Gini-Simpson index, our measure of species diversity that incorporated both richness and evenness, peaked in September in the lower Bay and decreased in November (Figs. 2C, 3C). This peak was again largely attributable to the migration of the biomass-dominant species, M. undulatus, out of the Bay after July. While composition and distribution of biomass remained fairly consistent among the remaining species during this time, the loss of such a large amount of biomass with M. undulatus manifested as an increase in evenness. It is for this reason that the Gini-Simpson index was lower than would be expected from species richness in the lower Bay during July, and higher than would be expected in September based on residuals from this relationship (Fig. 4A).
Functional, phylogenetic, and taxonomic diversity showed largely identical patterns. All three of these dimensions peaked through November in the lower Bay, which is later in the year than for richness (July) or Gini-Simpson diversity (September, Figs. 2D-F, 3D-F). The later peak in FD, PD, and TD was attributable to the more even distribution of biomass among functionally and evolutionarily divergent species. In July, M. undulatus accounted for nearly 50% of the biomass in the lower Bay. By September, 35% of the biomass was dominated by two species in v www.esajournals.org the family Sciaenidae, M. undulatus and L. xanthurus, which are also very similar functionally (Appendix A: Fig. A2). By November, 60% of the biomass was dominated by L. xanthurus, summer flounder (Paralichthys dentatus), and clearnose skate (Raja eglanteria), all of which are more distant functionally and phylogenetically than the dominant species in earlier assemblages (Appendix A: Figs. A2, A3). To that end, weighting by presence-absence produced patterns that mirrored species richness (Appendix B: Fig. B7). Thus, incongruences between specieslevel indices of diversity and those that incorporate functional, phylogenetic, or taxonomic differences are largely driven by the distribution of biomass among divergent species. It is also important to note that the scale of FD, PD, and TD are all much lower than the Gini-Simpson index or richness, indicating that there was considerable functional, phylogenetic, and taxonomic redundancy within these assemblages, especially during the summer months (Figs. 2, 3).
There were some interesting but comparatively minor deviations among functional, phylogenetic, and taxonomic diversity. FD was higher than expected based on PD and TD during March at the mouth of the Bay (Fig. 4B, C). This mismatch was due to spiny dogfish (Squalus acanthias) being more distant functionally than evolutionarily from two co-occurring elasmobranchs, smooth dogfish (Mustelus canis) and clearnose skate (Raja eglanteria) (Appendix A: Figs. A2-A4). FD was also higher than expected from TD, and to a lesser degree PD, in the lower Bay during September and throughout the Bay during November (Fig. 4B, C). This pattern was due to the large contribution to biomass by two ray species in the same genus, Gymnura altavela and G. micrura, which are two-and five-times more distant functionally than they are phylogenetically and taxonomically, respectively (Appendix A: Figs. A2-A4). Finally, FD was lower than expected from both PD and TD during May and July at the head of the Bay (Fig. 4B, C), as the three dominant species in terms of biomass at those times (L. xanthurus, M. undulatus, and M. Fig. 3. Plot of mean local diversity (i.e., tow-level) for each region of the Chesapeake Bay. Regions corresponds to each of the strata used in the sampling design and capture the estuarine salinity gradient (see legend), with near marine salinity at the mouth and near freshwater at the head. The scales for panels A and C-F are in units of effective numbers of species; the scale for panel B is units ranging from 0 (maximally uneven) to 1 (maximally even). Gini-Simpson, functional, phylogenetic, and taxonomic diversity indices are weighted by relative biomass. Error bars are 61 SE.
v www.esajournals.org americana) are all closer functionally than they are phylogenetically or taxonomically (Appendix A: Phylogenetic diversity was less than expected from TD at the head of the Bay during May (Fig.  4D), due to the presence of three species of catfishes-brown bullhead (Ameiurus nebulosus), channel catfish (Ictalurus punctatus), and white catfish (Ictalurus catus)-which are more distant taxonomically from gizzard shad (Dorosoma cepedianum) than we recovered in our phylogeny (Appendix A: Figs. A3, A4). Finally, PD was greater than expected based on TD in the middle Bay during November (Fig. 4D). This outcome was due to the five dominant species (L. xanthurus, M. undulatus, M. americana, P. dentatus, and C. regalis) being more distant phylogenetically than taxonomically (Appendix A: Figs. A3, A4). The differences among FD, PD, and TD were exaggerated when calculated using presence-absence (Appendix B: Fig. B7), indicating that mismatches among these three aspects of Fig. 4. Kriging interpolation of the residuals from generalized additive models of (A) Gini-Simpson diversity against richness, (B) functional diversity against phylogenetic diversity, (C) functional diversity against taxonomic diversity, and (D) phylogenetic diversity against taxonomic diversity, for each month of the trawl survey and across all years. Warm colors indicate higher than excepted diversity, cool colors indicate lower than expected. All indices were weighted by relative biomass.
v www.esajournals.org diversity in particular were driven by functional and evolutionary differences among species, and not the distribution of biomass.
Patterns in functional, phylogenetic, and taxonomic diversity weighted by relative abundance (Appendix B: Fig. B4) were muted compared to those weighted by biomass (Fig. 2). This dampening was likely because the most functionally and phylogenetically unique species in our system were elasmobranchs, which are on average much larger than teleosts and also tended to occur in relatively low abundances. Thus, their contribution to diversity was upweighted when considering relative biomass, but down-weighted when considering relative abundance. A notable exception was FD and TD, which remained consistently high in the lower Bay throughout the year when weighted by relative abundance (Appendix B: Figs. B3, B5), as opposed to peaking in September and November when weighted by relative biomass (Figs. 2D, 3D). This result was due to tows with low abundance of M. undulatus that also contained species whose distribution was more even in terms of individuals than in terms of biomass: P. dentatus, R. eglanteria, windowpane flounder (Scophthalmus aquosus), and spotted hake (Urophycis regia). (Additional results corresponding to the different weighting schemes-biomass, abundance, and presence/absence-are presented in Appendix B.)

Drivers of observed patterns in diversity
Using generalized additive models, we were able to explain the largest proportion of deviance for richness (56%) and smallest proportion for evenness (9.9%). Roughly the same proportion of total deviance was explained for Gini-Simpson, functional, phylogenetic, and taxonomic diversity (33, 39, 40, and 39% respectively). We divided predictors into suites of variables representing space (latitude and longitude), time (year and month), and environment (salinity, temperature, dissolved oxygen, and depth), and used a deviance partitioning approach to identify the unique contributions of each suite of predictors to diversity. Space explained approximately 24%, time approximately 10%, and environment ap- Fig. 5. The partial contributions of space, time, and environment to the total explained deviance from generalized additive models fit to each diversity index. Space was a smoothed combination of latitude and longitude by month, time was a linear combination of month and year, and environment was a smoothed combination of salinity, temperature, dissolved oxygen, and depth. Units are in percentage of total explained deviance, and error bars denote 61 SE on partial deviance estimates (see text for calculations). Gini-Simpson, functional, phylogenetic, and taxonomic diversity indices are weighted by relative biomass.
v www.esajournals.org proximately 22% of the total explained deviance in richness (Fig. 5). These values were similar in magnitude for Gini-Simpson diversity, FD, PD, and TD, indicating that space, time, and environment have approximately equivalent influence across the various aspects of diversity (Fig.  5). The one exception was evenness, where space accounted for the largest proportion of the total explained deviance at 6.4%, and environment explained the least at ,1%. (Values for partial deviances corresponding to the different weighting schemes are given in Appendix B: Table B1.)

DISCUSSION
Across ten years of survey data on Chesapeake Bay demersal fishes, we found strong agreement between measures of functional, phylogenetic, and taxonomic diversity (Fig. 1), in contrast to what has been shown for birds, mammals, and other marine fishes (Devictor et al. 2010, Mouillot et al. 2011a, Safi et al. 2011. Our results were due to strongly conserved functional traits and evolutionary distances between species within this system. This finding also supported our original hypothesis that measures of diversity derived from taxonomic hierarchy are suitable proxies for those derived from molecular phylogenies (Ricotta et al. 2012).
Incorporating data on the functional traits and evolutionary relatedness of species shifted peaks in diversity to later in the year, as compared to those based on traditional metrics of species richness and diversity (Figs. 2, 3). This result was due to the more even distribution of biomass across functionally, phylogenetically, and taxonomically distinct species, supporting our prediction that well-documented patterns of species dominance in this system would modify the interpretation of FD, PD, and TD. Finally, we found that space, time, and environment explained approximately equivalent proportions of deviance across diversity indices (Fig. 5), suggesting that these extrinsic factors do not differentially affect the various aspects of community diversity.
While we were able to identify several specific regions and months where functional, phylogenetic, and taxonomic diversity differed (Fig. 4), overall these differences were not as striking as have been shown in previous studies (Devictor et al. 2010, Mouillot et al. 2011a, Safi et al. 2011. There are several possible reasons for this finding. First, our analysis was done on a smaller regional scale (10 5 km 2 ), whereas studies addressing similar questions have been conducted on much larger (10 7 km 2 ; Devictor et al. 2010, Mouillot et al. 2011a or even global scales (Safi et al. 2011). A previous study showed that 'hotspots' of avian species richness shifted in space depending on the scale at which the analysis was performed, which was a consequence of the rate at which species accumulated in areas of varying size (Lennon et al. 2001). Extending this idea to functional and phylogenetic diversity, increasing scale may encourage the more rapid accumulation of functionally or evolutionarily unique species. Thus, larger scales may increase the opportunity for mismatch between various dimensions of diversity (Devictor et al. 2010). However, we point out that the scale of our analysis is one that is relevant to management and conservation for this system, particularly within an ecosystem context (Latour et al. 2003).
Second, and related to the issue of scale, the number of species in our study was an order of magnitude lower than in similar investigations using French bird assemblages (S ¼ 229; Devictor et al. 2010) and Mediterranean fishes (S ¼ 282; Mouillot et al. 2011a), and two orders lower than for global mammalian assemblages (S ¼ 4536; Safi et al. 2011), although it is comparable to a study of FD in a community of tropical estuarine fishes (S ¼ 62; Villéger et al. 2010). Whereas scale may drive mismatch between different aspects of diversity through mechanisms like increased environmental heterogeneity leading to accumulation of species with contrasting phenotypes and/or evolutionary histories, a larger species pool may do so purely through a sampling effect. In other words, having more species increases the probability that an assemblage will contain a functionally or phylogenetically unique species. This idea argues for a consideration of both scale and the size of the species pool, which are not necessarily independent, when characterizing diversity. However, values of species richness in the Chesapeake Bay are comparable to similar estuarine systems across North America and Europe (Lotze et al. 2006), and other coastal ecosystems are facing threats to local diversity from pollution, climate change, and overfishing (Sala and Knowlton 2006). Our analysis provides an example of how different dimensions of diversity might behave in similarly depauperate or increasingly species-poor systems.
Beyond issues of scale, the relatively low diversity of our study stems from several additional sources. First, the Bay can be considered a relatively harsh environment in that much of its volume is brackish, and the abiotic environment fluctuates strongly on seasonal or shorter timescales. Many fishes are not adapted to endure the osmoregulatory demands of fluctuating salinities (Day et al. 2012), and this physiological stress imposes strong spatial restrictions on local species richness. This effect is most evident in the lack of overlap in species composition between the near freshwater head of the Bay and the near marine mouth (with the exception of a few cosmopolitan species). Adaptation to varying levels of stress may also be partly responsible for the strong evolutionary conservation of traits revealed by the high Mantel correlations between our functional and phylogenetic matrices. Specifically, 9 of the 50 species in our data set belong to a single family, Sciaenidae, and members typically made up around a quarter to half of the biomass at any given station. Second, many species reside seasonally in the Bay, principally to feed, before migrating to overwintering grounds in the coastal waters of the Atlantic Ocean (Murdy and Musick 2013). As a result, local diversity in the Bay is limited temporally by active dispersal and life history. Finally, we restricted our analyses to demersal fishes to minimize potential bias caused by gear-selectivity (e.g., Villéger et al. 2010). Thus, our findings do not include pelagic or vagrant fishes, but the standardized approach of the survey sampled a large proportion of the dominant species in the Bay (Jung andHoude 2003, Buchheister et al. 2013). Consequently, our results are representative of spatial and temporal patterns in demersal fish diversity in the Chesapeake Bay.
Our analysis also provides strong evidence for the use of taxonomic distances as a proxy for functional and phylogenetic differences among species. While a recent study has proposed the use of taxonomy as a surrogate for phylogeny, they did so using plant assemblages numbering in the hundreds of species (Ricotta et al. 2012).
We demonstrate that this substitution is valid for assemblages with many fewer species. The strong relationship between TD and PD in our study is likely due to the fact that species were taxonomically diffuse, and thus evolutionary distances between species were sufficiently broad to be accurately captured by taxonomy. Had our assemblages been more closely related (e.g., within the same genus), then it is possible that taxonomic hierarchy would be too coarse to delineate the fine-scale differences revealed by molecular methods. Furthermore, our (gene tree proxy for) molecular phylogeny was very close to existing taxonomy. Phylogenetic hypotheses are constantly shifting in light of new (molecular) evidence, and taxonomy is slower to reflect these changes. In our case, evolutionary relationships suggested by phylogeny were already welldelineated by taxonomy, but in cases where phylogenetic evidence conflicts with existing information, it is up to the investigator whether to proceed with potentially dated taxonomic information, or with newer phylogenetic evidence.
Nevertheless, the use of taxonomic hierarchy over functional traits or molecular phylogeny to characterize diversity has important practical implications. Taxonomic data are easily accessible, and calculations of taxonomic diversity are relatively straightforward compared to constructing and extracting indices from molecular phylogenies. Furthermore, while molecular methods are improving in cost and efficiency, there will presumably always be organisms for which robust phylogenies are missing (Crozier et al. 2005, Ricotta et al. 2012. As a result, a common practice is to use taxonomy to resolve incomplete phylogenies (Wiens 2004). In such cases, and in cases where morphological taxonomists have delineated evolutionary differences reflecting phylogenetic signal among species, taxonomy and phylogeny are likely to lead to similar outcomes. Thus, we recommend using taxonomy when robust phylogenetic data are simply not available, when assemblages are taxonomically diffuse, and when taxonomically inferred evolutionary relationships are well-supported.
Our analysis also explicitly shows that the distribution of biomass among functionally or evolutionarily unique species can alter seasonal diversity patterns, specifically in sustaining high v www.esajournals.org FD, PD, and TD after peaks in species richness and Gini-Simpson diversity have tapered off (Figs. 2, 3; Appendix B). This finding supports the conclusions from previous analyses of groundfish diversity which showed that richness, evenness, and taxonomic differences provide unique information about community structure in the Mediterranean Sea (Mérigot et al. 2007a, 2007b). Our investigation emphasizes the interchangeability of functional or phylogenetic information with taxonomic data within this descriptive framework.
The relationship between evenness and functional and phylogenetic identity may also have implications for ecosystem processes. Theory predicts that the more even distribution of individuals or biomass among functionally complementarity species should enhance ecosystem functioning, although experimental evidence, mostly from grassland plants, is conflicted (Hillebrand et al. 2008). Moreover, recent evidence from tropical reefs and forest systems suggests that less abundant species can have traits that support unique or vulnerable functions (Mouillot et al. 2013). Reconciling the ecosystem implications of abundance-or biomass-distribution among species based on their traits or evolutionary history is a key issue moving forward.
Our use of generalized additive models revealed relatively consistent effects of space, environment, and time across different dimensions of diversity, suggesting that these three components do not independently affect different aspects of the Chesapeake Bay demersal fish community. However, we did find a larger role for factors relating to space and environment versus those relating to time (Fig. 5). This result once again emphasizes the role of spatially conserved environmental gradients, principally salinity, in driving diversity patterns in the Bay (Buchheister et al. 2013). It also suggests that temporally dependent migration and dispersal may have a smaller influence on the same patterns. While most fishes included in our study migrate in and out of the Bay at some point during their lives, their location within the Bay (and thus their contribution to local diversity) is still largely restricted by their physiological tolerances and preferences (Wagner 1999). The strong influence of environment in mediating dispersal, and thus spatial patterns in communi-ty diversity, is consistent with findings across many other systems (Cottenie 2005). This result is also consistent with the predictions by Stirling and Wilsey (2001) who posit that diversity indices, such as richness and evenness, will be strongly correlated when dispersal and migration are structuring communities, as is the case in the Chesapeake Bay. However, our analysis also indicated that a large proportion of deviance in diversity was independently explained by environment, suggesting that a role for non-spatially structured environmental variation exists in driving community composition in the Bay. One possible explanation lies with the extent of areas of low dissolved oxygen, which are not as spatially conserved as salinity (although do tend to occur in deep channel of the middle Bay), but which have a strong documented effect on local diversity (Buchheister et al. 2013).
Functional and phylogenetic diversity are increasingly applied in the ecological and conservation literature (Cadotte et al. 2011), yet our results suggest that their usefulness varies among systems. We show for Chesapeake Bay demersal fishes that measures of diversity incorporating taxonomic distances reflect spatial and temporal patterns as well as those that incorporate functional trait or phylogenetic information, and thus may provide a convenient and reliable proxy when phylogenetic data are redundant or unavailable. We also show that the distribution of biomass among functionally and evolutionarily divergent species can significantly alter patterns of diversity. Characterizing large-scale patterns in diversity, especially in relatively species poor systems, will be potentially more accurate and useful when it combines information on species composition, relative abundance or biomass, and existing taxonomic information before moving on to more intensive collection of functional trait and phylogenetic data.

ACKNOWLEDGMENTS
We thank all current and past ChesMMAP personnel, especially Christopher Bonzek, James Gartland, and the crew of the R/V Bay Eagle, for their diligent work in collecting and maintaining all the survey data used in this study. Support for ChesMMAP was provided by the NOAA Chesapeake Bay Office, the Virginia Environmental Endowment, the U.S. Fish and Wildlife Service, and the Virginia Marine Resources v www.esajournals.org Commission. We thank the staff at the Smithsonian Museum Support Center, especially Jeff Williams and Jerry Finan, and the VIMS Ichthyology Collection for allowing us to measure specimens. We also thank two anonymous reviewers for their comments. This work was supported by the National Science

Trait selection and measurement
We assembled information on 28 functional traits relating to biomass accumulation, distribution, food acquisition, locomotion, and reproduction for the 50 most abundant fish species in the ChesMMAP survey during the study period (Table A1). The 28 traits, their functional interpretation, and selected references are provided in Table A2. After examining relationships among continuous traits using Spearman rank correlations, we omitted three traits (indicated with in Table A2) that were strongly correlated and represented similar ecological functions to other traits (e.g., mean and maximum biomass).
Fourteen morphoanatomical measurements (Fig. A1, Table A3) were used to derive most ecomorphological traits (see legend, Table A2). The selection of ecomorphological traits was guided largely by the work of Villéger et al. (2010) and references therein, but also by Gatz Jr. (1979a) and Mason et al. (2008). Most morphoanatomical measurements were made directly from approximately 5 individuals of each species (range: 1-7) on museum fish specimens at the Virginia Institute of Marine Science Fish Collection and the Smithsonian Museum of Natural v www.esajournals.org History Fish Collection using calipers with 0.1 mm precision. Exceptions include caudal fin surface area, which was measured digitally from scientific illustrations depicted in Murdy et al. (2002), and caudal fin height, which was measured on both museum specimens and from illustrations in Murdy et al. (2002). Because most morphoanatomical measurements were made on preserved museum specimens, measurements requiring destructive manipulation of mouth or gill parts were considered, but ultimately not conducted.

Functional dendrogram
To calculate functional diversity, we constructed a functional dendrogram using hierarchical agglomerative clustering (Petchey and Gaston 2002), which groups species based on their   (Gatz 1979b) v www.esajournals.org  et al. 1997 positions in multivariate trait space. Because several different algorithms exist for hierarchical agglomerative clustering (Legendre and Legendre 1998), we used the approach outlined by Mouchet et al. (2008) of building multiple dendrograms using different algorithms, and choosing the one that best preserves the original dissimilarities between species. First, we computed a dissimilarity matrix between each pair of species from the functional trait data using Gower's dissimilarity index (Gower 1971). The Gower index is an ideal distance measure because it can incorporate continuous, categorical, and ordinal trait information (after applying the correction from Podani 1999), account for missing trait values, and apply weighting to individual traits. Our trait database represents all three types of data, and contains some missing values (e.g., the absence of morphological characters relating to the caudal fin in batoids, or skates and rays). Additionally, we chose to downweight the five traits representing different aspects of diet (percent crustaceans, percent fish, etc.), because these traits are not independent. Knowing the percent of the diet composed of molluscs provides some information about the percent of the diet composed of worms; for instance: % worms 100 À % molluscs. As a result, we weighted each diet trait by 1/5 so that collectively, they had equal weight to each of the other 24 traits in the trait database. Gower dissimilarities were calculated using the function gowdis in the package FD (Laliberté and Shipley 2011).
Second, we constructed functional dendrograms from the Gower dissimilarity matrix using  Nammack et al. 1985;Reichert 1998;Ross et al. 1995;Springer 1960;Struhsaker 1969;Waring 1984;White et al. 2003 Spawning location River, estuary, ocean Spatial spawning distribution and where young are likely to recruit The area in which spawning occurs, not necessarily within Chesapeake Bay. Daiber and Booth 1960;FLMNH 2012;García-Abad et al. 1999;Hildebrand and Cable 1938;Murdy et al. 1997;Packer et al. 2003;Stokes 1980;Wourms 1977 Daiber and Booth 1960;FLMNH 2012;García-Abad et al. 1999;Hildebrand and Cable 1938;Murdy et al. 1997;Packer et al. 2003;Stokes 1980;Wourms 1977 Notes: Eleven traits were considered ecomorphological (indicated with *), nine of which and were calculated from morphoanatomical measurements described in Table A3 (''dentition'' and ''swimming mode'' were scored as categorical and thus did not require measurements). Three traits were dropped from the final analysis due to high Spearman rank correlations (indicated with ). Traits were directly measured if source is blank. Full references in the 'Source' column are provided below.
§ Other primary sources are Bowman et al. 2000, Chao and Musick 1977, Cruz-Escalona et al. 2005, Ellis and Musick 2007, Hansen 1969, Hettler 1989, Johnson and Seaman 1986, Link et al. 2002, Martins et al. 2005, Ross 1989, Scharf and Schlicht 2000, Steimle and Shaheen 1999 v www.esajournals.org the following agglomeration methods: single linkage, complete linkage, Ward's minimum variance, unweighted arithmetic average clustering (UPGMA), and weighted arithmetic average clustering (WPGMA). We omitted two other algorithms suggested by Mouchet et al. (2008)-unweighted centroid clustering (UPGMC) and weighted pair group centroid clustering (WPGMC)-because they are not appropriate for non-metric distances. When used with Gower dissimilarities, these two methods can lead to reversals (or inversions) in the dendrogram that violate the ultrametric property (Legendre and Legendre 1998), and we required the dendrogram to be ultrametric (Pavoine et al. 2005). Dendrograms were built using the function hclust in the R base package. We also built a consensus dendrogram, which utilized all four clustering algorithms to produce a dendrogram that emphasized areas of agreement between the different algorithms (Mouchet et al. 2008). The consensus dendrogram was built using the function cl_consensus in the clue package (Hornik 2013).
Third, we computed the ultrametric distances from each dendrogram, in which the root-to-tip distance of the dendrogram was scaled to 1 using the function cl_ultrametric in the clue package (Hornik 2013). We compared the original distances to the ultrametric distances using a matrix norm, the 2-norm . We then extracted the pairwise distances from the ultrametric dendrogram that best preserved the original distances (i.e., lowest 2-norm value), which was the UPGMA dendrogram, using the function cl_dissimilarity in the clue package (Hornik 2013). Finally, we standardized the values so the maximum distance between species equaled 1 (Devictor et al. 2010). The final functional dendrogram is presented in Fig. A2 and is also provided as a Newick-format file.

Phylogenetic tree
To maximize data coverage, we queried GenBank to locate sequences available for vouchered specimens included in this study (Table A1). We found that cytochrome oxidase I (COI ) was available for 34 of our 50 target species. Vouchered tissue samples were available in the Virginia Institute of Marine Science Ichthyology  Table A2.
v www.esajournals.org Collection for three of the missing species (Urophycis regia, Astroscopus guttatus, and Etropus microstomus). For the 13 species for which both tissue samples and sequence data were unavailable, we substituted COI sequences for closely related congeners (indicated with an asterisk in Table A4). Only one congener was unavailable for Larimus fasciatus, in which case we substituted Cynoscion nebulosus based on the phylogeny of Sasaki (1989). In addition to the 50 species, we obtained 5 sequences to serve as outgroup taxa (indicated with a dagger in Table A4). The outgroup taxa were a mixture of Myxiniformes and Pertomyzontiformes, which are basal to Chondrichthyes, and included: Caspiomyzon wagneri, Eudontomyzon lanceolata, Ichthyomyzon bdelli-  We extracted total genomic DNA from vouchered tissue plugs for U. regia, A. guttatus, and E. microstomus using QIAGEN's DNAEasy tissue kit. We completed PCR amplification using universal primers given in (Folmer et al. 1994). PCR conditions were as follows: an initial 4 minute denaturation at 948C, followed by 45 cycles of 948C for 1 min, 458C for 1 min, and 728C for 1 min, with a final extension at 728C for 4 min. We cleaned up successfully amplified products using the QIAquick PCR Purification kit and eluted to a final volume of 25 ll. We used an ABI Sequencing Reaction Protocol and sequenced amplified COI in both directions on an Applied Biosystems 3130xl ABI Prism Genetic Analyzer   2006). This program searched for the tree with the best likelihood scores from 200 random maximum parsimony start trees, with GAMMAþP-Invar model of rate heterogeneity and the ML estimate of the alpha-parameter. Node supports were obtained from 999 bootstrapped trees using the rapid hill-climbing algorithm. The tree was rooted using the five outgroups, and we used the gene tree as a proxy for phylogenetic distance among all non-outgroup species.
To convert the tree to ultrametric, we used the function chronopl in the ape package (Paradis et al. 2004), which estimated divergence times using the penalized likelihood algorithm; lambda (the smoothing parameter) was determined through a data-driven cross-validation procedure (Sanderson 2002). We rooted the ultrametric tree, pruned the outgroups, extracted the cophenetic distance matrix from the ultrametric tree, and standardized the values so the maximum distance between species equaled 1. The final phylogenet-ic tree is presented in Fig. A3 and is also provided as a Newick-format file.

Taxonomic tree
We first obtained classification information for all 50 species using the latest accepted taxonomy from the Integrated Taxonomic Information System (http://www.itis.gov), which are presented in Table A1. The taxonomic tree was constructed from this classification table using the function as.phylo in the ape package (Paradis et al. 2004). Ultrametric distances were computed as the number of nodes separating each pair of species on the hierarchy (Clarke and Warwick 1998), weighted by the proportional decrease in taxon richness at each successive level on the hierarchy (Clarke and Warwick 1999). We obtained the distances using the function taxa2dist with the argument varstep ¼ TRUE in the vegan package (Oksanen et al. 2013). The final taxonomic tree is presented in Fig. A4 and is also provided as a Newick-format file.

Expanded analyses using abundance-and presence/absence-weighted indices
We generated scatterplot matrices to explore the associations between diversity indices-richness, evenness, Gini-Simpson, functional, phylogenetic, and taxonomic diversity-weighted by relative abundance (Fig. B1) and presence/absence (Fig. B2). There were overall weaker associations between components of diversity weighted by relative abundance (Fig. B1) compared to relative biomass (Fig. 1, main text), as indicated by Spearman's rank correlation coefficient (r s ). As noted in the main text, this result was likely due to rare but large species driving stronger dominance in terms of biomass than in terms of abundance. However, we emphasize that functional, phylogenetic, and taxonomic indices were all still highly associated (0.92 r s 0.97; Fig. B1), with the strongest relationship again occurring between phylogenetic and taxonomic diversity (r s ¼ 0.97). We also note that, when weighting by presence/absence, the Gini-Simpson index collapses to species richness, hence the perfect correlation (Fig. B2).
We mapped each diversity index weighted by relative abundance (Fig. B3) and presence/absence (Fig. B4), and constructed corresponding line graphs of mean local diversity by month and region of the Bay (Figs. B5, B6). As noted in the main text, indices weighted by abundance show largely congruent patterns with richness, with the exception of abundance-weighted functional diversity peaking in March as well as November, and taxonomic diversity remaining relatively high throughout the year in the lower Bay (Figs. B3, B5). Again, weighting by presence/absence showed patterns nearly identical to richness (Fig.  B4), indicating that mismatches identified in the main text were the function of uneven distribution of biomass (Fig. 2, main text). We also mapped the residuals of all pairwise combinations of each index against the other, fitted using generalized additive models (GAMs). Combinations of biomass-weighted indices are given in Figs. B7-B9; we note a subset of this figure is represented in Fig. 4 in the main text. Abundance-weighted indices are given in Figs. B10-B12, and presence/absence-weighted indices are given in Figs. B13-B15.
Finally, we applied a deviance partitioning approach to determine the relative contributions of space, time, and the environment in driving patterns in diversity. Partial deviances and their standard errors are given for all indices and weighting approaches in Table B1. We also provide barplots of the partial deviances for abundance-weighted (Fig. B16) and presence/ absence-weighted indices (Fig. B17) for comparisons with Fig. 5 in the main text.               Table B1. Partial deviances for each diversity index fitted against indicators of space (latitude and longitude), time (month and year), and environment (temperature, salinity, dissolved oxygen, and depth) using generalized additive models. Values are percentages of the total explained deviance 6 SE.

SUPPLEMENT
R script containing all data analyses and functional, phylogenetic, and taxonomic trees in Newick format (Ecological Archives C005-001-S1). Fig. B16. The partial contributions of space, time, and environment to the total explained deviance from generalized additive models fit to each diversity index. Gini-Simpson, functional, phylogenetic, and taxonomic diversity indices are weighted by relative abundance. Error bars denote 61 SE on partial deviance estimates. Interpretation given in the legend of Fig. 5, main text. Fig. B17. The partial contributions of space, time, and environment to the total explained deviance from generalized additive models fit to each diversity index. Gini-Simpson, functional, phylogenetic, and taxonomic diversity indices are weighted by presence/absence. A model could not be fit to evenness since the index collapses to 1 using presence/absence data. Additionally, the Gini-Simpson index collapses to richness when using presence/absence data and is likewise excluded. Interpretation given in the legend of Fig. 5, main text.