Species Diversity Associated with Foundation Species in Temperate and Tropical Forests

Aaron M. Ellison 1,* , Hannah L. Buckley 2 , Bradley S. Case 2, Dairon Cardenas 3, Álvaro J. Duque 4, James A. Lutz 5 , Jonathan A. Myers 6 , David A. Orwig 1 and Jess K. Zimmerman 7 1 Harvard University, Harvard Forest, Petersham, MA 01366, USA; orwig@fas.harvard.edu 2 School of Science, Auckland University of Technology, Auckland Central 1010, New Zealand; Hannah.Buckley@aut.ac.nz (H.L.B.); Bradley.Case@aut.ac.nz (B.S.C.) 3 Instituto Amazónico de Investigaciones Científicas Sinchi, Leticia, Colombia; dcardenas@sinchi.org.co 4 Universidad Nacional de Colombia Sede Medellín, Medelliín, Colombia; ajduque@unal.edu.co 5 Utah State University, S.J. & Jessie E. Quinney College of Natural Resources and the Ecology Center, Logan, UT 84322, USA; james.lutz@usu.edu 6 Washington University in St. Louis, Department of Biology, Saint Louis, MO 63130, USA; jamyers@wustl.edu 7 University of Puerto Rico – Rio Piedras, Department of Environmental Sciences, Rio Piedras 00925, Puerto Rico; jesskz@ites.upr.edu * Correspondence: aellison@fas.harvard.edu; Tel.: +1-978-756-6178


Introduction
Foundation species (sensu [1,2]) define and structure ecological communities and entire ecosystems through bottom-up control of species diversity and non-trophic modulation of energy and nutrient cycles [3]. Foundation species tend to be common and abundant and generally receive less attention from ecologists, regulatory agencies, or conservation biologists who are otherwise focused on the study, management or protection of rare, threatened, or endangered species [4,5]. However, because foundation species are likely to control the distribution and abundance of such rare species, it has been argued that foundation species should be protected before their populations decline to non-functional levels or disappear entirely [6].  The three temperate plots are all in the USA, and include Wind River in southwest Washington State, Harvard Forest in Massachusetts, and Tyson Research Center in Missouri ( Figure 1). The 25.6-ha Wind River plot ("WR"; 45.82 • N) is located in a late-seral Pseudotsuga menziesii/Tsuga heterophylla forest of the T.T. Munger Research Natural Area of the Gifford Pinchot National Forest. WR is at ≈365 m a.s.l. and receives ≈2500 mm of precipitation yr −1 , much of which falls as snow, leading to a 2-3 month summer dry season [26]. Half of the aboveground biomass at WFDP is concentrated in trees ≥ 93 cm DBH [27]. Regeneration of most species occurs preferentially in gaps created by windthrow [28], and large-diameter trees strongly inhibit recruitment [29].
The 35-ha Harvard Forest plot ("HF"; 42.54 • N) is a mix of deciduous and evergreen trees located at the Harvard Forest Long Term Ecological Research Site within the Worcester/Monadnock Plateau ecoregion [30] of transition hardwood/Pinus strobus/Tsuga canadensis forests [31]. HF is at ≈350 m a.s.l. and receives ≈1150 mm of precipitation (as rain and snow) yr −1 [32]. The HF plot is part of the Prospect Hill tract at Harvard Forest and has a mixed history of land-use and disturbance that includes selective logging and partial clearing for pasturage through the late 1800s, various outbreaks of nonnative pathogens and insects in the 1900s, and periodic hurricanes (in 1815 and 1938) and ice storms (most recently in 2008) [33][34][35].
The 20-ha Tyson plot ("TY"; 38.52 • N) is a late-successional temperate deciduous oak-hickory forest located at Washington University in St. Louis' Tyson Research Center on the northeast edge of the Ozark plateau. TY is at ≈200 m a.s.l. (elevational range 172-233 m a.s.l.), receives ≈ 1000 mm precipitation (as rain and snow) yr −1 , and contains strong edaphic and topographic gradients typical of the Missouri Ozarks [36,37]. Disturbance history at TY includes moderate grazing and selective logging in the early 1900s [38] and extreme droughts in 1988 and 2012. TY has remained relatively undisturbed by humans for ≈80 years [38].
The three tropical plots are in the Caribbean, Central America, and South America. The 16-ha Luquillo plot ("LFDP": 18.20 • N) is located within the University of Puerto Rico's El Verde Research Area of the Luquillo Experimental Forest on the Caribbean island of Puerto Rico. This wet tropical forest (≈3500 mm rainfall yr −1 ) includes a mix of well-conserved and secondary stands (logged and used for subsistence agriculture approximately a century ago) located at ≈380 m a.s.l. [39]. Hurricanes strike the region on average ≈ every 50 years. The most recent large hurricanes (Category 2 or higher) that affected the plot before the census data analyzed here were Hurricanes Hugo (1989) and Georges (1998) [40]. These caused defoliation and branch breakage but relatively little damage to stems, resulting in widespread recruitment of shrubs and sapling trees followed by subsequent thinning as the canopy re-establishes [41].
The 50-ha plot on Barro Colorado island, Panama ("BCI"; 9.15 • N) is on a hilltop that became an island when the Panama Canal was filled and Gatun Lake was created in 1914; the associated Barro Colorado Nature Monument has been managed since 1923 by the Smithsonian Institution. BCI is at ≈140 m a.s.l. and receives ≈2500 mm rainfall yr −1 ; five-six months of the year make up a pronounced dry season that receives little rain [42,43]. Hurricanes have never been recorded at BCI, and gaps created by single-to-multiple treefalls are the primary type of disturbance observed there [44]. At scales from individual gaps to the entire BCI plot, species diversity is more strongly controlled by recruitment limitation than gap dynamics [44].
The 25-ha Amacayacu plot ("AM": 3.8 • S) is located within the Amacayacu National Natural Park in Colombia, near the joint border of Colombia, Peru, and Brazil. AM is at ≈93 m a.s.l. within a wet tropical (≈3200 mm rainfall −1 ) terra-firme forest in the Colombian Amazon [45]. Although topographic relief at Amacayacu varies only over 20 m, it is more tightly associated with woody plant distributions than soil chemistry alone or in combination with topography [45].

