A framework for evaluating the influence of climate, dispersal limitation, and biotic interactions using fossil pollen associations across the late Quaternary

Environmental conditions, dispersal lags, and interactions among species are major factors structuring communities through time and across space. Ecologists have emphasized the importance of biotic interactions in determining local patterns of species association. In contrast, abiotic limits, dispersal limitation, and historical factors have commonly been invoked to explain community structure patterns at larger spatiotemporal scales, such as the appearance of late Pleistocene no-analog communities or latitudinal gradients of species richness in both modern and fossil assemblages. Quantifying the relative influence of these processes on species co-occurrence patterns is not straightforward. We provide a framework for assessing causes of species associations by combining a null-model analysis of co-occurrence with additional analyses of climatic differences and spatial pattern for pairs of pollen taxa that are significantly associated across geographic space. We tested this framework with data on associations among 106 fossil pollen taxa and paleoclimate simulations from eastern North America across the late Quaternary. The number and proportion of significantly associated taxon pairs increased over time, but only 449 of 56 194 taxon pairs were significantly different from random. Within this significant subset of pollen taxa, biotic interactions were rarely the exclusive cause of associations. Instead, climatic or spatial differences among sites were most frequently associated with significant patterns of taxon association. Most taxon pairs that exhibited co-occurrence patterns indicative of biotic interactions at one time did not exhibit significant associations at other times. Evidence for environmental filtering and dispersal limitation was weakest for aggregated pairs between 16 and 11 kyr BP, suggesting enhanced importance of positive species interactions during this interval. The framework can thus be used to identify species associations that may reflect biotic interactions because these associations are not tied to environmental or spatial differences. Furthermore, temporally repeated analyses of spatial associations can reveal whether such associations persist through time.

