Spatiotemporal distribution of Holocene populations in North America

Significance We provide the first maps to our knowledge of spatiotemporal paleodemographic growth following human migration into the Americas for the past 13,000 y, using a statistical approach that simultaneously addresses sampling and taphonomic biases. The Canadian Archaeological Radiocarbon Database is sufficiently complete in many areas, demonstrating high correspondence between continental-scale 14C-inferred population estimates and generally accepted archaeological history. Increases in population density seem robust for eastern and western North America, as well as central Alaska and the region surrounding Cahokia. These results are the first step toward being able to understand continental-scale human impacts on the North American ecosystem during the Holocene as well as demographic growth and migrations in relation to environmental changes. As the Cordilleran and Laurentide Ice Sheets retreated, North America was colonized by human populations; however, the spatial patterns of subsequent population growth are unclear. Temporal frequency distributions of aggregated radiocarbon (14C) dates are used as a proxy of population size and can be used to track this expansion. The Canadian Archaeological Radiocarbon Database contains more than 35,000 14C dates and is used in this study to map the spatiotemporal demographic changes of Holocene populations in North America at a continental scale for the past 13,000 y. We use the kernel method, which converts the spatial distribution of 14C dates into estimates of population density at 500-y intervals. The resulting maps reveal temporally distinct, dynamic patterns associated with paleodemographic trends that correspond well to genetic, archaeological, and ethnohistoric evidence of human occupation. These results have implications for hypothesizing and testing migration routes into and across North America as well as the relative influence of North American populations on the evolution of the North American ecosystem.