Identifying Candidate Foundation Species
At HF, four trees numerically dominate the assemblage: Acer rubrum, Quercus rubra, Pinus strobus and Tsuga canadensis. Decades of observational work throughout New England and experimental work at Harvard Forest and elsewhere have consistently supported the hypothesis that T. canadensis is a foundation species whereas the other three are not [16,17]. We started, therefore, by asking whether commonly used graphical assessments of forest structure could distinguish T. canadensis from the other three dominant species. Specifically, we first plotted a basic size-frequency plot for the entire 35-ha Harvard Forest plot; we used DBH of the living woody stems as the measure of tree size and plotted the species-specific mean DBH against the total number of living stems of each of the 51 living woody species in the plot ( Figure 2). This size-frequency plot showed a typical "reverse-J" shape [46], but the aforementioned four dominant species were "outliers", well to the right of the "reverse-J". To explore the size-frequency relationship of these four species in a more fine-grained way-at the scale at which interactions among individual trees is stronger-we rasterized the plot with a 20-m grid [25]. We then plotted the total basal area of each of the four dominant tree species in each of the 20 × 20-m contiguous subplots within HF ( Figure 3). This diameter-abundance plot showed that the foundation species T. canadensis dominated each subplot both numerically and in total basal area, whereas in any given subplot, the other three dominant species either had many individuals or large total basal area, but not both.  Table 1) fall well away from the expected "reverse-J" distribution.   Table 1) fall well away from the expected "reverse-J" distribution.  Table 1) fall well away from the expected "reverse-J" distribution.  Table 1.  Table 1. We then used these two graphical "fingerprints" to identify dominant species (in terms of abundance and size) in the other five plots (Figures 4 and 5). This analysis identified distinctive species in terms of departure from the "reverse-J" size-frequency distribution at the three temperate plots (Tsuga heterophylla and possibly Pseudotusga menziesii and the subcanopy tree Acer circinatum at WR; three Quercus spp. at TY), and at two of the tropical plots (the canopy tree Dacryodes excelsa and the understory palm Prestoea acuminata at LFDP; and Eschweilera coriacea at Amacayacu) ( Figure 4). Other than P. menziesii, A. circinatum, and E coriacea, these species also had basal area-abundance distributions similar to that seen for Tsuga canadensis at HF ( Figure 5). For the tropical sites, we also included in subsequent analyses species identified by the site PIs, who established these forest dynamics plots and are co-authors of this paper, as being of known importance for forest dynamics ( Table 1). Among these, Oenocarpus mapora at BCI had a basal area-abundance distribution similar to that seen for T. canadensis at HF ( Figure 5), but it fell within the "reverse-J" cloud of points in the diameter-abundance distribution ( Figure 4). We also included in our analysis several additional species at WR and TY that were identified by the site PIs as ecologically important. These included the canopy species Abies amabilis at WR and three Carya species at TY; and the understory or lower canopy species Taxus brevifolia at WR; and Asimina triloba, Cornus florida, Lindera benzoin, and Frangula caroliniana at TY. We then used these two graphical "fingerprints" to identify dominant species (in terms of abundance and size) in the other five plots (Figures 4 and 5). This analysis identified distinctive species in terms of departure from the "reverse-J" size-frequency distribution at the three temperate plots (Tsuga heterophylla and possibly Pseudotusga menziesii and the subcanopy tree Acer circinatum at WR; three Quercus spp. at TY), and at two of the tropical plots (the canopy tree Dacryodes excelsa and the understory palm Prestoea acuminata at LFDP; and Eschweilera coriacea at Amacayacu) ( Figure 4). Other than P. menziesii, A. circinatum, and E coriacea, these species also had basal area-abundance distributions similar to that seen for Tsuga canadensis at HF ( Figure 5). For the tropical sites, we also included in subsequent analyses species identified by the site PIs, who established these forest dynamics plots and are co-authors of this paper, as being of known importance for forest dynamics ( Table 1). Among these, Oenocarpus mapora at BCI had a basal area-abundance distribution similar to that seen for T. canadensis at HF ( Figure 5), but it fell within the "reverse-J" cloud of points in the diameter-abundance distribution ( Figure 4). We also included in our analysis several additional species at WR and TY that were identified by the site PIs as ecologically important. These included the canopy species Abies amabilis at WR and three Carya species at TY; and the understory or lower canopy species Taxus brevifolia at WR; and Asimina triloba, Cornus florida, Lindera benzoin, and Frangula caroliniana at TY.     Figure 4. Size-frequency distributions of the species in each of the six studied ForestGeo plots. Panels are ordered from the northernmost temperate (top left) to equatorial (bottom right). Species that do not lie on the "reverse-J" line or that were otherwise are thought to be "important" species (in the tropical plots) are identified in red. Species abbreviations are given in Table 1.

Figure 4.
Size-frequency distributions of the species in each of the six studied ForestGeo plots. Panels are ordered from the northernmost temperate (top left) to equatorial (bottom right). Species that do not lie on the "reverse-J" line or that were otherwise are thought to be "important" species (in the tropical plots) are identified in red. Species abbreviations are given in Table 1.  Figure 4 as "important" in the six studied ForestGeo plots. Plots are ordered from the northernmost temperate (top left) to equatorial (bottom right). Species abbreviations are given in Table 1.  Figure 4 as "important" in the six studied ForestGeo plots. Plots are ordered from the northernmost temperate (top left) to equatorial (bottom right). Species abbreviations are given in Table 1.