Environmental conditions, dispersal lags, and interactions among species are major factors structuring communities through time and across space. Ecologists have emphasized the importance of biotic interactions in determining local patterns of species association. In contrast, abiotic limits, dispersal limitation, and historical factors have commonly been invoked to explain community structure patterns at larger spatiotemporal scales, such as the appearance of late Pleistocene no-analog communities or latitudinal gradients of species richness in both modern and fossil assemblages. Quantifying the relative infl uence of these processes on species co-occurrence patterns is not straightforward. We provide a framework for assessing causes of species associations by combining a null-model analysis of co-occurrence with additional analyses of climatic diff erences and spatial pattern for pairs of pollen taxa that are signifi cantly associated across geographic space.
We tested this framework with data on associations among 106 fossil pollen taxa and paleoclimate simulations from eastern North America across the late Quaternary. Th e number and proportion of signifi cantly associated taxon pairs increased over time, but only 449 of 56 194 taxon pairs were signifi cantly diff erent from random. Within this significant subset of pollen taxa, biotic interactions were rarely the exclusive cause of associations. Instead, climatic or spatial diff erences among sites were most frequently associated with signifi cant patterns of taxon association. Most taxon pairs that exhibited co-occurrence patterns indicative of biotic interactions at one time did not exhibit signifi cant associations at other times. Evidence for environmental fi ltering and dispersal limitation was weakest for aggregated pairs between 16 and 11 kyr BP, suggesting enhanced importance of positive species interactions during this interval. Th e framework can thus be used to identify species associations that may refl ect biotic interactions because these associations are not tied to environmental or spatial diff erences. Furthermore, temporally repeated analyses of spatial associations can reveal whether such associations persist through time. aspects of communities through time and across space (Ackerly 2003, Blois et al. 2013a, Dalsgaard et al. 2013. However, interactions among species also have played a major role in structuring community composition and functioning (Jablonski 2008, Blois et al. 2013c, Wisz et al. 2013. Given recent interest in understanding how climate change may lead to new biotic interactions and unexpected ecological dynamics (Zarnetske et al. 2012, Blois et al. 2013c, there is a critical need to disentangle the joint eff ects of abiotic and biotic factors on community dynamics. Previous work on assemblage structure has quantifi ed community pattern as a single index -such as the number of species or the number of checkerboard pairs -that is then subject to null model analysis. Even a moderatelysized assemblage (e.g. a dataset with multiple sites and multiple species at each site) contains many potential species pairs, however, each of which may exhibit positive, negative, or random associations. In many cases, single metrics that summarize an entire assemblage can be deceptive , and it is more instructive to analyze individual pairs of species (Sfenthourakis et al. 2006). Gotelli and Ulrich (2010) use an empirical Bayes approach (Efron 2005) to control for the potentially large number of false positives that can emerge with the analysis of many species pairs. Th is kind of analysis allows for a determination of the relative frequency of positively, negatively, and randomly associated species pairs. However, non-random species associations are not necessarily caused by species interactions. Th us, a central dilemma is how to distinguish non-random species associations produced by actual species interactions from those produced by environmental fi ltering or dispersal limitations. All three processes can operate singly or in concert to generate both positive and negative species associations.
Ecologists working with modern faunas often explicitly or implicitly limit comparisons to a set of environmentally similar and spatially adjacent sites for which dispersal limitation is unlikely to be important (Phillips et al. 2003, Zhang et al. 2011). In such systems, it is reasonable to attribute non-random species associations to species interactions. Th e eff ects of species interactions certainly can be scaled up to larger spatial and temporal domains (Jablonski 2008, Gilman et al. 2010, Baiser et al. 2012, Blois et al. 2013c), and MacArthur (1972 argued explicitly for this scaling in his fi nal book, Geographical ecology. However, the eff ects of environment, dispersal, and history become progressively more important at larger spatial and temporal scales, and it is diffi cult to untangle them from the eff ects of species interactions (Ricklefs 2004).
Th is dilemma is illustrated clearly in fossil records. Th ese records usually encompass timescales at which environment, dispersal, and biotic interactions are all potentially important controls on species distributions and the associations among species, yet usually their eff ects cannot be directly observed. Often we have information only about species occurrences across space and through time, and perhaps information about past environments. Only rarely can we infer actual biotic interactions in fossil systems (Wilf et al. 2001, Kowalewski 2002, Currano et al. 2010, Pe ñ alver et al. 2012, Blois et al. 2013c, making it diffi cult to confi dently attribute the causes of past species associations to the infl uence of environmental similarity, interactions with other species, or other factors. Here, we provide a framework for inferring the importance of biotic interactions, dispersal limitation, and abiotic eff ects on positive and negative species associations. Th is framework can be applied to species associations measured at any spatial or temporal scale, but we illustrate it in an analysis of eastern North American plant assemblages based on fossil pollen data from the past 21 000 yr. Previous work on both individual species and communities has demonstrated that changes in fossil pollen assemblages across space and time are tightly linked with climate, particularly in the latest Pleistocene and early Holocene (Grimm et al. 1993, Williams et al. 2002, Shuman et al. 2004, Yu 2007, Blois et al. 2013a. Indeed, the tight linkages between vegetation and climate make fossil pollen data an excellent proxy for reconstructing past climates (Viau et al. 2006, Bartlein et al. 2011. Additionally, the recognition of individualistic species responses to deglaciation, the resulting formation of noanalog communities during the Pleistocene -Holocene transition, and the attribution of these communities to no-analog climates (Williams et al. 2004) suggest that species interactions should not be dominant drivers of community patterns, especially during the height of the no-analog period from 17 -11 kyr BP. Other evidence suggests that not all changes in vegetation can be attributed exclusively to climate. For example, at several sites in the Great Lakes region, the loss of megaherbivores and their associated species interactions at the end of the Pleistocene may have contributed to the formation of no-analog plant communities (Gill et al. 2009(Gill et al. , 2012. Additionally, rapid changes in climate during this time period led to dispersal lags in some European tree species Skov 2005, 2007), potentially creating transient species associations. In eastern North America, the infl uence of climate on fossil pollen assemblages was lower during the late Pleistocene than during the Holocene, indicating that biotic interactions or dispersal limitation may have been relatively more important (Blois et al. 2013a). In sum, environmental fi ltering, dispersal limitation, and biotic interactions all likely contributed to structuring fossil pollen taxa and assemblages across the past 21 000 yr, but the eff ects of each process have not been explicitly teased apart.
In this paper, we classify which associations among pollen taxa are likely to be signals of direct biotic interactions versus indirect associations due to dispersal limitation or environmental fi ltering. To do this, we 1) determine with pairwise null model analysis  which taxon associations are signifi cant across space and whether these associations are persistent through time, and 2) develop methods to tease apart the relative infl uences of climate, dispersal limitation, and biotic interactions in explaining those associations. We fi rst illustrate the method for a hypothetical presence -absence data matrix for a single time point. We then use paired fossil pollen and paleoclimate datasets from the last 21 000 yr in eastern North America to test this approach and examine changes in taxon associations through time. Th is framework applies to any paired datasets of species presence -absence and environmental characteristics at any time scale, whether paleontological or not.

Null model analyses of species co-occurrence
A large literature on null model analyses of species cooccurrence has accumulated over the past 80 yr (Harvey et al. 1983, Gotelli andGraves 1996). Th e initial impetus for these analyses was to ask whether co-occurrence patterns in replicated local assemblages were any diff erent than might be expected by chance (Connor and Simberloff 1979). Deviations from this simple null hypothesis were thought to refl ect species interactions (primarily interspecifi c competition), but might also refl ect environmental fi ltering or dispersal limitation. More recent analyses have emphasized the contribution of individual species pairs to assemblage-level patterns of species co-occurrence (Cardillo and Meijaard 2010). Th e typical dataset for a null model analysis is organized into a binary presence -absence matrix (McCoy and Heck 1987) in which rows are taxa, columns are sites or samples, and entries indicate the presence (1) or absence (0) of a particular species in a site. Th ese matrices are commonly collected by community ecologists who sample contemporary assemblages, but replicated fossil data based on standardized sampling methods can be represented and analyzed in exactly the same way.

Quantifying species co-occurrence
Th e strength of association between two species can be quantifi ed by the C-score (Stone and Roberts 1990). For a particular pair of species in the matrix, the C -score is calculated as C ab ϭ ( R a -S )( R b -S ). Here R a is the row total for species a (number of occurrences of species a ), R b is the row total for species b , and S is the number of sites that contain both species a and species b . Relatively large values of C ab indicate that the species pair is segregated, with few or no sites that contain representatives of both species. Relatively small values of C ab indicate that the species pair is aggregated or nested, with many sites containing representatives of both species. Th e minimum value of C ab is zero, which occurs when at least one species always co-occurs with the other ( R a ϭ S or R b ϭ S ). Th e maximum value of C ab is ( R a )( R b ), which occurs when both species never co-occur ( S ϭ 0). Th e average C -score calculated for all possible species pairs has traditionally been used as a single metric of community-wide aggregation or segregation. However, a matrix with a significantly large C -score (segregated matrix) will contain both aggregated and segregated species pairs . For this reason, we fi rst tested each data matrix for signifi cance of the overall C -score, and then proceeded to determine which individual species pairs were signifi cantly aggregated or segregated.

Null models
C -scores for individual species pairs were compared to the statistical expectation for a set of stochastic matrices in which species interactions, environmental associations, and dispersal limitations are not important. Th ese matrices were not generated by a mechanistic model of species dispersal and colonization, such as the neutral model (Hubbell 2001). Instead, they were generated by a simple statistical randomization that captures the general eff ects of these processes without specifying a particular mechanistic model and its parameter values . We used a ' fi xed-fi xed ' null model (Gotelli 2000) in which the sample space consists of random matrices that have the same dimensions as the original matrix, the same percentage matrix fi ll, and the same row and column totals. Th e use of fi xed row and column totals preserves inherent heterogeneity that is commonly observed in the number of species occurrences (row totals), and the number of species observed per site (column totals). Preserving these features of the data matrix is especially important for paleoecological matrices because taphonomic biases can cause some species to be preserved more often than others, and can cause some sites or samples to yield more species than others (Behrensmeyer et al. 2000). By preserving these features in the original data matrix, the fi xed-fi xed algorithm ensures that only patterns of species co-occurrence above and beyond those generated by the margin totals of the matrix will be detected. Extensive benchmark testing for co-occurrence analysis (Gotelli 2000) and nestedness analysis  has demonstrated that the fi xed-fi xed algorithm has good type I error properties, and will infrequently reject the null hypotheses for heterogeneous matrices generated by simple stochastic sampling. At the same time, the fi xed-fi xed algorithm has reasonably good statistical power for detecting patterns of species segregation in structured matrices that have been deliberately seeded with random noise (Gotelli 2000).
We generated the C -score and implemented the fi xed-fi xed algorithm in the FORTRAN program Pairs (Ulrich 2008). Pairs uses an effi cient swapping algorithm, in which the elements of 2 ϫ 2 submatrices are randomly chosen and swapped. Th is algorithm creates a sequence of matrices that preserve row and column total sums. However, sequential matrices created by swapping this way are serially correlated and not truly independent. To break this dependence, we used an independent swap algorithm (Gotelli and Entsminger 2003), in which each random matrix is generated by a separate swapping sequence that always begins with the original matrix. To create each random matrix, we swapped n submatrices, where n equals 10 times the number of matrix elements (e.g. 10 ϫ (number of species ϫ number of sites)). A fresh sequence of swaps was used to generate each random matrix. One thousand random matrices were created this way for each empirical matrix, although consistent results can be obtained with as few as 100 random matrices (Ulrich 2008).

Testing individual pairs
In a standard null model analysis, a single community-wide metric (such as the average C -score for a matrix) is compared to the histogram of C -scores generated for a sample of random matrices (H o ), and the probability of the observed C -score|H o is estimated from the tail area of the simulated distribution (Manly 1995). Th e same approach to estimating p values can be applied to the C -score measured for each individual pair of species in the matrix. However, because there are m ( m -1)/2 such non-independent pairs of species in a single matrix (where m ϭ number of species in the matrix), the potential rate of false positives with so

Site classifi cations
For any particular pair of non-random species, each site can be assigned to one of four mutually exclusive classes, based on the presence of one species (1,0) or (0,1), both species (1,1) or neither (0,0) ( Fig. 1). Th e average characteristics of sites assigned to these 4 classes will generate additional patterns that can be used to distinguish among diff erent causes of non-randomness. If the site characteristics are measured as continuous variables (such as average annual precipitation), sets of sites can be compared with ANOVA or t-tests. If the site characteristics are measured as discrete variables (such as depositional environment for fossil materials), sets of sites can be compared with a 2-way contingency table analysis. In each case, the analysis will pinpoint whether characteristics of sites vary systematically based on the presence or absence of each species member in the pair.
For segregated species pairs, the critical comparison is between sites that have one species (1,0) and sites that have the other species (0,1) (i.e. allotopic sites). If species interactions are the critical factor in producing segregation, these two classes of allotopic sites should not diff er systematically in either their environmental characteristics or their spatial arrangement. Such species interactions might include pairwise competition or predation, but also might refl ect indirect eff ects of other species. For aggregated species pairs, the critical comparison is between sites that have both species (1,1) (i.e. syntopic sites) and sites that have neither species (0,0). If species interactions are important in producing aggregation, these syntopic and empty sites should not diff er systematically in either their environmental characteristics or their spatial arrangement. Such pairwise interactions might include pairwise mutualism or commensalism (or even predation), but also might refl ect indirect eff ects of other species.

Environment tests
For each site within a data matrix, we have diff erent measures of environment, either continuous (e.g. annual precipitation or temperature, as in this study; see Climate and distance data) or categorical (e.g. soil type or depositional environment). In the case of continuous measures, a one-way ANOVA can be used to compare the environment between the allotopic sites of segregated pairs ((1,0) vs (0,1)) and between the syntopic sites and empty sites of aggregated pairs ((1,1) vs (0,0)). Th e null hypothesis is that site characteristics do not diff er systematically between these pairs of site classifi cations. For the categorical measures, a two-way contingency table can be used to classify the sites. For the segregated species pairs, we counted the frequency of each environmental type for the two kinds of allotopic sites ((1,0) and (0,1)). For the aggregated species pairs, we counted the frequency of each environmental type for the syntopic sites (1,1) and the empty sites (0,0). For both kinds of two-way data tables, we used a chi-square test of association. Th e null hypothesis was that the frequencies of diff erent environmental types did not diff er among the site classes. If this null hypothesis is rejected, a parsimonious interpretation is that environmental associations are at least partly responsible for segregated or aggregated patterns of species occurrence (Fig. 1, 2). many statistical tests is quite high. A similar problem arises in the analysis of microarrays, in which the expression levels of thousands of potentially non-independent genes are assayed with parametric or non-parametric statistical tests (Kammenga et al. 2007). For null model analysis of the cooccurrence of individual species pairs , we adapted an empirical Bayes approach originally proposed by Efron (2005) for this problem of screening large numbers of non-independent tests. In brief, C -scores for each species pair are rescaled to a [0,1] range and binned into histogram categories. Next, the simulated data are binned in a similar way, and the mean and 95% confi dence interval of the C -scores of simulated species pairs in each bin is calculated. Finally, the original C -score values within each bin are ordered from smallest to largest C -scores. For the Pairs analysis, pairs of species are retained whose C -scores are above the simulated mean for the bin (Bayesian mean criterion), and which would be statistically signifi cant if the species pair was treated as an independent test. Th is ' double screen ' reduces some of the false positives that would arise by simply retaining all species pairs for which the uncorrected association (aggregated or segregated) yielded p Ͻ 0.05. For the bins that are near 0.0, these largest C scores will represent aggregated species pairs. For the bins that are near 1.0, these largest C scores will represent segregated species pairs. Th is is less conservative than a cut-point based on the 95% confi dence interval for bin deviations, but more conservative than an unadjusted count of signifi cant pairs, and usually more conservative than a sequential Bonferroni correction in which the pairs are ordered by their p values and a cutoff is imposed that is determined by both the individual p-value and its rank. Benchmark tests of the Pairs algorithm show that it is eff ective (though not perfect) at controlling for false positives while still allowing for detection of a relatively small subset of non-random species pairs from a binary presence -absence matrix . We ran the Pairs analysis for each data matrix to identify the subset of species pairs that exhibited strong aggregation or segregation.

Identifying the causes of non-randomness
Null model analysis has been a successful tool for identifying non-random patterns of species associations. But the analysis cannot, by itself, point to the causes of such segregated or aggregated patterns. Here we consider explicitly two major classes of mechanisms that might lead to non-random associations of species pairs: dispersal limitations and habitat or climate (environmental) fi ltering of species into groups with similar environmental niches. All signifi cant pairs that did not show signals of signifi cant environmental variation or dispersal limitation may provide evidence of a signifi cant species interaction, though it is also possible that environmental factors not considered in this analysis could contribute to nonrandom associations. To infer the roles of environmental factors and dispersal limitation, we move beyond the results of the standard null model tests with additional analysis of the characteristics of the sites. As described below, we focus specifi cally on subsets of sites that diff er signifi cantly in either environmental characteristics or spatial location.

Spatial overlap
In addition to environmental characteristics of the sites, we also often have geographic coordinates for each site. We treated the geographic coordinates as a two-element response vector, and then used a one-way MANOVA to compare the distance in coordinate space between the group centroids of the diff erent site types (e.g. Supplementary material Appendix 1, Fig. A1). Note that this analysis does not simply compare the geographic separation of sites with species A versus species B. Rather, for segregated pairs, it tests for the spatial overlap of allotopic sites ((1,0) and (0,1)), and, for aggregated pairs, it tests for the spatial overlap of syntopic (1,1) and empty (0,0) sites. Th e null hypothesis is that the diff erent classes of sites are not spatially segregated from one another. If this null hypothesis is rejected for segregated pairs, a parsimonious interpretation is that some form of historical dispersal limitation is at least partly responsible for the pattern of aggregated or segregated occurrence of a particular species pair. We decided to not impose a false positive correction to the set of spatial overlap and environmental tests, because we were already working with a subset of species pairs that are non-random, and did not wish to be overly conservative. However, this decision does aff ect the outcome of the classifi cation scheme (see Discussion).

Logical outcomes
Th e above framework leads to a forked logic tree with a total of 9 possible outcomes (Fig. 2). Th e fi rst level of branching is a tripartite division of all species pairs into random, aggregated, or segregated. Within the aggregated and segregated pairs, each fork can be further divided by whether or not sites show signifi cant patterns of spatial proximity, and then further divided by whether or not sites diff er in environmental characteristics (or vice versa; the order of division after the fi rst level of branching does not matter). Th e logic outcomes in the tree are not mutually exclusive. For example, ecologists working in the MacArthurian tradition would argue that competitive interactions can force inferior competitors into substandard environments and reduce spatial overlap, and that these kinds of interactions can scale up to the biogeographic scale. But for fossil data, whose total length spans timescales of centuries to millennia and beyond (and a temporal resolution ranging from decadal to multi-centennial or longer), it is more parsimonious to attribute such patterns to diff erences in environment, depositional environment, or dispersal limitation, if there is evidence for these processes. Th us, in this framework, species interactions are invoked only for cases in which species pairs are signifi cantly aggregated or segregated, but there is no evidence of environmental associations or non-random spatial associations (Fig. 2).

Fossil pollen data
To test our framework with paleontological data, we relied on fossil pollen data from the late Quaternary of eastern North America, which provide a spatially averaged description of vegetation composition with a resolution that can vary from ca 10 m to 10 km, depending on the site and vegetation (Prentice 1988, Webb 1993. Th ese data have been used previously for biogeographic analyses of community dissimi-larity (Blois et al. 2013a, b). Briefl y, the fossil pollen data are drawn from sites in the Neotoma Paleoecology Database ( Ͻ www.neotomadb.org Ͼ ), with site age models revised and standardized (Blois et al. 2011). Previous work indicates that the absolute temporal precision among sites in the dataset is about 500 yr throughout the latest Pleistocene, although relative temporal precision within sites can be on the order of decades to centuries (Blois et al. 2011). Th us, the data allow confi dent inferences of ecological responses to millennialscale climate change despite uncertainty in the base data. Overall, 527 sites contributed to the dataset, though not all sites were included at all times (Supplementary material Appendix 1, Fig. A2). Most sites were lacustrine environments, and the majority of the remaining sites were bogs or mires. Full details of site selection are given in Blois et al. (2011). We use the term ' taxon associations ' throughout the paper when referring to our pollen data because taxonomic resolution in this paper is usually but not strictly at the genus level. Many pollen types naturally correspond to genus resolution, though in some cases the taxa represent closely related and palynologically indistinguishable genera ( Ostrya and Carpinus , Juniperus and Th uja , Rumex and Oxyria , and Ambrosia -type pollen grains). Any available species-level data were aggregated to the generic level. Th e original fossil pollen data represent the relative abundance of each fossil pollen taxon at each site, calculated relative to the total upland pollen sum for the site, which included taxa identifi ed to taxonomic levels higher than genus. Relative abundances were then linearly interpolated to 1000-yr time steps, for every 1000 yr from 21 to 0 thousands of calibrated years before present (kyr BP). Th e program Pairs (Ulrich 2008), used for later analyses, considers any value greater than zero as ' present ' . Th is conversion may classify some taxa that are particularly prone to long-distance transport (e.g. Pinus ) as ' present ' when they were not actually present in the plant community surrounding the lake (Webb et al. 1981). Additionally, some of the rare taxa may be classifi ed as absent when they were present on the landscape. Th us, the inferred number of signifi cant taxonomic associations is likely to be an underrepresentation of the actual number of signifi cant pairs through time. In sum, the base data consist of a presence -absence matrix for all 106 fossil pollen taxa at all sites across eastern North America, for a snapshot of time every 1000 yr from 21 kyr BP to the present. Th ese 22 presence -absence matrices (rows ϭ genera, columns ϭ sites) and 22 site -environment matrices (rows ϭ sites, columns ϭ environment, described below) have been deposited in FigShare (doi: 10.6084/m9.fi gshare.841756).

Climate and distance data
To quantify environmental characteristics at each site, we relied on paleoclimate simulations from the National Center for Atmospheric Research Community Climate System Model ver. 3 (CCSM3) developed by Liu et al. (2009). CCSM3 provides transient paleoclimate simulations that are independent of the pollen data starting at 22 kyr BP, with seasonally averaged model outputs saved at a decadal time step (Liu et al. 2009). Th ese data were downscaled to a 0.5 ϫ 0.5 degree grid by Veloz et al. (2012) for every 1000 yr from 21 kyr BP to the present (available at Ͻ purl. org/climate Ͼ ). For each site and available time slice, we extracted the mean annual precipitation, mean winter (December -February) temperature, and mean summer (June -August) temperature. Because our primary concern was teasing apart the relative infl uence of climate or dispersal limitation on generating signifi cant taxon associations rather than looking at the infl uence of individual climate variables on pairwise associations, the three variables were fi rst scaled and centered, and then transformed to a principal components axis (PCA). We used the scores from the fi rst two principal components as measures of ' climate ' at each site. Th e fi rst principal components axis explained from 65.7 -84.1% of the variation, depending on the time period (Supplementary material Appendix 1, Table A1). While axis 1 had the strongest correlations with winter temperature, all three climate variables (precipitation, winter temperature, summer temperature) loaded relatively evenly at all time points (Supplementary material Appendix 1, Table A1). Th e second principal components axis primarily captured variation in precipitation and secondarily variation in summer temperature, and explained 13.9 -31.3% of the variation. Together, the fi rst two axes explain 93.7 -98.3% of the total climate variation and thus represent suitable proxies for the overall ' climate ' experienced at each site through time. Geographic coordinates are in North America Albers Equal Area Conic projection (origin: 40N, 96W, standard parallels: 20N, 60N).

Pairs analysis
We used the program Pairs (Ulrich 2008) to detect significant positive and negative taxon associations across space for each time period separately, implementing the Bayesian mean criterion for assessing signifi cance as described above.

Environment and dispersal analyses
Because climate variables were continuous, we used a MANOVA with the fi rst two principal component scores for each site as a two-element response vector to test for the infl uence of climate between the allotopic sites of segregated pairs ((1,0) vs (0,1)) and between the syntopic and empty sites of aggregated pairs ((1,1) vs (0,0)). Similarly, we used a MANOVA test with the geographic coordinates of each site to test for dispersal limitation, as described above. In this case, the test boils down to whether the geographic centroid of sites in one category (e. g. (1,0) for segregated pairs) is signifi cantly diff erent from the geographic centroid of sites in the other category (e.g. (0,1) for segregated pairs) (Supplementary material Appendix 1, Fig. A1). We then tallied the number and proportion of aggregated pairs and segregated pairs in four categories: 1) taxon pairs that show a signal of dispersal limitation (signifi cant value for the distance test); 2) taxon pairs that show a signal of climatic diff erences (signifi cant value for environment test); 3) taxon pairs that show signals of both dispersal limitation and climatic diff erences; and 4) taxon pairs that show signals of neither process. Th is fi nal category represents cases in which taxon associations can be more confi dently attributed to biotic interactions because there is no evidence of climatic eff ects and/or dispersal limitation.

Pairs
Over the past 21 000 yr, the number of unique taxon pairs ranged from 1176 (at 21 kyr BP) to 3570 (since 2 kyr BP). Table 1. Summary statistics for the Pairs analyses through time. From left to right, the columns show the time period of analysis (in thousands of calibrated years before present, kyr BP), number of sites in the dataset, the total number of taxon pairs, the proportion of random taxon pairs (relative to the total), the number of signifi cantly segregated and aggregated taxon pairs, and the proportion of segregated and aggregated pairs, relative to the total number of non-random pairwise associations.   Table 1 for sample sizes for each time period. (B) Attribution of the causes of signifi cantly segregated pairs, e.g. the relative proportion of signifi cantly segregated pairs that can be attributed to climate (blue), dispersal limitation (red), either or both processes (purple), or neither climate nor dispersal processes (e.g. a potential biotic interaction; gray). (C) same as (B), but for aggregated pairs.
In all time periods, nearly all taxon pairs were randomly associated (98 to 100% of possible taxon pairs; Table 1, Fig. 3), refl ecting the conservative nature of the Pairs tests for significantly associated taxon pairs. Of the non-random taxon pairs, the proportion of signifi cantly segregated versus aggregated pairs varied through time. For most of the past 21 000 yr, the aggregated pairs were more frequent than segregated pairs, but the proportion of segregated pairs steadily increased through Table 2. The proportion of segregated and aggregated pairs that can be attributed to climate, dispersal limitation, both, or neither process through time. Note that the total number of segregated or aggregated pairs per time slice may be different than the totals in Table 1 because statistical signifi cance could not be assessed for some taxon pairs.   Fig. 3).

Environment and dispersal analyses
Th e vast majority of sites with signifi cantly aggregated or segregated taxon pairs diff ered in climate or geographic distance, and usually diff ered in both aspects (Table 2, Fig.  3). Th ere were only a few cases in which neither climate nor geographic diff erences among sites were detected for signifi cant taxon associations. In such cases, we infer that biotic interactions may have played a role. Th e proportion of aggregated pairs that were potentially caused by biotic interactions increased further back in time (Table 2, (Table 3).

Discussion
Environmental fi ltering, dispersal limitation, and biotic interactions are three intertwined processes that have been commonly invoked to explain patterns of species associations across the landscape and through time (Tuomisto et al. 2003, Kraft et al. 2008. Th ey are often examined separately (but see Jim é nez et al. 2012), however, and no unifi ed framework exists for disentangling the three processes. Here, we provide such a framework and examine the role of these three processes in explaining spatial patterns of taxon associations through time for late Quaternary fossil pollen assemblages. While it may still be impossible to cleanly disentangle processes (e.g. when signals are consistent with both dispersal limitation and environmental fi ltering), our framework narrows down the sets of taxa for which biotic interactions are most likely to be important and leads to an explicit consideration of the kinds of patterns that can be used to distinguish environmental correlates from the eff ects of biotic interactions or dispersal limitation (Fig. 1, 2). Within the subset of 449 signifi cant pairs, we found that the potential role of biotic interactions in explaining signifi cant taxon associations at these spatial and temporal scales was minimal. Instead, the correlational patterns reported here are most consistent with climatic fi ltering or dispersal constraints across space as the key processes controlling generic coexistence of pollen taxa in eastern North America appear to be either climatic fi ltering or dispersal constraints across space; most often both processes were consistent with the data (Fig. 3). Th ese two processes have been commonly invoked to explain patterns of community structure or species changes Skov 2007, Blois et al. 2013a) through time. Th e pairwise approach taken here emphasizes that pairwise associations detected in assemblages of species may be transitory and are often tied to climatic or environmental diff erences, or spatial eff ects.

Patterns of signifi cant pairs through time
From 21 kyr BP to the present, between 98 and 100% of taxon pairs were randomly associated among sites (Table 1, Fig. 3). Th e large proportion of random taxon pairs is unsurprising. Th e method of testing individual pairs imposes a strong screen for type I error, but more importantly, most taxa occur in relatively few sites. For good statistical reasons, it is diffi cult to assert that segregated pairs are non-random when both members of the pair are relatively rare, though we can detect aggregation more easily in this case. But without any additional evidence, the most parsimonious interpretation of the observation that two rare species do not co-occur frequently is that the pattern is due to chance. Th us, our approach is an inherently conservative method to begin with, but avoids falsely attributing biological processes to patterns that are more parsimoniously accounted for by simple sampling properties of the data.
Of the non-random subset, there were more aggregated than segregated pairs in most time periods (Table 1, Fig.  3). Th ese results for fossil assemblages form an interesting contrast with a recent meta-analysis of pairwise associations in 272 presence -absence matrices for modern assemblages . In modern assemblages, most species pairs also showed random associations, although these tended to be concentrated in data matrices from a relatively small number of studies. However, the nonrandom fraction for modern assemblages was dominated by segregated species pairs, with a 4-fold excess of perfectly segregated checkerboard pairs compared to the most conservative null expectation. In the fossil assemblages examined here, the segregated taxon pairs show a trend of increasing frequency over the last 8000 yr, but exceed the frequency of aggregated taxon pairs only in the most recent (modern) sample. One reason for the overall dominance of aggregated versus segregated taxon pairs in these data matrices (relative to modern data matrices) may be that our analyses were conducted at the genus level, which means that taxa are more widespread across space than the individual species within each genus and thus more likely to be aggregated. Additionally, diff erential species richness within genera may aff ect the patterns of aggregation and segregation, and sensitivity tests regarding the magnitude of this eff ect would be an interesting avenue for future research with modern data. Nevertheless, these explanations do not explain the trend through time in the numbers of segregated and aggregated taxon pairs. Th is time series of matrices contains a sample size bias, with number of sites increasing in more recent assemblages. With fewer sites there will typically be fewer taxa present and/or detected, which obviously aff ects the absolute number of signifi cant taxon pairs that are identifi ed. Th ere is no obvious reason why diff erences in the number of sites should aff ect the relative frequency of segregated versus aggregated taxon pairs, however. One possibility is that this could relate to time-averaging within samples, especially if the older samples encompass more time than a sample taken from closer to the surface and thus artificially cause the appearance of aggregation. However, previous studies have shown that time-averaging is minimal in fossil sediment cores from lakes in eastern North America (Davis and Deevey 1964). Instead, the increasing proportion of segregated taxon pairs through time may also be a real pattern and refl ect the stronger zonation of plant communities along latitudinal gradients during the Holocene compared to the late Pleistocene (Williams et al. 2004).

Are signifi cant taxonomic associations due to biotic interactions?
Th e proportion of signifi cantly segregated or aggregated pairs that were potentially caused by biotic interactions was higher in the latest Pleistocene, during deglaciation, than in the Holocene (Table 2). Th ere were very few signifi cant taxon pairs during this time to begin with (Supplementary material Appendix 1, Table A2), however. Additionally, Blois et al. (2013a) noted that the time slices 21 -15 kyr BP were unreliable for inferring spatial dissimilarity patterns due to sample size diff erences through time (Table 1). Although there is no obvious reason that attribution of the cause of signifi cant pairs would be biased by the low number of sites in the earlier times within the dataset, we regard the results from 21 -15 kyr BP as tentative. If this is a real pattern, however, it could indicate increased importance of biotic interactions during a period of rapid environmental change and add to other work indicating that biotic interactions   (Table 1). If there is no color for a given time period (e.g. whitespace), the taxon pair was classifi ed as a random pair at that time or the cause could not be statistically attributed to any potential outcome. Colors are the same as in Fig. 3. See Supplementary material Appendix 1, Table A3 to match the index values with a taxon pair.
were signifi cant drivers of community dynamics (Gill et al. 2009(Gill et al. , 2012 under these circumstances. Only three time periods showed signifi cant taxon segregations that could not be attributed to either dispersal limitation or environmental diff erences. All of these were for time periods at or prior to 15 kyr BP -20, 16, and 15 kyr BP -and the segregations were transitory and involved diff erent taxa in each case (Table 2, 3; Fig. 4A). If these associations do refl ect the infl uence of biotic interactions, the interaction is transitory and does not persist through time.
If we had imposed a false discovery screen for both the environmental and spatial overlap tests, we might have found more pairs that would be classifi ed as examples of biotic interactions. However, most taxon pairs showed strong eff ects of both spatial and environmental segregation, so there would not be that many pairs classifi ed as nonsignifi cant by both tests after screening for false positives.
Taxon aggregations were also infrequent, but more common than segregations. In 12 of 22 time slices there were taxon pairs for which aggregations could not be attributed to climate or large-scale spatial overlap between syntopic versus empty sites. Th ese were most common between 16 and 11 kyr BP. Such pairs might refl ect positive biotic interactions such as direct mutualisms, indirect eff ects such as shared pollinators or exclusion from the same sites due to a shared competitor or predator, or unmeasured habitat or climatic associations. Th e occurrence of potential positive biotic interactions in the latest Pleistocene could also provide support for the stress-gradient hypothesis, which posits that positive interactions buff er species in physically harsh environments (Bertness andCallaway 1994, He et al. 2013).
Similar to segregated pairs, many of the aggregated taxon pairs were impermanent through time. However, in the case of signifi cantly aggregated pairs, particular taxa often were involved in signifi cant aggregations at several time slices (Table 3). For example, Bidens and Xanthium were significantly aggregated through time, and those diff erences were not attributable to climate or distance, from 14 -11 kyr BP. Th is association arises during the height of climate variability at the transition from glacial to interglacial periods, a time which also corresponds to the start of the ' no-analog ' period at many of the sites in eastern North America (Gill et al. 2012). Neither genus is considered to be characteristic of the ' no-analog ' period, however. Sample size is low for both Bidens and Xanthium (i.e. they are found at only a handful of sites in each time period), partly because both pollen types are often classifi ed as undiff erentiated Asteraceae. Taxonomic classifi cation should be consistent throughout each individual core, so interanalyst and intersite variability in taxonomic identifi cation should not bias one particular time period over another. Both genera often occur today in wetter sites, so they may prefer the same soils and hydrological environments (Chadde 2002) and this biotic association may refl ect an unmeasured habitat association (such as soil type, which was not included as a habitat variable in this study because we do not have appropriate data for most sites and times). Life history characteristics rather than biotic interactions may also play a role -both genera contain zoochorous plants, so perhaps the same suite of dispersers was contributing to their aggregated pattern. For both potentially unmeasured habitat preferences and life history characteristics, it is unclear why the aggregation should arise only from 14 -11 kyr BP. Later in the Holocene, the pair Vaccinium and Chamaedaphne (likely C. calyculata ) was signifi cantly aggregated at both 5 and 2 kyr BP. Although our framework classifi es this association as a biotic interaction, it more likely represents an unmeasured habitat association because many species in both Vaccinium and Chamaedaphne are found in peat bogs. Finally, the pair Sambucus and Hypericum represented a potential biotic interaction at two time slices (11 and 8 kyr BP), although in this case there are no obvious unmeasured habitat characteristics or life history traits that explain the signifi cant aggregation.

Conclusions
We have presented a parsimonious framework for analyzing species associations and inferring the roles of biotic interactions, habitat fi ltering, and dispersal limitation as causal agents. For fossil pollen assemblages of eastern North America, most taxon associations across space were random rather than signifi cantly aggregated or segregated. Th e relatively few positive and negative taxon associations in space could be attributed to climatic fi ltering and/or dispersal limitation, and not to direct biotic interactions. For those rare taxon pairs that potentially represented biotic interactions, the interaction was not persistent through time. Our identifi cation of biotic interactions is inherently conservative. Although environmental or geographic diff erences among sites were the primary processes associated with signifi cant taxon associations, biotic interactions were likely operating as well; they simply were not the dominant drivers of associations at these scales of space and time. Th is suggests the possibility that species interactions may be more important in shaping community dynamics on short-term rather than paleoecological time scales. Th ese methods can be applied to other assemblages of modern and fossil plants and animals for additional insights into the causes of species associations across a range of spatial and temporal scales.