Tree community data

Preindustrial forest composition was characterized using early land survey data aggregated to a 25 km2 grid across southern Québec. Early land survey data were extracted from logbooks reporting the original surveys of 302 townships between 1790 and 1900 (Supplementary Figure 1). Surveys were conducted along the boundaries of the townships and subdividing range lines within townships (every 1.6 km, Supplementary Figure 1), where surveyors described the forest composition in the form of lists of taxa (e.g., hemlock, beech and maples) (note that these data were collected in a different way than the General Land Office data widely used to reconstruct past forests characteristics in the northeastern United States52). In total, 103,011 lists of taxa from 1790 to 1900 were precisely georeferenced using historical and modern digital cadastral maps and used in this study (Supplementary Figure 1). These observations were divided into two geometric types: point observations (69% of total observations, usually spaced 100–200 m apart) and line descriptions (31% of total observations and usually pertaining to one lot boundary of ~260 m). In order to incorporate these two data types in the same database, a weight was assigned to each observation that represented the length of line observations and the mean spacing of point observations (the mean distance to the previous and next observations53). Data were then aggregated to a 25 km2 grid; in each grid cell we estimated a taxon’s frequency (Supplementary Figure 2) as the cumulative weight of observations in which the taxon was mentioned divided by the total weight of observations. A previous study validated the high accuracy of this method to estimate the nineteenth-century forest composition by comparison with the results obtained using early forest inventories available for a restricted portion of the study area54.

Modern forest composition was quantified using the Québec government’s forest inventories since 1980. These inventories are based on 400 m2 circular plots distributed across different types of productive forest through stratified random sampling (Supplementary Figure 1). Only plots that were located within a maximum distance of 3 km from historical observations were retained, resulting in a total of 59,359 plots. The vast majority of plots have been continuously forested during the last century, with <1% of plots having developed on abandoned farmlands. Within these plots, all stems >9 cm diameter at breast height of each species were measured and inventoried and were used to compute species basal area per plot. Surveyors of the nineteenth century did not systematically distinguish all phylogenetically close species and thus several species had to be grouped at the genus level to match the taxa that were recorded by historical surveyors (Picea spp., Pinus spp., Acer spp., Populus spp.). Because the nineteenth-century surveyors only specified the most abundant taxa in their observations, we also controlled the number of taxa per modern plot in order to obtain more comparable datasets. Taxa that represented <5% of the total basal area of a modern plot were removed; this threshold was chosen to obtain roughly the same number of taxa per observation in both historical and modern datasets (Supplementary Figure 1). Modern taxa frequencies were finally aggregated to the 25 km2 grid by calculating number of plots where a taxon is present divided by the number of plots in the corresponding grid cell (Supplementary Figure 2).

Taxa and community indices

Changes in functional composition were quantified using indices that characterize the affinities of taxa within communities with temperature (CTI), drought (CDTI), light (CSTI) and disturbance (CDI). These community indices were calculated as follows. First, indices that quantify taxon-specific relationships with temperature, drought tolerance, shade tolerance and disturbance were assigned to the main 17 taxa present in the historical and modern datasets (Supplementary Tables 1 and 2). The “taxon temperature index” (TTI) was calculated as the median annual temperature across the geographic range of a given taxon7,15,55 using interpolated annual surface temperature maps56 and continental scale distribution maps57. To assess the robustness of our analyses to the decision to use median temperatures, we also calculated variants of TTI using the mean, 10th percentile, and 90th percentile of temperatures across the geographic range of a given taxon (see Supplementary Methods; Supplementary Tables 4 and 5), with the latter two representing the cold and warm limits of distributions, respectively. Because none of these variants significantly changed the results (Supplementary Figure 7 and Supplementary Table 6), we report only results for the median in the main text. Taxon drought tolerance and shade tolerance indices (TDTI and TSTI, respectively, Supplementary Tables 1 and 2) were taken from a database of northern hemisphere trees and shrubs31,32. These indices were derived from a meta-analysis of published literature ranking taxa according to their drought and shade tolerance on a common scale (from low tolerance; 1, to high tolerance; 5); they have been frequently used to quantify tree functional composition32,58. The rankings obtained with our taxon climatic indices (TTI and TDTI) are also consistent with recent analyses of growth–climate relationships of restricted pools of taxa considered in our analysis37,38. Concerning disturbance-related indices, since TSTI mostly accounts for adaptions to severe stand-replacing disturbances (e.g., land clearing, clearcutting), we also developed for this study a more comprehensive taxon disturbance index (TDI) that considers adaptations to both severe and partial disturbances (e.g., partial cutting). TDI was calculated with a point system involving 11 key traits (Supplementary Tables 1) taken from various sources including the PLANTS and TRY databases59,60 and references on tree species ecology61,62. Five traits were considered primarily as adaptations to stand-replacing disturbances: shade tolerance, growth rate, longevity, age at sexual maturity and capacity of vegetative reproduction. Six additional traits were considered to reflect adaptations to both stand-replacing and partial disturbance: fruiting frequency, seed abundance, effective seed dispersal, germination substrate requirement, seedling vigor, and response to release. For each trait, the range of values across species was first standardized to vary from 0 to 1 (giving them equal weight) and then summed across traits for a given taxon to obtain TDI values.