Metrics of Forest Structure and Species Diversity
For each of the focal species (Table 1), we calculated the total number of stems, the total basal area and the mean basal area at both the plot scale (Table 1) and within cells of the 20-m raster (i.e., each of the contiguous 20 × 20-m subplots within each plot). In each subplot, for the remainder of the woody plant community (i.e., all species other than the focal species), we calculated the total abundance of stems, the total number of species, the Hill-number equivalents of Shannon's diversity and the inverse Simpson's indices [47], and the mean Bray-Curtis dissimilarity (of each subplot relative to each of the other subplots in the raster [48]; the latter is a measure of how "modal" or "outlying" each subplot is in its species composition relative to all other grid cells. For the calculation of Bray-Curtis dissimilarity, empty grid cells (i.e., those with no species other than the focal species) were assigned a dummy-value ("empty") to allow the calculation of all pairwise dissimilarities; this means that all empty cells appear to be maximally different from occupied cells, but maximally similar to each other. For all datasets, we used only the records for the main stem of each species, thus removing any additional stems; this eliminated the problem in spatial analyses of having multiple stems at the same point location. Initial exploration of the data showed that the effect of removing additional stems (primarily of small shrubs) on relative patterns in stem abundance and basal area was negligible for most species across all datasets. The diversity() and vegdist() functions in the R vegan package [49] were used for calculating each diversity metric.

Codispersion Analysis
Codispersion analysis quantifies the spatial association between two variables measured at individual spatial locations or on a grid [19]. Here, we examined codispersion between the total basal area of each focal species in each subplot of the 20-m raster versus each of the community metrics summarizing the diversity and composition of the other species in same subplot [20]. The spatial association quantified by the codispersion coefficient ranges from −1 to 1, where positive values represent a positive spatial association and negative values represent a negative spatial association. The results are presented as a hemispherical codispersion graph, which shows how the coefficient changes with distance (i.e., spatial lags) and direction (angles around a hemisphere). Thus, differences observed at different spatial lags indicate the scale(s) over which spatial processes are operating within the forest plots, whereas differences in the strength or sign of the correlation in different directions are diagnostic of anisotropic spatial processes [20]. To ensure sufficient sample size when computing codispersion coefficients between measurements made in 20-m grid cells (the size of each subplot), we used spatial lags in 20-m intervals from 20 to one-quarter the size of the minimum dimension of each plot [50]. We calculated codispersion coefficients using code custom-written and compiled in C (to reduce computation time), but with a link to R that allows for easy manipulation of input and output datasets [51].

Significance Testing with Null Model Analysis
For each dataset, focal species, and community metric combination, we recalculated the codispersion 199 times using the observed focal species total basal area raster and both (1) community metric rasters calculated from the point-pattern of the other species under a random shift around a torus ("toroidal-shift" null model) and (2) a spatially randomized raster of the community metric for the other species in the community ("CSR" null model) [21]. For each null model, the observed codispersion value for each cell was then compared to the distribution of the 199 null values and deemed significant if it fell outside 95 percent of the values (i.e., a two-tailed test). The toroidal-shift null model breaks the spatial association between the focal species and the other species in the community but retains any larger-scale spatial patterns, making this test slightly more conservative because it accounts for spatial patterns caused by the environment. In contrast, the CSR null model breaks the spatial association between the focal species and the other species in the community and simultaneously randomizes any larger-scale spatial patterns caused by environmental patchiness in the plots.

Data and Code Availability
Each of the ForestGeo plots were established at different times and have been censused every five years. To maximize comparability among datasets, we used data collected within a fifteen year period: the first censuses at Wind River All datasets are available from the ForestGEO data portal https://ctfs.si.edu/datarequest). R code for all analyses is available on GitHub (https://github.com/buckleyhannah/FS_diversity).

Forest Structure and Species Diversity
Across the entire plots and within the 20 × 20-m subplots, abundance (number of main stems) of the hypothesized foundation species and other focal species was highest in the three temperate plots (Table 1, Figure 6). The average and total basal area of these same species were highest at Wind River (Tsuga heterophylla) and Harvard Forest (Tsuga canadensis) in the temperate region, at least 2-4-fold greater than the most abundant important species in the tropical plots (Dacryodes excelsa at Luquillo and Alseis blackiana at BCI) (Table 1, Figure 6). Note that focal species were absent in some of 20 × 20-m subplots at both BCI and AM (white cells in the bottom two rows of Figure 6). In contrast, abundance, richness, Shannon diversity and Simpson's diversity of associated species were substantially higher in the three tropical plots than in any of the temperate plots ( Figure 7). Average subplot-wise Bray-Curtis dissimilarity of associated species was higher in the three temperate plots and at Amacayacu than at Luquillo or BCI (Figure 7). Note that only focal species were present in many of the 20 × 20-m subplots at WR and in a handful of 20 × 20-m subplots at HF and TY (white cells in Figure 7).  Figure 6). Note that focal species were absent in some of 20 × 20-m subplots at both BCI and AM (white cells in the bottom two rows of Figure 6). In contrast, abundance, richness, Shannon diversity and Simpson's diversity of associated species were substantially higher in the three tropical plots than in any of the temperate plots ( Figure 7). Average subplot-wise Bray-Curtis dissimilarity of associated species was higher in the three temperate plots and at Amacayacu than at Luquillo or BCI (Figure 7). Note that only focal species were present in many of the 20 × 20-m subplots at WR and in a handful of 20 × 20-m subplots at HF and TY (white cells in Figure 7). Mean DBH Figure 6. Abundance, total basal area, and mean basal area of all hypothesized foundation species and other focal species (Table 1) in the six 20-m rasterized forest dynamics plots. Sites are ordered top-to-bottom by decreasing latitude, and each plot is scaled relative to the 50-ha BCI plot. White squares in the BCI and AM panels indicate none of the focal species were present in that subplot. Figure 6. Abundance, total basal area, and mean basal area of all hypothesized foundation species and other focal species (Table 1) in the six 20-m rasterized forest dynamics plots. Sites are ordered top-to-bottom by decreasing latitude, and each plot is scaled relative to the 50-ha BCI plot. White squares in the BCI and AM panels indicate none of the focal species were present in that subplot.

Codispersion Between Focal Species and Associated-species Diversity at Harvard Forest
As with our initial screen for candidate foundation species, in which we used a graphical "fingerprint" observed for Tsuga canadensis at HF to screen for important species at other sites (Figures 2-5), we first examined the codispersion statistics and graphs illustrating spatial associations between metrics of woody plant diversity and total basal area of Tsuga canadensis relative to the other three dominant tree species at HF within 20 × 20-m subplots.
At all spatial lags (computed from 20-125 m) and directions, the strongest and most consistently significant measures of codispersion between total basal area and metrics of associated-species diversity at Harvard Forest were observed for the foundation species, T. canadensis ( Figure 8, Table A1). Total basal area of T. canadensis was negatively spatially associated at all spatial lags and directions with the total abundance, species richness, Shannon diversity, and inverse Simpson's diversity of associated-woody species, but was positively spatially associated at all spatial lags and directions with average Bray-Curtis dissimilarity ( Figure 8, Table A1). All of these codispersion coefficients-at virtually all spatial lags and directions-were significant relative to those generated with the more conservative toroidal-shift null models ( Figure 9). The few codispersion coefficients that were not significant when tested with the more conservative toroidal-shift null model were significant when tested with the CSR null model ( Figure 10). In contrast, codispersion between total basal area of the other three focal species at HF and all metrics of associated-species diversity were weaker and inconsistently significant (Figures 8-10, Table A1). Only the consistently positive codispersion between Q. rubra and either abundance of associated woody species or Bray-Curtis dissimilarity were always significant relative to null expectation (Figures 9 and 10). Overall, the absolute value of the mean or median codispersion between total basal area of T. canadensis and metrics of associated-species diversity ranged from 2-10-times greater than the mean or median codispersion between the total basal area of any of the other three dominant tree species and measures of associated-species diversity (Table A1).  (Table 1) in the six 20-m rasterized forest dynamics plots. Sites are ordered top-to-bottom by decreasing latitude, and each plot is scaled relative to the 50-ha BCI plot. White squares in the WR, HF, and TY panels indicate that only the focal species were present in that subplot.

Codispersion Between Focal Species and Associated-species Diversity at Harvard Forest
As with our initial screen for candidate foundation species, in which we used a graphical "fingerprint" observed for Tsuga canadensis at HF to screen for important species at other sites (Figures 2-5), we first examined the codispersion statistics and graphs illustrating spatial associations between metrics of woody plant diversity and total basal area of Tsuga canadensis relative to the other three dominant tree species at HF within 20 × 20-m subplots.
At all spatial lags (computed from 20-125 m) and directions, the strongest and most consistently significant measures of codispersion between total basal area and metrics of associated-species diversity at Harvard Forest were observed for the foundation species, T. canadensis ( Figure 8, Table A1). Total basal area of T. canadensis was negatively spatially associated at all spatial lags and directions with the total abundance, species richness, Shannon diversity, and inverse Simpson's diversity of associated-woody species, but was positively spatially associated at all spatial lags and directions with average Bray-Curtis dissimilarity ( Figure 8, Table A1). All of these codispersion coefficients-at virtually all spatial lags and directions-were significant relative to those generated with the more conservative toroidal-shift null models ( Figure 9). The few codispersion coefficients that were not significant when tested with the more conservative toroidal-shift null model were significant when tested with the CSR null model ( Figure 10). In contrast, codispersion between total basal area of the other three focal species at HF and all metrics of associated-species diversity were weaker and inconsistently significant (Figures 8-10, Table A1). Only the consistently positive codispersion between Q. rubra and either abundance of associated woody species or Bray-Curtis dissimilarity were always significant relative to null expectation (Figures 9 and 10). Overall, the absolute value of the mean or median codispersion between total basal area of T. canadensis and metrics of associated-species diversity ranged from 2-10-times greater than the mean or median codispersion between the total basal area of any of the other three dominant tree species and measures of associated-species diversity (Table A1).  Figure 10. Summaries of codispersion values are given in Table A1).   Figure 10. Summaries of codispersion values are given in Table A1).  Figure 10. Summaries of codispersion values are given in Table A1).

Trees at temperate sites
Observed patterns of codispersion between total basal area of focal species and metrics of associated-species diversity at WR and TY were qualitatively similar to those observed for the non-foundation species at HF. At WR, codispersion between total basal area of Tsuga heterophylla and all measures of diversity was consistently negative ( Figure 11) but only significant at all lags and directions for abundance and richness of associated species ( Figure A1). Values of codispersion at WR between T. heterophylla and associated-species diversity were much lower than those observed for T. canadensis but similar to those observed for Q. rubra at HF (Table A1). Codispersion between total basal area of Pseudotsuga menziesii and Bray-Curtis similarity was positive, and was significant for about two-thirds of the lags and directions when assessed with the toroidal-shift null model ( Figure 12) but not the CSR null model ( Figure A1).
At TY, codispersion between associated-species abundance and richness and basal area of two of the oaks-Quercus alba and Quercus rubra-was negative and consistently significant (Figures 13  and 14 and Figure A2). Codispersion between Shannon diversity of associated-species and basal area of Q. rubra was also consistently and significantly negative (Figures 13 and 14 and Figure A2). Codispersion between Bray-Curtis dissimilarity and basal area of Q. alba and Carya ovata were significantly positive and negative, respectively (Figures 13 and 14 and Figure A2). As at WR, the range of observed, significant, codispersion values measured at TY paralleled those measured for Q. rubra at HF (Table A1). A notable exception were the values of codispersion between basal area of Q. rubra at TY and associated-species richness, which were nearly identical to that observed for T. canadensis at HF (Table A1).

Trees at temperate sites
Observed patterns of codispersion between total basal area of focal species and metrics of associated-species diversity at WR and TY were qualitatively similar to those observed for the non-foundation species at HF. At WR, codispersion between total basal area of Tsuga heterophylla and all measures of diversity was consistently negative ( Figure 11) but only significant at all lags and directions for abundance and richness of associated species ( Figure A1). Values of codispersion at WR between T. heterophylla and associated-species diversity were much lower than those observed for T. canadensis but similar to those observed for Q. rubra at HF (Table A1). Codispersion between total basal area of Pseudotsuga menziesii and Bray-Curtis similarity was positive, and was significant for about two-thirds of the lags and directions when assessed with the toroidal-shift null model ( Figure 12) but not the CSR null model ( Figure A1).
At TY, codispersion between associated-species abundance and richness and basal area of two of the oaks-Quercus alba and Quercus rubra-was negative and consistently significant (Figures 13  and 14 and Figure A2). Codispersion between Shannon diversity of associated-species and basal area of Q. rubra was also consistently and significantly negative (Figures 13 and 14 and Figure A2). Codispersion between Bray-Curtis dissimilarity and basal area of Q. alba and Carya ovata were significantly positive and negative, respectively (Figures 13 and 14 and Figure A2). As at WR, the range of observed, significant, codispersion values measured at TY paralleled those measured for Q. rubra at HF (Table A1). A notable exception were the values of codispersion between basal area of Q. rubra at TY and associated-species richness, which were nearly identical to that observed for T. canadensis at HF (Table A1).  Table A1).  Figure 11). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A1 for significance testing using a CSR null model.  Table A1).  Table A1).  Figure 11). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A1 for significance testing using a CSR null model.  Figure 11). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A1 for significance testing using a CSR null model.   Figure 14 and Figure A2. Summaries of codispersion values are given in Table A1).  Figure 13). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A2 for significance testing using a CSR null model.  Table A1).  Figure 14 and Figure A2. Summaries of codispersion values are given in Table A1).  (Figure 13). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A2 for significance testing using a CSR null model.  (Figure 13). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A2 for significance testing using a CSR null model.

