A key issue in evolutionary biology is whether selection acting at levels higher than the individual can cause evolutionary change. If it can, then conceptual and empirical studies must consider how selection operates at multiple levels of biological organization. Here, we test the hypothesis that estimates of broad-sense community heritability, H²_C, can be used to predict the evolutionary response by community-level phenotypes when community-level selection is imposed. Using an approach informed by classic quantitative genetics, we made three predictions. First, when we imposed community-level selection, we expected a significant change in the average phenotype of arthropod communities associated with individual tree genotypes [we imposed selection by favoring high and low NMDS (nonmetric multidimensional scaling) scores that reflected differences in arthropod species richness, abundance and composition]. Second, we expected H²_C to predict the magnitude of the community-level response. Third, we expected no significant change in average NMDS scores with community-level selection imposed at random. We tested these hypotheses using three years of common garden data for 102 species comprising the arthropod communities, associated with nine clonally replicated Populus angustifolia genotypes. Each of our predictions were met. We conclude that estimates of H²_C account for the resemblance among communities sharing common ancestry, the persistence of community composition over time, and the outcome of selection when it occurs at the community level. Our results provide a means for exploring how this process leads to large-scale community evolutionary change, and they identify the circumstances in which selection may routinely act at the community level.

ESTIMATES OF BROAD-SENSE COMMUNITY HERITABILITY (H²_C)

Keith et al., (2010) designed an experiment to examine the effect of genotypic variation in plants on arthropod community organization using an 18-yr-old common garden with replicated clones of nine different Populus angustifolia genotypes. All trees planted within the common garden were collected from a single interbreeding population (Martinsen et al., 2001). Trees were identified using molecular markers that allowed exclusion of hybrids and inclusion of genetic variants characteristic of “pure” P. angustifolia. Genotypes represented in the common garden had been haphazardly selected from trees growing along the Weber River in northern Utah, USA and were planted in a haphazard design. Nine tree genotypes with four to seven replicates each were selected from existing stocks, yielding a total of 44 trees, whose average height was 10-15m.

Following Wimp et al., (2005), Keith et al., (2010) censused arthropod communities on trees in each of three years (2004-06) and summarized the 102 species (67 families in 12 orders) using nonmetric multidimensional scaling ordination (NMDS). This approach evaluated arthropod community composition as a quantitative trait (Bradshaw & Stettler, 1995), as in studies of diverse multivariate plant traits including phenology, phytochemistry, morphology, sink-source relationships (e.g., Holeski et al., 2012) and interactions with each other that result in different interaction networks for different tree genotypes (e.g., Lau et al., 2016; Keith et al., 2017). Censuses of arthropod community composition for each tree provided individual trait values for this quantitative character according to standard quantitative genetics methods (e.g., in this case, among lineages of clones grown in a common environment; Falconer & McKay, 1996; Shuster et al., 2006).

COMMUNITY-LEVEL SELECTION BASED ON NMDS SCORES

HOW NMDS SCORES DESCRIBE COMMUNITY PHENOTYPE

We assumed that differences in arthropod abundance arose from genetic interactions that had occurred between arthropod genotypes and the genotypes of their tree hosts that either favored their persistence or disfavored the persistence of other individuals or species not found among the sampled arthropod community (c.f., Shuster et al., 2006). This assumption of genetics-based interactions has been experimentally confirmed for several herbivorous species used in our analyses. For example, transfer experiments showed pronounced differences in resistance among tree genotypes, including arthropod preferences for trees where their survival was greatest and avoidance of tree genotypes where their survival was lowest (e.g., aphids, Whitham et al. 1989; mites, Evans et al. 2008). Experiments have also shown intraspecific genetic differences in mites and aphids in which some genotypes do best on some tree genotypes, but not others (e.g., Evans et al., 2008; Smith et al., 2020). Because some insects, such as aphids can affect 100s of other species including insects, spiders, fungi and birds, their genetics-based interactions with individual tree genotypes can directly and indirectly affect whole communities of organisms (e.g., Dickson & Whitham, 1996; Smith et al., 2011; Keith et al., 2017) including ecosystem processes such as nutrient cycling (Schweitzer et al., 2005).

Thus, the above and numerous other studies (reviews Whitham et al., 2012, 2020), confirm the use of NMDS scores to quantify a multivariate phenotype arising from the genetic interactions of arthropod symbiont and tree genotypes. For this reason, we expected that directional selection favoring communities expressing particular NMDS scores would change the distribution of NMDS scores in the next year, thereby simulating an evolutionary response to community-level selection.

CONTROLLING FOR TEMPORAL AUTOCORRELATION AMONG YEARS

To control for potential autocorrelation of NMDS scores among study years, possibly due to shared maternal environments, position within the common garden, and persistent induced or epigenetic responses among years, we performed a 2-way ANOVA of the NMDS data with tree genotype, year, and their interaction as factors, as well as a Durbin-Watson test for autocorrelation on these data.

DIRECTIONAL COMMUNITY-LEVEL SELECTION ON NMDS SCORES

To simulate community-level selection based on variation in community phenotype, we selected the NMDS scores with the 10 largest and 10 smallest values from the array of 44 one-dimensional NMDS scores generated by the 2004 communities within the common gardens described in Keith et al., (2010). This procedure identified two groups of trees whose abundances of arthropod species led to the highest and lowest NMDS ordinations for the 2004 sample. Analogous to the statistical analysis of quantitative traits, this procedure provided no ecological information on the reasons for these community phenotypic similarities or differences, only that they occurred. We consider such ecological anonymity one of the strengths of this approach.