For species grouped at the genus level (Picea spp., Pinus spp., Acer spp. and Populus spp.), values of TTI, TDTI, TSTI and TDI were averaged across the main species present in study area (our results were robust to various combinations of species included in these taxon-level values; see Supplementary Methods; Supplementary Tables 7 and 8, Supplementary Figures 8 and 9). For each time period separately, community indices were calculated for each 25 km2 grid cell as the weighted mean index value across taxa found in that grid cell, with weights determined by taxon frequency. In order to obtain values comparable across the four community-level indices, despite the different scales of taxon-level indices, all taxon-level indices were first standardized to vary from 0 to 1 among taxa prior to calculations of community-level indices (Supplementary Table 2). Thus, each community index potentially ranges from 0 to 1 (Supplementary Figure 3), and changes in community indices (the difference between the 1980–2010 and 1790–1900 periods) potentially range between −1 and 1.

Climate and population density data

Climate data were obtained at the 25 km2 grid scale for the 1901–1980 period using BIOSIM 10 software63. The year 1901 corresponds to the earliest date for which accurate climate records are available in our study area, while 1980 correspond to the earliest modern inventories used in the analysis. Changes in mean annual temperature were calculated as the linear slope of mean annual temperature per calendar year for each cell separately. To assess changes in moisture regime, we used the SPEI, which is based on the difference between monthly precipitation and temperature-induced potential evapotranspiration64. Negative values of SPEI indicate dry events, with values < −1.5 commonly considered as severe drought, while positive values indicate wet events. The 24-month SPEI between January 1901 and December 1980 were calculated in R using the SPEI package65 and changes in SPEI were calculated with the linear slope of SPEI per month for each cell separately.

In the absence of precise and extensive data on land-use history over the relevant spatial and temporal scales in our study area, we used changes in human population density since the early nineteenth century as a proxy for the intensity of anthropogenic disturbances over this time period. Sub-regions of our study area experienced their most intense periods of land use at different times. For example, some areas reached their peak in agricultural development and rural population as early as the end of the nineteenth century28, shortly after surveyors opened the territory for settlement (Supplementary Figures 1 and 5). Other regions reached such peaks between the middle and the end the twentieth century29. In order to account for this heterogeneity, we computed changes in mean population density between 1831 and, for each cell separately, the maximum mean population density recorded in one of the three subsequent censuses in 1871, 1951 or 2001 (see details in Supplementary Figure 5), which represent three key periods of land-use history in our study area28,29. Data from the 1831 and 1871 censuses have been made available by the Centre Interuniversitaire d’Études Québecoises (CIEQ; Projet GÉORIA, version 2003) and data from the 1951 and 2001 censuses are freely available on the Canada’s Century Research Infrastructure (CCRI) website (https://ccri.library.ualberta.ca).

Statistical analysis

Relationships between changes in community indices and potential predictors were tested with linear mixed effects models, which were constructed in R with the nlme package66. In all models, we accounted for natural biogeographic variation in forest composition by first conducting a spatially constrained clustering of 25 km2 cells based on preindustrial composition and then using the six resulting groups (see Supplementary Figure 6) as a random factor in all models. Spatial autocorrelation was taken into account within cluster groups using an exponential spatial correlation structure constructed with the corExp function contained in the nlme package66. We first tested for individual effects of potential predictors on community indices: we tested effects of ΔTemperature on ΔCTI, of ΔSPEI on ΔCDTI and of ΔPopulation on both ΔCSTI and ΔCDI. Second, because taxon-level affinities with disturbance, light and drought were correlated to some degree (Supplementary Table 3), we tested effects of potential predictors after controlling for changes in other community indices32. Effects of ΔTemperature on ΔCTI and of ΔSPEI on ΔCDTI were tested after accounting for changes in disturbance-related indices (i.e., ΔCSTI and ΔCDI). Conversely, effects of changes in population density on ΔCSTI and ΔCDI were tested after accounting for changes in climate change-related indices (i.e., ΔCTI and ΔCDTI). Finally, we tested correlations between changes in climate- and disturbance-related community indices, with ΔCTI or ΔCDTI as dependent variables (separately) and ΔCSTI and ΔCDI as predictors.