Trees at Tropical Sites
Codispersion between tree size or abundance and measures of associated-species diversity at the three tropical sites were much weaker than at the temperate sites (Figures 15-17), and were less frequently statistically significant for most spatial lags or directions (Figures 18-20 and Figures A3-A5). Notable exceptions included: positive codispersion between basal area of Dacryodes excelsa and Bray-Curtis dissimilarity at LFDP ( Figure 15); negative codispersion between basal area of Oenocarpus mapora and associated-species richness and abundance at BCI ( Figure 16); and positive codispersion between basal area of O. mapora and associated-species Simpson's diversity at BCI (Figure 16). In all cases, however, values of codispersion were much less than that observed for Q. rubra at HF (except for codispersion between basal area of O. mapora and associated-species abundance, which was similar to that of Q. rubra at HF) (Table A2). At Amacayacu, the most species-rich site, total basal area of the dominant tree species had small and rarely significant spatial associations with metrics of diversity of associated woody species (Figures 17 and 20 and Figure A5). the three tropical sites were much weaker than at the temperate sites (Figures 15-17), and were less frequently statistically significant for most spatial lags or directions (Figures 18-20 and Figures A3-A5). Notable exceptions included: positive codispersion between basal area of Dacryodes excelsa and Bray-Curtis dissimilarity at LFDP ( Figure 15); negative codispersion between basal area of Oenocarpus mapora and associated-species richness and abundance at BCI ( Figure 16); and positive codispersion between basal area of O. mapora and associated-species Simpson's diversity at BCI (Figure 16). In all cases, however, values of codispersion were much less than that observed for Q. rubra at HF (except for codispersion between basal area of O. mapora and associated-species abundance, which was similar to that of Q. rubra at HF) (Table A2). At Amacayacu, the most species-rich site, total basal area of the dominant tree species had small and rarely significant spatial associations with metrics of diversity of associated woody species (Figures 17 and 20 and Figure A5).