D atabases of radiocarbon dates obtained from archaeological sites are used as a source of information about past demographic changes of Holocene populations (1)(2)(3)(4)(5). These are typically analyzed using summed probability distributions for a given region, which are interpreted as time series of site occupancy and inferred population tendency. For example, frequency distributions have been used to estimate continental-scale population growth in North America (North America here means the United States and Canada, the extent of our database) over the past 15 ka using the Canadian Archaeological Radiocarbon Database (CARD) (2,3). Regional studies in North America have highlighted the impacts of environmental change on population size using these methods and documented the impact of human activities on the vegetation (1,(6)(7)(8)(9)(10). Although some studies have attempted to investigate between-region demographic changes for particular time periods (11) or map the point distribution for certain time intervals (3), the spatiotemporal distribution of demographic growth through the Holocene in North America after the arrival of humans has not been tracked at a continental scale.
The underlying assumption of the "dates as data" approach is that the frequency of 14 C dates is proportional to population density. This relation improves when aggregated data have a sufficient sample size to apply methods to overcome archaeological sampling bias and taphonomic and calibration effects (12). Summed probability distributions have been used to assess paleodemographic trends in the Americas, Europe, Eurasia, and Australia (e.g., refs. 1, 2, 11, and 13-20) and methodological issues have been discussed at length (see refs. 17, 21, and 22 and Supporting Information for examples). However, to determine and compare demographic patterns in both space and time, a way is needed to estimate the density of archaeological records at a particular time.
We build upon previous work in North America by using an alternative method with data from the CARD and provide the first estimates to our knowledge of the spatiotemporal distribution of the population of North America during the past 13 ka. The kernel method is used to convert the spatial distribution of 14 C dates into estimates of population density at 500-y intervals with a kernel radius (or bandwidth, terms we use interchangeably) of 600 km, a bin size and radius that optimize the macrotemporal and -spatial patterning at continental and geological scales (Supporting Information). The radius of 600 km was chosen based on the optimal bandwidth (595 km) calculated following Scott's rule for bandwidth selection (23). This bandwidth does not have a simple physical meaning (such as the range of a hunter-gatherer); it is meant to represent where humans could have been (in a probabilistic sense) during the 500-y interval based on the distribution of dates, because each interval represents many generations of humans.
The CARD is unique because it contains 35,905 14 C dates, 33,756 of which have geographical, chronological, and descriptive information permitting spatial and temporal analysis. The frequency of dates increases exponentially through time, with a maximum frequency around 0.7 ka calibrated B.P. (all future references to dates are in thousands of years calibrated B.P., denoted simply as ka). Previous results (2), as well as our own cross-validation

Significance
We provide the first maps to our knowledge of spatiotemporal paleodemographic growth following human migration into the Americas for the past 13,000 y, using a statistical approach that simultaneously addresses sampling and taphonomic biases. The Canadian Archaeological Radiocarbon Database is sufficiently complete in many areas, demonstrating high correspondence between continental-scale 14 C-inferred population estimates and generally accepted archaeological history. Increases in population density seem robust for eastern and western North America, as well as central Alaska and the region surrounding Cahokia. These results are the first step toward being able to understand continental-scale human impacts on the North American ecosystem during the Holocene as well as demographic growth and migrations in relation to environmental changes. exercises (Supporting Information), indicate that this amount and distribution of dates are sufficiently representative to begin studying past populations at a continental scale (2).
Based on our method, we hypothesize that the resulting intensity distributions (referred to as radiocarbon frequency population estimates, or RFPEs) capture paleodemographic trends because (i) the RFPEs show clear patterns when methods to reduce taphonomic and sampling bias are applied and (ii) the patterns correlate well with preexisting archaeological interpretations of human cultural change across the continent that were derived partly from 14 C dates but also other data sources and inferences, thus representing a partial independent test. The results of this study can be used to model human demographic growth in relation to environmental change and cultural innovation and to determine human impacts on the landscape or fauna at a continental scale, particularly in relation to the early anthropogenic hypothesis (24).

Distribution of 14 C Dates
At the continental scale of this study, the distribution of 14 C dates is the result of human activity that is proportional to population density, although data incompleteness, sampling intensity, and natural taphonomic processes have an effect as well. There are too few dates in the CARD before 13 ka to create reliable population estimates before this time using the methods used in this analysis. Other known issues with the database include a variable temporal and geographic representation, and in particular a lack of data from the southern United States, which limits interpretation of results in this area.
As shown in Fig. 1A, some areas are extensively represented (e.g., Wyoming) and others less so (e.g., Texas and Florida). Areas with greater sampling intensity are commonly the result of heightened archaeological interest, mapping objectives, and available funding, whereas others relate to issues involved with the construction of the database, including variable effort in data gathering from available literature (3). This creates a spatial bias that can confound continental-scale maps of population density. To account for this bias and normalize the series of RFPE maps, a spatially varying sampling intensity was computed using 14 C dates from all sites. Using kernel density estimation, a sampling intensity map was then made that treated multiple dates from the same site in each time interval as one and standardized the influence of sites with higher sampling intensities ( Fig. 1B; details in Materials and Methods).
There is a tendency in the CARD for a greater number of more recent archaeological sites because geological deposits tend to disappear over time due to normal taphonomic processes (3). We appreciate that taphonomy may vary spatially but apply a spatially constant taphonomic correction to the maps (2,25) because the undertaking of how, where, and by how much this should be varied requires more study at regional scales. We propose that areas along some coastlines may require an even greater correction due to the loss of accessible archaeological sites as a consequence of sea level rise, or due to reduced sampling of higher elevation stranded shorelines. Similarly, for the Arctic, the kernel function radius of 600 km used to obtain the RFPEs (which is indicative of the spatial extent of the influence of a site on the surrounding intensity) may be too large due to geographic constraints (see Results, Canadian Arctic and Supporting Information).

Data Verification
The only time period where patterns produced using the CARD data may be verified with independent data are just before European contact, when ethnohistoric estimates of the spatial distribution of population have been made (26). Although the absolute numbers have been questioned (27), the relative demographic and spatial patterns should be broadly accurate (26). The RFPE patterns estimated from 0.4 to 0.6 ka ( Fig. 2A) generally correspond to those based on precontact estimates of population size (Fig. 2B), with half of the distinct regions in the maps differing by only one level of density. Across a band of the northern United States (and British Columbia), where we have high confidence in the data, the estimates from both sources are similar. In the south where the database is known to be incomplete, and in New England, our estimates of population density are much lower. In the north, an area where the ethnohistoric data were known to be incomplete, our estimates differ by only one or two levels. Despite calibration, taphonomic, or sampling effects, a broad-scale pattern emerges from the data and a reasonable hypothesis is that it is recording a demographic trend. This is the first verification of the accuracy of the CARD data to our knowledge and suggests that the current number of archaeological records in the database, although not exhaustive, is sufficient to produce a general representation of the demographic history of North America (2), if the regional incompleteness of the database is accounted for. Larger datasets and refinements to the method will enhance this kind of modeling, as will updates of the CARD itself.

Results
Alaska. RFPEs fluctuate steadily for the entire study period ( Fig.  3 and Fig. S1), which is realistic considering Alaska was the location of repeated migrations into the Americas (28). Estimates are relatively high in northern and central Alaska in comparison with other parts of North America between 7.5 and 13 ka, with a peak at 11.5 ka. This coincides with an abundance of sites associated with fluted points, characteristic of technologies used by Paleoindian groups after 13.5 ka (29). Starting around 6 ka, a second increase is centered on the Aleutian Islands and slowly moves east over the next 1.5 ka. This agrees with archaeological evidence suggesting Aleut peoples colonized the central islands ca. 6 ka (30). The stable but relatively low estimates in Alaska just before the Aleut occupation may be associated with a well- known "hiatus" in human presence (31). Beginning at 4.5 ka, a 1,000-y-long increase is observed, composed of two centers along the northern and southwestern coasts, coinciding with the suspected migration of early Paleoeskimo from Siberia around 4.5 ka (32) and the Arctic Small Tool tradition in the eastern Aleutians (33). Estimates show a decrease for several hundred years and then begin to increase at the northern coast at 2.5 ka and along the west coast (2 ka), until finally the majority of the Alaskan region seems occupied. The development of Thule culture in Alaska (beginning ca. 1 ka) coincides with the increasing RFPE patterns (32,34).
Canadian Arctic. The maps show a rapid increase in RFPEs in coastal northwestern Canada, including the Mackenzie Delta, Banks Island, and western Victoria Island around 6.5 ka, and an increase shortly thereafter on Ellesmere Island (Fig. S1). This predates the Independence I (35) and Early Pre-Dorset (36) occupations by several thousand years and may be due to known dating problems from Arctic sites (37) and temporal smoothing of the maps. However, these dates support recent genetic evidence suggesting there was human activity in the general region around 6,000 y ago (38). By 4 ka, there are signs of occupation across the entire Canadian Arctic, and at 2.5 ka, site density increases to the south. This is due to data from Dorset culture sites that have been found from Victoria Island in the west to Greenland and Newfoundland in the east (39). RFPEs on Victoria Island as well as the area surrounding Foxe Basin are likely underestimated, because Pre-Dorset archaeological evidence for the latter suggests a population was centered on this region around 2.5 ka (40). An increase in RFPEs along the southernmost Thule migration route from Alaska to western Hudson Bay beginning at 0.75 ka is expected in the context of history of the Thule people (41). A site in Nunavut at 9.5 ka, in an area that was glaciated at this time, causes an isolated population center in the central Arctic, which persists until 8.1 ka, at which point it is joined with the population further west (Fig. S1) and by which time the ice had retreated. The spatial extent of this population is a product of the kernel radius chosen (600 km); a larger radius would merge this population with its western neighbor (Fig. S2). Although these dates seem inconsistent with the ice extent (which is, however, based on radiocarbon data), it is reasonable to hypothesize that humans were found along the edge of the ice at this time.
British Columbia and the Ice-Free Corridor. Before and during the Last Glacial Maximum, human migration into the Americas via an interior route was not possible due to the presence of ice sheets. At ∼13.5 ka an ice-free corridor opened between the Laurentide and Cordilleran Ice Sheets (42,43). Genetic and cultural evidence suggests a nomadic pre-Clovis population would have been the first to migrate northward into the corridor from further south around 10.5 ka, following a previous migration south via a Pacific Ocean coastal route (29). In the maps, a high-density population is observed in the ice-free corridor two millennia before 10.5 ka (Fig.  S1). RFPEs remain high until 12 ka, when the Cordilleran Ice Sheet retreated and westward human expansion into this new terrain likely occurred. Around 11.5 ka, RFPEs increased along the entire west coast of British Columbia. This timing is in line with archaeological findings of maritime communities living on the outer coast of Haida Gwaii (44) and is consistent with the hypothesis that Pleistocene-Early Holocene sites in outer coastal regions would today be underwater and thus poorly sampled archaeologically. A mirror isostatic effect, in which glacial loading depresses areas under ice sheets while raising offshore areas beyond the ice margin, strands shorelines along the mainland at much higher elevations, which leads to undersampling (45). At 7 ka, a second increase at the coast is observed, followed by an increase in the corridor. After 5.5 ka, populations seem to increase in both locations, potentially due to the stabilization of relative sea level on the coast and the development of new trap and tool technologies east of the Rocky Mountains (46). In the Late Holocene, the entire province of British Columbia was densely populated and remained so until 1 ka, reflecting the demographic success of complex hunter-gatherer cultures.
Eastern United States. Beginning around 11 ka, there was an increase in RFPEs in the southeastern United States, which grew larger and expanded northward then westward after 9.5 ka (Fig.  S1). This interpretation is plausible in light of numerous archaeological studies indicating a strong presence of Paleoindians in the southeast before 10 ka (47). RFPEs decreased east of the Appalachians after 9.5 ka and were centered in Missouri at 8.5 ka, where the presence of Paleoindian populations has been confirmed (48). Estimates fluctuated until 4.9 ka when populations grew to the east and west of the Appalachians as well as in the Middle Atlantic and New England regions.  but the dates from this site are consistent with one another and are accepted in archaeological literature (50). The RFPEs decreased following this, as temperatures and dryness increased during the Younger Dryas (51). RFPE values fluctuated and began to increase around 10.5 ka, consistent with archaeological evidence from Debert (50). At 8.5 ka, RFPEs tripled in Newfoundland and Labrador (Fig. 3) with the onset of early Maritime Archaic settlement along the coast [5.5-8 ka (52)]. Estimates are relatively low in the Maritime Provinces because many sites were presumably inundated by sea level rise. There is a slight increase in RFPEs during the Middle Archaic in eastern Quebec and central Labrador, and at 4.5 ka there was a migration into these provinces. A larger population seems to have spread across eastern Hudson Bay and the coastline of Labrador between 3.2 and 4 ka during the late Archaic. The influence of Arctic-adapted Paleoeskimo populations (e.g., Groswater and Dorset sites), the Meadowood complex, and Late Woodland populations are observed during the last 3 ka, although these values may be overestimated due to a higher number of sites being divided by a low sampling intensity in more recent times. RFPEs decreased in central Quebec but remained high in Newfoundland, Labrador and the area surrounding New Brunswick until 1.5 ka, when estimates increase over Atlantic Canada.
Central and Western United States. The oldest maps contain two centers of relatively high estimates in Arizona and on the Colorado-Kansas border, areas where paleoecological records confirm the presence of Paleoindians (53). At 12 ka, these centers shifted to Texas, although this area remains anomalous for the entire record due to the relatively large number of dates in the Earlyand Mid-Holocene compared with the more recent past (Fig. S1). By 10 ka, an increase appeared in California and persisted until 7.5 ka, which is likely associated with offshore communities exploiting marine resources (54). Also during this time there was a slight increase in RFPEs in Montana, which became joined with larger populations living in the ice-free corridor. The latter half of the Holocene was characterized by an intensification in Idaho, which shifted toward the coast at 4 ka.

Discussion
If known issues with the database are accounted for, RFPEs correspond with archaeological, paleoenvironmental and genetic evidence for human migration and dispersal into and across North America (Table 1). The maps (Fig. 3, Fig. S1, and Movie S1) show the importance of proximity to the coast for human settlement as well as areas where human activity seems to have been continuous throughout much of the Holocene (northeastern North America, Alaska, and western Plains). The maps also show examples where sampling bias (i.e., the quality of the CARD data) affects the RFPEs. For example, the Kivalliq Region (west of Hudson Bay) and the Ungava Region show an overestimation beginning ca. 4 ka. This is connected to selective archaeological work in the region and an extensive research program into the local precontact history (56)(57)(58) combined with comparatively few dates in the CARD before this interval. We have shown that aggregated continental data exhibit macrodemographic patterns comparable to known regional histories; in addition, use of the CARD permits a continental-scale analysis. At 12.5, 9, 4.5, and 1 ka, the overall pattern is close to generally accepted culture-historical reconstructions (29,41,44,49,52,54). There was an overall increase in continental population estimates from 12.5 ka, as expected, followed by a decrease between 0.5 and 1 ka due partly to the taphonomic correction, but also reflecting the decreased use of 14 C dating in the context of historical methods and a lessened archaeological interest in recent periods. The distribution at ca. 0.8-1.2 ka provides a broad example of both current and ancestral population distribution because the majority of currently understood large-scale migrations had occurred by this period, and comparatively little migration has occurred since, with the exception of the Athapaskan migration into the southwest.
These results illustrate the value in applying advanced statistical methods to aggregate 14 C data from archaeological databases. The accumulation of known 14 C data into facilities such as the CARD should be a priority, because this offers an opportunity to summarize previous research and provide a context for regional and local studies. "Dates as data" (7, 12) is a relevant means of investigating and modeling demographic, and thus historical, trends across long periods of time and continental areas. The kernel density method, in combination with statistical approaches to reduce sampling and taphonomic bias, holds potential for accommodating known issues in 14 C datasets that cannot be handled by other statistical analyses developed so far (e.g., averaging the number of observations). These results suggest that the CARD is a highly useful archive of paleodemographic data that can be used to investigate subjects such as migration routes into and across North America as well as a valuable tool for studies linking anthropogenic impacts with post-ice age faunal extinctions, ecosystem decline, and changing environmental and climatic conditions.

Materials and Methods
The CARD (3) contains 35,905 radiocarbon dates derived from cultural and paleoenvironmental material collected from 9,149 geographically distinct archaeological sites. The selection of dates for this study depended on completeness of the database (locations for 77 dates are missing), location (61 dates are from Russia), dating information (2,071 entries are missing normalized 14 C ages, normalized 14 C errors, or both), and classification (cultural, paleoenvironmental, etc.). Descriptions of the cultural association and type of material dated and general comments about the entry were also considered, although dates were not eliminated from the study if this information was missing. This resulted in 29,609 cultural dates that we used to create the RFPE maps and 4,087 dates from paleoenvironmental (or unknown) material (Fig. 1A). The 14 C ages calculated from the CARD (1) have been converted to calendar years (based on CALIB v.5.0.2 and IntCal04) using the median of the 2 sigma ages of the calibration curve (3). Because the maps are based on intervals of 500 y, any error when using the median as an estimate is much less than the sample interval.
To reduce sampling bias in the CARD, which is a result of inconsistencies in regional-(e.g., Canada vs. United States) and local-scale (e.g., Wyoming) sampling strategies and intensities, a sampling intensity map (Fig. 1B) was created using kernel density estimation. Only geographically distinct sites (n = 7,754) were used to create Fig. 1B; the amount of sampling (intensity) is estimated by considering only the geographic location of each site so that multiple dates at the same site have no influence. All archaeological sites were used in this step regardless of the age of the dated material or its classification (i.e., cultural, paleoenvironmental, or unknown). This map is interpreted as the density of sites, in other words, the distribution of sampling sites in space, and not the number of 14 C dates per site. If this remained unaccounted for, regions with higher sampling intensity would always be associated with elevated population estimates.
The underlying assumption of this approach is that a greater number of dates is indicative of a higher population density (3,59). To produce the RFPE maps, a grid of points every 0.1°between 23-85°N and 173-51°W was created to define the area of interest (North America) in ESRI ArcMap 10.1. This grid, in combination with ArcGIS shapefiles of the extent of glacial ice at 1,000-y intervals (60), was used to differentiate between glaciated and inundated areas and land. The period from 0.5 to 13 ka was chosen due to the lack of CARD dates before and after this date range. To visualize patterns in older time periods we would need to use a larger smoothing radius, but this tends to overgeneralize the resulting patterns (Fig. S2). The data were examined in 500-y intervals as a compromise to avoid aggregating too many dates and having too few to produce a reliable estimate, with consecutive intervals overlapping every 100 y.
To temporally smooth the data within each interval and account for 14 C dating errors, dates that do not occur directly within the interval are considered if they occur within 400 y before or following the interval. The dates within an interval are given a weight (w) of 1, whereas the dates within ±400 y of an interval are given a weight of <1. The weighting is achieved by taking a date that occurs outside of an interval and computing the temporal distance to the interval boundary and dividing the result by 400 (i.e., dates that occur at a greater distance from an interval boundary will have a smaller weight). Dates within 400 y before or following an interval are considered because they can still have an influence on population intensity (i.e., the estimated density of people at a given time and location is influenced by how many people were there before them) and to account for possible dating errors. When multiple dates exist for a single site within an interval, each date is considered separately.
Kernel density estimation, a method previously applied to archaeological data from Europe (61,62), is used to produce the RFPE maps. An Epanechnikov kernel with a bandwidth fixed at 600 km (Fig. S3) is chosen, although the Gaussian kernel method and adaptive bandwidth options were also examined (Fig. S4). First, raw population estimates are produced. Generally, these estimates are highly correlated to the sampling intensity because a greater sampling effort produces a larger number of dates, which results in higher population estimate (2). Therefore, producing RFPE maps based on raw estimates does not truly reflect human activity due to this sampling bias. A more reliable estimate is obtained when the raw population estimates are divided by the sampling intensity ( Fig. 3 and Fig. S1). In doing so, RFPEs can be interpreted as the relative frequency of dates per site, which should better reflect population intensity than only considering the number of dates. A similar approach has previously been used to account for a temporal bias (2). The values on the maps have been rescaled to 1 and are plotted in a logarithmic scale, with each color being approximately three times that of the previous color. This is done to maintain a consistent scale across all time intervals, and to make temporal changes more apparent, especially during earlier intervals. The resulting maps were corrected using a taphonomic formula that had previously been applied to the CARD data (1,25).
To assess the quality and completeness of the CARD, we ran two additional analyses (Figs. S5-S8). In the first, individual dates within the database flagged as "anomalous" (i.e., too young or old) were removed (n = 1,419, or 5% or the original data points) and a new series of maps were made. This is representative of a worst-case scenario in which the dates with the highest probability of negatively influencing the RFPEs are removed, thus addressing the potential effect of "bad dates" on the quality of the maps. In the second analysis, a new series of maps is made after 50% of the dates (n = 17,953) had been randomly selected and removed, to determine whether or not the amount of data in the CARD is spatially representative. Figs. S6 and S8 show a very high degree of similarity between the original RFPE maps and the two new series despite the missing data, putting to rest our concerns regarding the state and usability of the database and its applicability in spatiotemporal analyses. A supplementary test was performed using two subsets of calibrated ages based on the IntCal04 and IntCal13 calibration curves (Fig. S9). This test was done to determine if changes made between IntCal04 and IntCal13 would have an effect on the calibration results. The resulting calibrated ages were nearly identical allowing us to use a previously calibrated version of the CARD based on CALIB v.5.0.2 and IntCal04 (1-3) (see Supporting Information for additional discussion).