We next identified the 20 trees in the common garden from which these 2004 scores were drawn, isolated these trees and the abundances of the species within these trees’ insect communities from the rest of the 2005 sample, and then performed a one-dimensional NMDS ordination on these 2005 communities. From this 2005 ordination, we then selected the NMDS scores with the 5 largest and 5 smallest values. We identified the 10 trees in the 2006 common garden from which these 2005 scores were drawn, we isolated these trees and the abundances of the species within these trees’ insect communities from the rest of the 2006 sample, and we again performed a one-dimensional NMDS ordination on these 2006 communities. We tabulated the genotypes and NMDS scores for each episode of selection in Appendix B.

RANDOM COMMUNITY-LEVEL SELECTION ON NMDS SCORES

Consistent with studies of selection acting on quantitative traits, we expected a random selection of NMDS community scores in each year to produce no evolutionary response, because a random selection of NMDS scores summarizing community phenotypes in each episode of selection would produce no distinguishable change in the average community phenotype. Thus, as a control for our community-level selection experiment on NMDS scores described above, we performed the same procedures on the 2004-06 samples of trees and their communities, except with the NMDS scores chosen at random without replacement using a random number generator. We replicated our control procedure 10 times to simulate five independent selection series each on high and low NMDS lineages. We plotted both sets of results.

MEASURING THE RESPONSE TO COMMUNITY-LEVEL SELECTION

We calculated the mean and 95% confidence limits for the NMDS scores in each of the following samples: (1) the original 2004 sample of scores [2004 initial; N=44 communities], (2) the sample of the 10 high and 10 low scores in the 2004 sample [2004 selected, N=20 communities], (3) the 20 scores of the 2005 communities on the trees identified by the 2004 selected scores [2005 response; N=20 communities], (4) the sample of 5 high and 5 low scores in the 2005 sample [2005 selected; N=10 communities], and (5) the 2006 communities on trees identified by the 2005 selected scores [2006 response; N=10 communities]. We also calculated the mean and 95% confidence limits for each of the five sets of randomly selected NMDS scores used as controls for our high and low NMDS lineages. These methods allowed us to simulate the effects of directional selection on the NMDS-quantified phenotypes of the arthropod communities inhabiting trees in the common garden, and to compare the results of that selection with randomly selected NMDS scores from the same communities in samples of similar size.

For each i-th episode of community-level selection, where i=2004 or 2005 (no data were available to document a response to selection after 2006 so this episode was not included), we estimated the community phenotypic mean, Z_Ci, and standard deviation, s_ZCi, of the initial distribution of NMDS scores, as well as the mean, Z_Ci*, of the distribution of selected NMDS scores from the trees in that selection episode, for our experimental and control communities. We calculated the community selection differential, S_Ci, for each i-th selection episode as the difference between the average NMDS phenotype of the selected and initial samples of NMDS scores, standardized by the standard deviation of the initial sample of scores, or,

            S_Ci = (Z_Ci* – Z_Ci) / s_ZCi.                                                            (1)

For comparison, we used the tabulated values for selection differentials in Becker, (1985, p. 161-174) based on the number of individuals in each selection episode. We estimated the cumulative selection differential over our two episodes of community-level selection as the sum of the two selection differentials estimated for the 2004 and 2005 samples. Note that the magnitude of the selection differential depends on the size of the population before selection and the number of selected individuals (Becker, 1985). Because our simulation progressively reduced the numbers of communities included within each episode of community selection, the community selection differentials were expected to become proportionally smaller.

We estimated the response to community-level selection in each i-th selection episode, R_Ci, as the difference between the average NMDS scores from trees in the year after selection was imposed, Z_Ci+1, and the average NMDS scores comprising the initial set of communities before selection, Z_Ci, or

                        R_Ci = Z_Ci+1 – Z_Ci.                                                                      (2)

Following the breeders’ equation (R = h²S; Falconer & McKay, 1996), and methods described in Wade et al., (1996) we estimated the realized community heritability in each i-th episode of selection (2004, 2005) as the ratio of the response to community-level selection, R_C(i), to the community selection differential, S_C(i) or

                        R_C(i)/S_C(i) = H²_{C(realized[i])}.                                                          (3)      

We estimated the realized community heritability overall, H²_{C(realized[total])}, as the ratio of the cumulative response to selection, ΣR_C(i), to the cumulative community selection differential, ΣS_C(i), over 2004 and 2005, which we compared with estimates of broad-sense community heritability, H²_C, from Keith et al., (2010).

The data from Keith et al. (2010, attached) consists of abundance data from an 18-yr-old common garden with replicated clones of nine different Populus angustifolia genotypes. All trees planted within the common garden were collected from a single interbreeding population. Trees were identified using molecular markers that allowed exclusion of hybrids and inclusion of genetic variants characteristic of “pure” P. angustifolia. Genotypes represented in the common garden had been haphazardly selected from trees growing along the Weber River in northern Utah, USA and were planted in a haphazard design. Nine tree genotypes with four to seven replicates each were selected from existing stocks, yielding a total of 44 trees, whose average height was 10-15m.

The community heritability calculator is an Excel spreadsheet that can be used to estimate broad sense community heritability from analysis of variance output using one-dimensional NMDS scores on community abundance data.