Codispersion Between Focal Understory Species and Associated-species Diversity
Understory species were included in our focal species list at two temperate sites (WR, TY) and two tropical sites (LFDP, AM) ( Table 1; Figures A6-A11; Tables A3, A4). Of all of these, consistently significant codispersion between basal area and any metrics of associated-species diversity was observed only for Acer circinatum at WR, Cornus florida at TY, and Rinorea lindeniana at AM (Figures A10  and A11). At WR, codispersion between basal area of A. circinatum and both Shannon and Simpson diversity was consistently negative, but was consistently positive with Bray-Curtis dissimilarity ( Figures A6-A8). At AM, the pattern was reversed. There, codispersion between basal area of R. lindeniana and both Shannon and Simpson's diversity were consistently positive, but was consistently negative with Bray-Curtis dissimilarity ( Figures A9-A11), and were of the same magnitude or lower than that observed for Q. rubra at HF (Table A4). At TY, codispersion between basal area of C. florida and both associated-species richness and Bray-Curtis similarity were consistently negative ( Figure A6 and of low magnitude (Table A3).

Codispersion Between Focal Understory Species and Associated-species Diversity
Understory species were included in our focal species list at two temperate sites (WR, TY) and two tropical sites (LFDP, AM) ( Table 1; Figures A6-A11; Tables A3 and A4). Of all of these, consistently significant codispersion between basal area and any metrics of associated-species diversity was observed only for Acer circinatum at WR, Cornus florida at TY, and Rinorea lindeniana at AM (Figures A10 and A11). At WR, codispersion between basal area of A. circinatum and both Shannon and Simpson diversity was consistently negative, but was consistently positive with Bray-Curtis dissimilarity (Figures A6-A8). At AM, the pattern was reversed. There, codispersion between basal area of R. lindeniana and both Shannon and Simpson's diversity were consistently positive, but was consistently negative with Bray-Curtis dissimilarity ( Figures A9-A11), and were of the same magnitude or lower than that observed for Q. rubra at HF (Table A4). At TY, codispersion between basal area of C. florida and both associated-species richness and Bray-Curtis similarity were consistently negative ( Figure A6 and of low magnitude (Table A3).  Figure 19 and Figure A4.    Figure 19 and Figure A4.    Figure 15). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A3 for significance testing using a CSR null model.  Figure 19. Statistical significance (red: P ≤ 0.05; blue: P > 0.05) of the codispersion coefficients calculated between measures of diversity and total basal area of the two focal tree species in 20 × 20-m subplots in the 50-ha BCI plot ( Figure 16). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A4 for significance testing using a CSR null model.  Figure 15). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A3 for significance testing using a CSR null model.  Figure 15). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A3 for significance testing using a CSR null model.  Figure 19. Statistical significance (red: P ≤ 0.05; blue: P > 0.05) of the codispersion coefficients calculated between measures of diversity and total basal area of the two focal tree species in 20 × 20-m subplots in the 50-ha BCI plot ( Figure 16). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A4 for significance testing using a CSR null model. Figure 19. Statistical significance (red: P ≤ 0.05; blue: P > 0.05) of the codispersion coefficients calculated between measures of diversity and total basal area of the two focal tree species in 20 × 20-m subplots in the 50-ha BCI plot ( Figure 16). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A4 for significance testing using a CSR null model.   Figure 17). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A5 for significance testing using a CSR null model.

Discussion
Foundation species in forests define and structure entire ecological communities primarily through non-trophic effects such as providing physical points of attachment for associated species (e.g., epiphytes), creating habitat for associated fauna (e.g., tree-holes), modulating hydrological flow of adjacent streams, or altering the chemical composition of soils [2,3,52]. Forest foundation species most frequently are common, abundant, and large, and have been identified more frequently in temperate forests [2,53] than in tropical ones [2]. Ellison et al. [2] suggested that foundation tree species would be more likely in comparatively species-poor temperate forests because functional redundancy among trees there would be much less common than in species-rich tropical forests (e.g., [54,55]). In partial support of this idea, Lamanna et al. presented data showing a peak in the size of functional-trait space for tree assemblages growing at mid-latitudes (where our temperate plots are located) but more narrow functional-trait space for tropical assemblages [56]; see also [57]. These analyses imply that foundation tree species should be less common in very high latitude ecosystems (examined in more detail by [56] than by [57]). In those ecosystems, low-growing perennial, cushion-and tussock-forming plants have been found to be foundation species (e.g., [58]).
Indeed, the goal of testing community ecology's neutral theory (which posits that differences among trophically similar species, such as trees, are irrelevant to their demographic success and subsequent spatial distribution), deriving from Hubbell's work in a 13-ha plot in Costa Rica [59,60], in part motivated the establishment of the 50-ha plot at BCI in 1980 and subsequently other tropical forest dynamics plots. Although neutral theory has been supported for some distributional patterns in many tropical forests, it fails in others [61]. Niche theory (which posits functional non-equivalence of co-occurring, trophically similar species) more strongly fits data observed in temperate forests [62] and many tropical ones, too [61]. Thus, there is no a priori reason to suspect on functional grounds that foundation species should be less likely to occur in tropical forests than in temperate ones. We note also that the definition of foundation species [1,2] implicitly assumes that they always have non-neutral effects on the systems they define. However, it is plausible that the high abundance of these species causes them to have neutral influences on co-occurring species via ecological drift (i.e., by keeping local  Figure 17). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a toroidal-shift null model. See Figure A5 for significance testing using a CSR null model.

Discussion
Foundation species in forests define and structure entire ecological communities primarily through non-trophic effects such as providing physical points of attachment for associated species (e.g., epiphytes), creating habitat for associated fauna (e.g., tree-holes), modulating hydrological flow of adjacent streams, or altering the chemical composition of soils [2,3,52]. Forest foundation species most frequently are common, abundant, and large, and have been identified more frequently in temperate forests [2,53] than in tropical ones [2]. Ellison et al. [2] suggested that foundation tree species would be more likely in comparatively species-poor temperate forests because functional redundancy among trees there would be much less common than in species-rich tropical forests (e.g., [54,55]). In partial support of this idea, Lamanna et al. presented data showing a peak in the size of functional-trait space for tree assemblages growing at mid-latitudes (where our temperate plots are located) but more narrow functional-trait space for tropical assemblages [56]; see also [57]. These analyses imply that foundation tree species should be less common in very high latitude ecosystems (examined in more detail by [56] than by [57]). In those ecosystems, low-growing perennial, cushion-and tussock-forming plants have been found to be foundation species (e.g., [58]).
Indeed, the goal of testing community ecology's neutral theory (which posits that differences among trophically similar species, such as trees, are irrelevant to their demographic success and subsequent spatial distribution), deriving from Hubbell's work in a 13-ha plot in Costa Rica [59,60], in part motivated the establishment of the 50-ha plot at BCI in 1980 and subsequently other tropical forest dynamics plots. Although neutral theory has been supported for some distributional patterns in many tropical forests, it fails in others [61]. Niche theory (which posits functional non-equivalence of co-occurring, trophically similar species) more strongly fits data observed in temperate forests [62] and many tropical ones, too [61]. Thus, there is no a priori reason to suspect on functional grounds that foundation species should be less likely to occur in tropical forests than in temperate ones. We note also that the definition of foundation species [1,2] implicitly assumes that they always have non-neutral effects on the systems they define. However, it is plausible that the high abundance of these species causes them to have neutral influences on co-occurring species via ecological drift (i.e., by keeping local populations of associated species at low abundance, contributing to their demographic stochasticity). This might help explain why the basal area of the potential foundation species we examined was often negatively correlated with local species richness but positively correlated with beta diversity of associated species.
Despite the ecological importance of foundation species, it has proven challenging to identify plausible foundation species from among groups of locally dominant species or those that are perceived by various criteria to have importance in particular systems. It has taken decades of observational and experimental work to support the hypothesis that Tsuga canadensis is a foundation species in eastern North American forests [16] and that Pinus albicaulis is a foundation species in high-elevation western North American forests [53]. In both cases, their foundational characteristics were identified only as these species were declining [2]. If foundation species really are of critical importance for associated biodiversity, it is crucial that they be identified and protected while they are still common, locally abundant, and present in their full size distribution [6,27]. Because experiments on forest trees can take decades to centuries (and are generally impractical with longer-lived species), alternative approaches to identifying foundation forest species are needed.
Our statistical screening of abundance and size data from six large forest dynamics plots suggested particular characteristics of potential foundation species from two different size-frequency and diameter-abundance distributions (Figures 2-5). The foundation species T. canadensis was a notable outlier from the "reverse-J" distribution of abundance versus mean diameter in the entire Harvard Forest 35-ha forest dynamics plot (Figure 2) and had a much greater range and spread of values of total basal area versus abundance in the 20-m rasterized plot (Figure 3). A similarly small set of species that had similar characteristics were identified in the other forest dynamics plots we analyzed (Figures 4  and 5), with more candidates identified in the temperate plots than the tropical ones. This result is not especially surprising, as the more even abundance distributions of most tropical forests suggest that any species with both of these screening characteristics would be less common in the tropics. Alternatively, groups of species with similar traits (e.g., the Eschweilera species [Lecythidaceae] at AM, Dipterocarpaceae in Southeast Asia, or Quercus species at TY) or particular characteristics of individuals themselves (e.g., very large trees [26,27]) could have disproportionate influence akin to foundational characteristics in the forests in which they occur. This hypothesis is in need of additional exploration.
We also note that departures from the "reverse-J" size-frequency distribution may depend on when in succession a stand is sampled [63]. For example, abundance of early seral species (such as Pseudotsuga menziesii at WR) may initially be very abundant but relatively small. Because P. menziesii has been decreasing in abundance for at least the last 200 years at WR [64], it might be identified as a candidate foundation species in a 300-year-old stand, but not in the older forest dynamics plot at WR. That both the subcanopy A circinatum and the late-seral T. heterophylla basically mirror the position of P. menziesii in the size-frequency plot (Figure 4) adds some support for this hypothesis. Conversely, the foundational characteristics of T. heterophylla suggested by Figures 4 and 5 might not yet be manifest this early in succession (cf. [63]). Related to this, some of the species that fall off the "reverse-J" have intermediate shade-tolerance and their role is likely to change during succession. Subsequent analysis of spatial relationships between these candidate foundation species, along with others of known ecological importance (Table 1), and metrics of diversity of associated species pointed to attributes of T. canadensis shared by only a few other species. These included codispersion values >|0.25| at virtually all spatial lags and directions for two or more metrics of diversity that were often negative for measures of abundance or alpha-diversity (within-subplot species richness, Shannon diversity, or inverse-Simpson's diversity) but positive for a measure of beta diversity (average subplot Bray-Curtis dissimilarity). That is, foundation species are negatively correlated with diversity of associated species at local (within-subplot) spatial scales where they dominate the vegetation but are positively associated with beta diversity at the larger (between subplot) spatial scale of the whole-plot. In virtually all cases (a notable exception being Pseudotsuga menziesii at WR), little effect of underlying environmental patchiness was observed: the toroidal-shift null model and the CSR null model gave nearly identical results, with the toroidal-shift null model being predictably more conservative.
By these criteria, candidates for additional, more detailed study for their foundational roles include Acer circinatum at WR, Dacryodes excelsa at Luquillo, and Oenocarpus mapora at BCI (Figures 15  and 16 and Figure A6). In the western Cascades of Oregon, A. circinatum is a major component of the subcanopy, where it can grow rapidly, have high biomass, and form broad canopies that suppress other species [65,66]. This suppression of recruiting species could lead to negative association between A. circinatum basal area and local diversity of associated species (as is also seen with T. canadensis in eastern North American forests). At Luquillo, adult D. excelsa dominates ridges where trees form root grafts that may, in part, explain their relatively low hurricane mortality (<1%) [67,68]. Dacryodes excelsa saplings also grow more rapidly than saplings of other late successional species [68], which may explain the lack of pole-sized individuals in their populations [67]. Finally, O. mapora alters litter quality and quantity and the local microenvironment, and limits establishment of small-seeded and shade-intolerant species at BCI [69].
Disturbance frequency also may play a role in both statistical "fingerprints" of foundation species. For example, intermediate disturbance regimes (moderately frequent disturbances of moderate size and severity, such as small understory fires, multi-tree blowdowns, or selective harvesting) could push mid-successional, shade-intermediate species such as Q. rubra or A. rubrum at Harvard Forest back into the "reverse-J" distribution. At the same time, successional dynamics of these species (especially recruitment and prolonged growth) could result in transient conditions whereby they have both relatively high basal area and large numbers of small individuals in any given subplot. Where these species fall within the fingerprints, therefore, may indicate how their life-history traits are aligned with the frequency and type of local disturbances. How foundation species exert control on disturbance regimes in forests has not been previously explored. It would be useful to explore how positive feedbacks between a forest foundation species and local disturbances (especially fire, small windstorms or microbursts, and treefalls) could affect size-frequency and diameter-abundance distributions of associated species. These processes also are likely to differ in temperate and tropical forests (e.g., [70][71][72]).
Although local diversity of associated woody species is lower in the presence of foundation species, local diversity of other taxa associated with foundation species is often higher or composed of unique assemblages. Non-woody vascular and especially nonvascular plant taxa, and fauna are rarely sampled at the individual tree scale or in every 20 × 20-m subplot in large forest dynamics plots, but a handful of other studies have examined associations between these taxa and either T. canadensis or our candidate foundation species. For example, T. canadensis hosts unique assemblages of arthropods [73][74][75][76]. Acer circinatum is associated with a diverse and abundant assemblage of epiphytes [77]. At Luquillo, higher numbers of invertebrate species are associated with Dacryodes excelsa than with other dominant canopy species (including Cecropia schreberiana) [78], but abundances and herbivory rates of invertebrates on D. excelsa are lower than on other co-occurring species [79]. Several understory palms, including O. mapora, recently have been described as "terrestrial litter trappers" [80] that collect leaf litter at the base of their foliage. Two other litter-trapping palms, Asterogyne sociale and A. spectabilis, have very large numbers of arthropods in the trapped litter [81]; comparable studies with O. mapora would provide a useful comparison with the known and potential foundation species discussed here.
Finally, species like T. canadensis, P. menziesii and D. excelsa "name" their communities: hemlock, Douglas-fir, and tabonuco forests, respectively. The forest stands they dominate are easily identifiable as distinctive patches on the landscape. When a foundation species like T. canadensis declines, the associated regional or landscape-scale beta diversity is likely to decline along with it. Such regional homogenization-as indicated by lower Bray-Curtis values-is more characteristic of the majority of our focal species and would be a likely consequence of the loss of our candidate foundation species. It is precisely this larger-scale heterogeneity that may be the greatest value gained by protecting foundation species while they are still common.

Conclusions
Foundation species define and structure ecological communities but are difficult to identify before they are declining. Comparative analyses of six large forest dynamics plots spanning a temperate-to-tropical gradient in the Western Hemisphere identified statistical "fingerprints" of potential foundation species. Specifically, known and potential foundation species are outliers from the common "reverse-J" size-frequency distribution, and have negative effects on alpha diversity and positive effects on beta diversity at most spatial lags and directions. Potential foundation species appear to be more likely in temperate forests, whereas in tropical forests, foundational characteristics may be more likely to occur in species groups within genera or families. As foundation species (or species groups) decline, local (alpha) diversity may increase but landscape-scale (beta) diversity is likely to decline along with them. Preservation of beta biodiversity may be the most important consequence of protecting foundation species while they are still common. Acknowledgments: An early version of this paper was presented at the 2014 International LTER meeting in Valdivia, Chile. We especially thank Ronny Vallejos for development of software to measure and plot codispersion, and for his continuing collegial collaboration with A.M.E., H.L.B., and B.S.C. Two anonymous reviewers provided helpful comments that substantially improved the final version of the paper. We acknowledge logistical support from the Gifford Pinchot National Forest, the USDA Forest Service Pacific Northwest Research Station, and the volunteers individually listed at http://wfdp.org for the WR plot. We also thank the many field technicians who helped census the HF plot; Jason Aylward for field supervision, data screening and database management at HF; and John Wisnewski and the Harvard Forest Woods Crew for providing materials, supplies, and invaluable field assistance with plot logistics. We thank the Tyson Research Center staff and more than 100 students and researchers that have contributed to tree censuses at TY. We sincerely thank the many volunteers who collected tree census data at LFDP. Vegetation data from BCI are part of the BCI forest dynamics research project founded by S.P. Hubbell and R.B. Foster and now managed by R. Condit, S. Lao, and R. Pérez through the Center for Tropical Forest Science (CTFS) and the Smithsonian Tropical Research Institute (STRI) in Panamá. We are very grateful for the assistance of coworkers from the Comunidad de Palmeras and the students of forest engineering from the Universidad Nacional de Colombia who collected tree census data at the AM plot. This paper is a contribution of the Harvard Forest and Luquillo LTER programs.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Table A1. Codispersion statistics for focal tree species in temperate forests (Figures 8-20 and Figures A1 and A2).              Figure A6). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a CSR null model.   Figure A6). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a CSR null model. Figure A8. Statistical significance (red: P ≤ 0.05; blue: P > 0.05) of the codispersion coefficients calculated between measures of diversity and total basal area of the focal understory species in 20 × 20-m subplots in the 25.6-ha Wind River (Acer circinatum, Taxus brevifolia) and 20-ha Tyson Research Center (Asimina triloba, Cornus florida, Frangula caroliniana, Lindera benzoin) plots ( Figure A6). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a CSR null model.      Figure A9). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a CSR null model.   Figure A9). Statistical significance was determined by comparing observed codispersion at each spatial lag with a distribution of 199 spatial randomizations of a CSR null model.