Maritime Living around the North Pacific Rim: a brief history

The North Pacific Rim (Fig. 1) is home to Indigenous and mixed settler communities with deep histories of interdependence with marine ecosystems. In the past, Indigenous communities from Hokkaido to the Gulf of Alaska subsisted on many of the same marine resources, captured in broadly similar ways using darts, spears, harpoons, hook and line, and nets (Fitzhugh and Crowell, Reference Fitzhugh and Crowell1988; Yamaura, Reference Yamaura1998; McCartney and Veltre, Reference McCartney and Veltre1999; Fitzhugh and Dubreuil, Reference Fitzhugh and Dubreuil1999; Steffian et al., Reference Steffian, Saltonstall and Kopperl2006). The smaller North Pacific islands mostly lack terrestrial resources, whereas on the larger islands and mainland coasts people supplemented marine and aquatic foods with terrestrial harvests of deer, caribou, moose, and small mammals captured with arrows and darts, spears, snares, and traps.

Figure 1. Map of N. Pacific showing region included in the current analysis. Previous comparisons (Fitzhugh et al., Reference Fitzhugh, Brown, Misarti, Krupnik and Crowell2020) focused only on regions 6 (Aleutians), 8 (Kodiak), and 9 (Kurils). The current analysis expands to neighboring regions of Alaska (1–8) and Hokkaido+Kurils (9–14).

Archaeologically, the evidence for specialized maritime technologies and practices appears in the Holocene, though indirect evidence suggests possible maritime living during the last glacial maximum as well (Fitzhugh, Reference Fitzhugh, Cassidy, Ponkratova and Fitzhugh2022a). Coastal settlement is evident by the early Holocene in the subarctic Northeast Pacific (southern Alaska), later in the subarctic Northwest Pacific (Sea of Okhotsk), and finally in the Arctic (Bering and Chukchi Sea) regions (Fitzhugh, Reference Fitzhugh, Friesen and Mason2016, Reference Fitzhugh, Cassidy, Ponkratova and Fitzhugh2022b) (Fig. 2).

Figure 2. Representative culture history sequences of North Pacific archaeological cultures with maritime components found in the regions included in this analysis.

For millennia, communities in the Gulf of Alaska and Aleutians designed hand-crafted skin boats to navigate stormy and hypothermia-inducing waters. They harvested sea mammals as substantial as gray and humpback whales and fished for cod, halibut, sculpin, and greenling in nearshore environments and on the open shelf, tens of kilometers offshore. They caught salmon and other anadromous fish by the thousands with spears, weirs, nets, and traps, and they harvested shellfish, octopus, seaweeds, and other vital resources in the intertidal zones (Knecht and Davis, Reference Knecht, Davis and Dumond2001; Fitzhugh, Reference Fitzhugh2003, Reference Fitzhugh, Friesen and Mason2016; Corbett et al., Reference Corbett, West and Lefevre2010; West et al., Reference West, Hatfield, Wilmerding, Lefevre and Gualtieri2012; Steffian et al., Reference Steffian, Saltonstall, Yarborough, Friesen and Mason2016).

Unangax and Sugpiaq communities have continued to harvest traditional foods through the period of Russian and American colonialism. They have also become vitally integrated into the expanding commercial fisheries markets for the export of salmon, cod, pollock, halibut, crab, and an increasing range of other products. As they have for over a century, Native and non-Native fishing crews continue to operate multi-generational fishing businesses out of villages and towns around the Gulf of Alaska and Aleutians (Reedy-Maschner, Reference Reedy-Maschner2004; Carothers, Reference Carothers, Colombi and Brooks2012; Reedy, Reference Reedy2020; Carothers et al., Reference Carothers, Black, Langdon, Donkersloot, Ringer, Coleman, Gavenus, Justin, Williams, Christiansen and Samuelson2021). The marine ecosystem provides food security through subsistence harvests to island and coastal villages, while commercial operations ensure that families and communities have access to cash-based essential resources and services.

On the other side of the North Pacific, the subarctic waters of northeast Hokkaido, Sakhalin, mainland coasts of the Okhotsk Sea, the Kuril Islands and Kamchatkan peninsula are occupied today by a mix of settlers (mostly Japanese, Russians, and Koreans) and Indigenous Ainu, Orok, Nivkhi, Even, Koryak, and Itel'men descent groups (Fitzhugh and Crowell, eds., 1988). In comparison to subarctic Alaskan coasts, the Asian Pacific Northeast shores experience stronger seasonal climate swings, cooled in the winter by bitter winds from Siberia and the circulation of chilled waters and pack ice along the eastern coast of the Bering Sea and from the northern Sea of Okhotsk to the shores of Hokkaido. Whale, walrus, seal and sea lion, salmon, cod and other marine fish, shellfish, and macroalgae were consumed along with terrestrial game and plants, where available. Historically, the Ainu used their intimate expertise of wild harvesting to supply products such as sea mammal furs, hawk and eagle feathers, and fish for trade to central Japan (Hudson, Reference Hudson1999). Okhotsk people also became active intermediaries in the trade of commodities between the wealthy states of China and Manchuria on the one hand and Japan on the other (Hudson, Reference Hudson1999; Walker, Reference Walker2001; Tezuka, Reference Tezuka2009). Indigenous and non-Indigenous communities from Hokkaido and the Russian Far East participate in extensive commercial fish, shellfish, and seaweed harvesting today, much of it bound for the expansive Japanese market as well as global consumption.

In contrast to the subarctic zones of Alaska and the southern Okhotsk Sea, the seasonally frozen northern Okhotsk, and around the Bering and Chukchi Seas, maritime economies were developed no earlier than 4000–5000 years ago (Tremayne and Rasic, Reference Tremayne, Rasic, Friesen and Mason2016), becoming specialized on marine mammal hunting and fishing after 3500 yr (Dumond and Bland, Reference Dumond and Bland1995; Fitzhugh, Reference Fitzhugh, Friesen and Mason2016, Reference Fitzhugh, Cassidy, Ponkratova and Fitzhugh2022b; Gusev, Reference Gusev, Cassidy, Ponkratova and Fitzhugh2022). On these coasts, people have since hunted ice-adapted seals, walrus, and whales from land-fast and pack ice, cracks in the ice (leads), and open water. They harvested (and continue to harvest) Arctic char and salmon from local rivers, terrestrial mammals (especially reindeer/caribou), and migratory birds as part of seasonal subsistence cycles. These animals were major staples for those living on the Bering and Chukchi Sea coasts, in the Brooks Range and up the rivers like the Anadyr, Kobuk, Noatak, Yukon, and Kuskokwim. Throughout these regions, mixtures of terrestrial and aquatic foods continue to be important for subsistence and food security of Iñupiat and Yupiit. In the post-contact era, commercial fisheries are much more recent and limited compared to the regions around the Gulf of Alaska and Aleutians.

As we see, residents of the North Pacific Rim adapted to and developed varying degrees of dependence on marine food resources through the Holocene. Those living in subarctic island environments were most fully dependent on the marine environment for their livelihoods and have the deepest histories of Indigenous fisheries, but even those communities living in the interior zones were linked to marine ecosystems through their reliance on anadromous fish like salmon and char. Where and when terrestrial game such as caribou, moose, deer, and migrating birds were plentiful, economies were less tethered to the productivity of the oceans. We expect populations closer to the coasts and especially on the islands of the subarctic Kodiak, Aleutians, and Kurils to have been most intimately influenced by changing marine ecological conditions.

Synoptic Feedbacks in the Climatology, Oceanography, and Marine Ecology of the North Pacific Subarctic Gyre and Marginals Seas Today

In the realm of natural resource policy, memories tend to be short and reference intervals for federal and state resource management even shorter, often based on trends observed over only a few years or decades. This leads to the handicap known as the “shifting baseline syndrome” (Jackson and Alexander, Reference Jackson, Alexander, Bolster, Chavez, Cournane, Erlandson, Field, Hardt, Lange, Leavenworth, Lotze, MacCall and McClenachan2011; McClenachan et al., Reference McClenachan, Ferretti and Baum2012; Klein and Thurstan, Reference Klein, Thurstan, Máñez and Poulsen2016; Eckert et al., Reference Eckert, Ban, Frid and McGreer2018). Historical analyses have revealed cycles of ecological variability otherwise missed in the annual to decade timeframes. Longer interval data, from multi-decadal to centennial to millennial could better inform managers on the ranges of variability that need to be considered in developing scenarios of the future that we know will not be similar to the past century.

At the shorter of these time scales, it is already well established that coupled North Pacific climate and ocean dynamics affect marine food webs and commercial harvests at multi-decadal scales. For example, the Pacific Decadal Oscillation (PDO) index tracks spatial anomalies in sea surface temperatures (SST) across the North Pacific (Schneider and Cornuelle, Reference Schneider and Cornuelle2005; Newman et al., Reference Newman, Alexander, Ault, Cobb, Deser, Di Lorenzo, Mantua, Miller, Minobe, Nakamura and Schneider2016). With annual records going back to the early 20th century, the PDO index reveals a multi-decadal oscillation in dominant and spatially discriminated warm and cold regimes. Importantly for our purposes, the PDO pattern was discovered through an effort to explain inversely correlated trends in salmon catch records from south to north along the North American West Coast (Francis et al., Reference Francis, Hare, Hollowed and Wooster1998; Mantua and Hare, Reference Mantua and Hare2002). When temperatures were lower than the century average for most years of a decade (or more) along the West Coast of North America, salmon catches increased for fishers in Washington, Oregon, and California but declined in the Gulf of Alaska. Catch trends reversed in warmer years.

The PDO is just one index used to understand complex climate-ocean-ecology variation in the North Pacific (Furtado et al., Reference Furtado, Di Lorenzo, Anderson and Schneider2012), and recent synthesis studies have concluded that PDO/PDO-like variability arises from a combination of drivers, situated both within the North Pacific and in the tropics (Newman et al., Reference Newman, Alexander, Ault, Cobb, Deser, Di Lorenzo, Mantua, Miller, Minobe, Nakamura and Schneider2016). Locally, North Pacific multi-decadal SST variability is closely connected to the dynamical strength and trajectory of the winter Aleutian Low (AL) pressure system. The AL is in turn influenced by the differential of continental vs. oceanic cooling in winter and atmospheric anomalies elsewhere in the northern hemisphere circulation system, including the El Niño Southern Oscillation (ENSO) and the East Asian Jet Stream (EAJS) (Chavez et al., Reference Chavez, Messié and Pennington2011; Newman et al., Reference Newman, Alexander, Ault, Cobb, Deser, Di Lorenzo, Mantua, Miller, Minobe, Nakamura and Schneider2016; Nagashima et al., Reference Nagashima, Addison, Irino, Omori, Yoshimura and Harada2021). PDO-like variability manifests as spatially patterned anomalies of warm and cold surface water across the North Pacific that reverse sign (flipping warm for cold and cold for warm) at some approximate interval. Accordingly, decades that are anomalously warm off the coast of western North America are anomalously cold around Japan, and vice versa (Chavez et al., Reference Chavez, Messié and Pennington2011). These differences are correlated with the productivity of commercial fisheries (e.g., salmon, anchovies, sardines) on both sides of the North Pacific basin. Chavez and colleagues (Reference Chavez, Ryan, Lluch-Cota and Ñiquen2003) see these patterns resulting from climate forcing effects on bottom-up processes tied to currents and upwelling, nutrient supply, primary production, and trophic amplification.

Because detailed, formalized scientific data collection started in this region only in the early 20th century, understanding of these patterns is limited in time and focused on seasonal, annual, and decadal variability and in the context of intensive commercial harvesting (Pennoyer, Reference Pennoyer1979; Hilborn, Reference Hilborn2012). A century of data on North Pacific ecological fluctuation is insufficient to clarify the potential range of variability that the integrated North Pacific marine ecosystem is capable of, especially under rapidly changing climate conditions. A variety of analyses have explored the time depth of North Pacific climate variability using proxy records such as tree rings (Biondi et al., Reference Biondi, Gershunov and Cayan2001; MacDonald and Case, Reference MacDonald and Case2005; D'Arrigo et al., Reference D'Arrigo, Mashig, Frank, Wilson and Jacoby2005; D'Arrigo and Wilson, Reference D'Arrigo and Wilson2006; Yasuda, Reference Yasuda2009; Barron and Anderson, Reference Barron and Anderson2011). They indicate that longer interval climate fluctuations have altered the intensity of decadal-scale (PDO-like) variability over century to millennial time scales, with implications for the spatial evolution of climate patterns in the past (e.g., Barron and Anderson, 2010; Nagashima et al., Reference Nagashima, Addison, Irino, Omori, Yoshimura and Harada2021).

Despite evidence linking climate and ecology in the late 20th century, the factors controlling ecosystem change in the North Pacific remain poorly understood (Beamish and Bouillon, Reference Beamish and Bouillon1993; Hare and Francis, Reference Hare and Francis1995; Hare et al., Reference Hare, Mantua and Francis1999; Overland and Stabeno, Reference Overland and Stabeno2004). In particular, the complexity of the relationship between long- and short-term ecosystem dynamics are highlighted by recent paleoecological research showing that sockeye salmon (Oncorhynchus nerka) populations in the Gulf of Alaska increased over century to millennial scales in times of enhanced storminess and anomalously cool sea-surface temperatures (SSTs; as during the Little Ice Age) and decreased in warmer time periods (such as the Medieval Climate Anomaly) (Misarti et al., Reference Misarti, Finney, Maschner and Wooller2009; Finney et al., Reference Finney, Alheit, Emeis, Field, Gutiérrez and Struck2010). This is in direct contrast to 20^th century PDO–salmon relationships. Such studies make clear that the historical record of ecosystem response to climate variability is too short, and that the paucity of long-term data impedes our ability to forecast future ecosystem changes or to discriminate between the compounded effects of natural variability and anthropogenic influences. Fortunately, the record of human interaction with marine ecosystems goes back thousands of years, giving us access to a deeper archive of information about both ecosystems in the past and human interaction with those systems (Misarti et al. Reference Misarti, Finney, Maschner and Wooller2009).

Synoptic Changes based on Holocene Proxies

Figure 3 shows long-term trends in proxies of climate, ecology and human populations based on marine, glacial, and lake cores, as well as archaeological radiocarbon accumulations. The atmospheric proxies depict century- to millennial-scale modalities in the North Pacific circulation from both sides of the North Pacific basin (Fig. 3a and 3b). These changes appear to correspond with late Holocene proxy Pacific salmon abundance trends from Kodiak Island (Fig. 3c; Finney et al., Reference Finney, Gregory-Eaves, Douglas and Smol2002) and human populations estimates (Fig. 3d; Mann et al., Reference Mann, Crowell, Hamilton and Finney1998; Fitzhugh et al., Reference Fitzhugh, Friesen and Mason2016, Reference Fitzhugh, Brown, Misarti, Krupnik and Crowell2020), which also appear to have fluctuated at century and millennial scales around the North Pacific Rim. Together, these data are consistent with the possibility that human populations were sensitive to fluctuations in climate and food availability. They also suggest that millennial-scale variability in marine resources has significantly higher amplitudes than 20th century decadal shifts would imply (Hare et al. Reference Hare, Mantua and Francis1999; Finney et al. Reference Finney, Gregory-Eaves, Douglas and Smol2002; see also Maschner et al., Reference Maschner, Finney, Jordan, Misarti, Tews, Knudsen, Maschner, Mason and McGhee2009b).

Figure 3. (a) Oxygen isotope record from NE Pacific (ice and lake cores): precipitation and Aleutian Low proxies (JBL from Anderson et al., Reference Anderson, Abbott, Finney and Burns2005; MTL from Fisher et al. Reference Fisher, Osterberg, Dyke, Dahl-Jensen, Demuth, Zdanowicz, Bourgeois, Koerner, Mayewski, Wake and Kreutz2008). (b) ESR dates from a Sea of Japan marine core: reflecting shifts of the East Asian jet stream and related to AL strength and direction (Nagashima et al. Reference Nagashima, Tada and Toyoda2013). (c) δ¹⁵N from Karluk Lake sediment core: proxy for sockeye salmon spawning population (Finney et al. Reference Finney, Gregory-Eaves, Douglas and Smol2002). (d-g) Archaeological radiocarbon temporal frequency distributions (TFDs), not taphonomically corrected (compare to Fig. 4), but cleaned to reduce sampling bias (Fitzhugh et al., Reference Fitzhugh, Brown, Misarti, Krupnik and Crowell2020; Maschner et al., Reference Maschner, Betts, Cornell, Dunne, Finney, Huntly, Jordan, King, Misarti, Reedy-Maschner and Russell2009a). See Table 1 for sample sizes and text for discussion.

Table 1. Archaeological radiocarbon total and effective sample sizes used to produce Figures 3(d-g), 4, and the TIMI analyses in Figure 5.

If the trends in archaeological proxy populations are substantially influenced by bottom-up constraints on food availability, the spatial asymmetry between the Alaskan and Kuril data series then implies inverse spatial structuring of the ecological conditions between the east vs. west sides of the North Pacific (Fig. 3; Fitzhugh et al., Reference Fitzhugh, Brown, Misarti, Krupnik and Crowell2020). At least for the Gulf of Alaska and Aleutians, these patterns also imply that the mechanisms forcing the marine ecosystems can be synchronous across large regions. Based on contemporary understandings of North Pacific climate variability, the intensity and position of the Aleutian Low (Rodionov et al., Reference Rodionov, Bond and Overland2007)—or some similarly manifested climate pattern—underly this oscillating pattern. In other words, local and regional near-shore ecosystems may respond to the spatially structured forcing of coupled atmospheric-oceanographic processes that affect the entire North Pacific Rim systematically, but not uniformly, at intervals of centuries to millennia. In Fitzhugh et al. (Reference Fitzhugh, Brown, Misarti, Krupnik and Crowell2020), we called this the Long Interval Oscillations in the North Pacific (the “LIONPac model”).

Evaluating the Possibility of Century to Millennial Scale Oscillations in Marine Ecodynamics and Paleodemography

The LIONPac model proposes the existence of multi-century to millennial-scale fluctuations in human populations caused by historically undocumented regime shifts in marine ecosystem. While unprecedented, researchers have only recently started to develop the techniques and accumulate the data to assess variability over such scales. Ultimately, evaluating this proposal will require coordinated interdisciplinary research to compile temporally resolved paleoecological proxy records (paleoclimate, paleoceanography, paleoecology, zooarchaeology, isotopic biochemistry, paleogenetics and archaeology) from around the North Pacific. Those data are not yet available at necessary spatial grain and temporal resolution, though we include preliminary comparisons of available evidence in what follows. A more immediately achievable approach is to validate the archaeological paleodemographic models over larger geographic regions to identify the scale(s) of coherence in the trends between regions.

Collaborations, initiated through the Paleoecology of Subarctic and Arctic Seas (PESAS) working group that inspired this special issue, have started assembling relevant evidence (see West et al., Reference West, Etnier, Barbeaux, Partlow and Orlov2020 and Nagashima et al., Reference Nagashima, Addison, Irino, Omori, Yoshimura and Harada2021, in this issue). This coordinated research program promises to open new understanding of the paleoecology of the North Pacific and to extend the dynamic baselines available to contemporary ecologists, managers, and planners. Much more work like this is needed before we can evaluate the ecological predictions of the LIONPac model or replace it with any alternative.

The data series that inspired the initial formulation of the model were limited to just three areas from the Northeast Pacific and one from the Northwest Pacific (Fig. 3d). While the Northeast patterns were suggestive, it is easy to imagine how the correlations could be incidental, and the lack of any adjacent data for the Kurils made it impossible to assess the robustness of the Northwest pattern. Here we present a significantly larger number of regional time series, assumed to represent relative population variation, and we evaluate a set of hypotheses derived from the LIONPac model. We argue that spatial coherence in proxy population time series can be used to strengthen interpretations about large-scale demographic processes and their potential causes (e.g., Crema et al., Reference Crema, Habu, Kobayashi and Madella2016; Tremayne and Winterhalder, Reference Tremayne and Winterhalder2017; Tremayne and Brown, Reference Tremayne and Brown2017).

Here we outline the hypotheses that frame our examination of the population models, and in the next section we examine the assumptions underlying the use of radiocarbon-based archaeological TFDs (temporal frequency distributions) as proxy population models. Where TFDs can be taken primarily as representations of population change through time, we hypothesize the following:

C1. Large-scale environmental or cultural changes that effect the productivity of ecosystems, economies, human health and wellbeing across adjacent regions should lead to positively correlated human populations trends through feedbacks on fertility, mortality, in-migration, and out-migration. These changes could include widespread ecological regime shifts, adoption of common technological innovations, or the spread of epidemic diseases or occurrence of catastrophic events (e.g., earthquakes, volcanic eruptions, tsunamis). Net in-migration from distant regions may not be linked to local causes, but the larger the scale of movement the more likely that adjacent regions will share the resulting growth trend. Such migration (in or out) might be related to shifts in the spatial distribution of productive ecosystems, the development/adoption of technologies that open new ecological niches to resident or immigrant communities, or access to or victimization by new methods of warfare.

C2. More local-scale environmental and cultural changes will more often result in inversely correlated population trends, for example when people leave one region in favor of an adjacent one. The scale of analysis matters in interpreting these patterns.

C3. Population movements that do not cross the boundaries of geographical units of analysis will not alter net population trends. Until recently, communities living primarily by hunting and collecting wild foods moved logistically if not residentially, and with some frequency—with longer moves less common than shorter ones. Different scales of spatial analyses then should capture different dynamics. Smaller spatial units should capture more frequent and less exceptional patterns of population movement, and population fluctuations in adjacent areas are more likely to be related. Larger spatial units, by contrast, should predominantly reflect endemic population variability, with transborder migration impacting trends only under the more extreme cases of large-scale migration events.

Unfortunately, the larger the scale of analysis, the more internal variability will be averaged in net population trends (see Timpson et al., Reference Timpson, Colledge, Crema, Edinborough, Kerig, Manning, Thomas and Shennan2014). This averaging reduces analytical sensitivity in spatial and temporal dimensions, while preserving dominant trends, but also sometimes accentuating unrelated but synchronous internal dynamics that may not be individually significant (i.e., causally unrelated but correlated patterns in different parts of the analytical region).

Expectations: C1 should be observable in positively correlated trends (up or down) in demographic (TFD) evidence. C2 should be detected in inversely correlated TFD trends in adjacent series. C3 simply acknowledges that there is a spatial scale dependency to population histories that will determine the sensitivity of the analysis and potential for interpreting observed patterns. Neutral growth in a portion of a time series does not imply that demographically meaningful events did not affect populations within the regions, only that they are not detectable at the resolution available.

Regions where people made their living primarily from the same ecosystem resources should have more similar population trends than regions of different ecotypes. For the purposes of these analyses, we define ecotype into marine (coastal) and terrestrial with the terrestrial divided into forested (boreal) and tundra zones. Boreal zones are found in the interiors of both Hokkaido and Alaska, while inland tundra zones are only found in regions of Alaska. Because the tundra zones also happen to be relatively close to the coasts, we expect them to potentially track coastal demographic features, so they constitute a class of ‘semi-coastal’ units when testing for coastal vs. interior relationships.

Expectations: Proxy population trends should track more closely between regions sharing the same ecological classification (coastal, tundra, or interior) than between them.

H3a. If the forces driving these patterns are predominantly marine (e.g., related to regime shifts in anomalously cold and warm water temperatures linked to the spatial shifts in the Aleutian Low system, as proposed by the LIONPac model), then coastal population trends on the same side of the Pacific should exhibit strong positive correlation, terrestrial and coastal population trends within macro-regions should be uncorrelated, and coastal population trends should be inversely correlated on either side of the North Pacific.

H3b. If the forces driving these patterns are predominantly terrestrial (e.g., temperature-linked impacts on caribou and deer survival), then non-coastal populations on the same side of the Pacific should exhibit strong positive correlations, terrestrial and coastal populations within macro-regions should be uncorrelated, and terrestrial population trends should be inversely correlated with those across the ocean. Based on synthetic paleoecological reconstructions for Alaska's interior regions showing highly asynchronous ecological variability over small distances (Kaufman et al., Reference Kaufman, Axford, Henderson, McKay, Oswald, Saenger, Anderson, Bailey, Clegg, Gajewski and Hu2016), we do not expect this hypothesis to be supported.

H3c. If the forces driving macro-regional ecological productivity cross marine and terrestrial ecozones (e.g., due to the integrating ecological effects and subsistence importance of anadromous fish or the greater integration of coastal and interior human populations), then population macro-regional patterns should be positively correlated within macro-regions and inversely correlated between them, irrespective of coastal vs. terrestrial classification.

Entanglement in expanding “world system” economies (Chase-Dunn and Grimes, 2002; Hudson, Reference Hudson2004) could have had positive and negative effects on population dynamics. Pursuit of resources for exchange and of trade partners with whom to trade could have stimulated migrations, social competition, and wealth (for successful actors) that might support population expansion. At the same time, aggressive competition could have led to excess mortality, balancing out gains in endemic population expansion, and expanding trade networks would have also supported transmission of pathogens with potentially significant impacts on previously more insular populations (Fitzhugh et al., Reference Fitzhugh, Phillips, Gjesfjeld, Whallon, Lovis and Hitchcock2011). Epidemic diseases are known to recur periodically in vulnerable populations until endemicity is established, which occurs more rapidly in dense and highly connected (e.g., urbanized) populations (Suzuki, Reference Suzuki2011). Smallpox was one of the first major diseases of colonial expansion into North Pacific communities. The disease hit central Japan by the 8^th century CE (Behbehani, Reference Behbehani1983) and could have reached Hokkaido through trade interactions soon after. Hokkaido and Kuril populations could have been impacted by multiple epidemics over the centuries of the low-level contact that prevailed from the late first millennium to mid second millennium CE (Hudson, Reference Hudson1999, 206–232). Russian-American conquest brought diseases first to the Gulf of Alaska region in the early 19^th century and later up the coast into central and northwest Alaska in the mid to late 19^th century (Fortuine, Reference Fortuine1989). Demographic impacts should then have been seen first on the south coast, then north coast, and lastly in the terrestrial interior.

In the Alaskan case, significant population growth is not anticipated (the numbers of Russian-American period colonizers were always very small relative to the Native populations), and the primary signal should be decline, based on historical documentation of the process (e.g., Lührmann, Reference Lührmann2008). Recent research has significantly undermined prior claims that virgin soil epidemics spread faster than the arrival of the colonial disease-bearers themselves (see Liebmann et al., Reference Liebmann, Farella, Roos, Stack, Martini and Swetnam2016; Feinman and Drake, Reference Feinman and Drake2021). If correct, then it is unlikely that synchronized population collapse in Alaska prior to Russian contact could be attributed to indirect transmission of diseases brought by European colonial incursions into eastern and southern North America and Mexico in the 16^th and 17^th centuries.

Visual interpretation of Proxy Population Curves

Figure 4 shows the archaeological radiocarbon TFDs for the last 10,000 years for 14 regional units across the North Pacific, including five from Hokkaido and one from the Kuril Islands (Fig. 4b) and eight from Alaska (Fig 4c) (see Figure 1 for geographical boundaries). The regions are defined as spatially contiguous zones, divided into coastal, near-coastal tundra (for Alaska), and interior boreal forest eco-regions to distinguish between regions supporting maritime vs. terrestrially dominated subsistence opportunities. Other than the Kuril Islands, we do not attempt to present data from the Russian Far East or Northeast Siberia because research histories and the resulting radiocarbon datasets for these regions are too small and fragmentary to produce meaningful TFDs at suitable spatial scales (Kuzmin, Reference Kuzmin2010).

Figure 4. Stacked, taphonomically ‘corrected’ TFD graphs of archaeological radiocarbon dates from regions around Alaska, Hokkaido and the Kuril Islands presented as proxy paleodemography models (see text for explanation and interpretation). Dark central trend lines are smoothed central tendency representations of the underlying data using the CRFPAB method (see Supplement). 95% confidence interval envelopes are shown as lighter grey upper and lower bounds. The left panel (4a) includes all regional data series; the central (4b) and right (4c) panels present greater detail of the Hokkaido+Kuril and Alaska series, respectively. Yellow and grey bars are manually placed highlights to show areas of notable trend increase and decrease, respectively. These highlights are placed wherever either the ending min-max envelope has moved out of the range of the starting min-max envelope (making a neutral growth interpretation unlikely) or where the central trend line at least doubles or halves during the excursion. Blue dashed lines mark 5000 yr, which is the approximate cut-off for visual analysis and absolute cut off for the TIMI spatial autocorrelation (earlier trends are insufficiently robust for comparison). The purple dashed arrow in panel 4b illustrates the time transgressive shift from North Hokkaido to the Kurils thought to track expansion of the Okhotsk Culture. The red dotted arrows on Figure 4b and 4c demarcate the largely synchronous increase in populations, starting ca. 1300 yr in Hokkaido and after ca. 750 yr in all Alaska. The narrow black dashed lines in the central (4b) and right (4c) panels depict the correlated population collapses indicated for Hokkaido/Kurils ca. 800-500 yr and in Alaska ca. 400 yr and continuing through colonial contact. Blue block arrows on the top right corners of 4b and 4c point to major late Holocene TFD growth and collapse trends across multiple time series for the Hokkaido/Kuril and Alaska data sets, respectively. These trends are negatively correlated in time between the two macro-regions (see text for discussion). The gray strip across the top of 4b delineates Hokkaido culture historical phases: IJ = Initial Jomon, EJ = Early Jomon, MJ = Middle Jomon, LJ = Late Jomon, FJ = Final Jomon, E-J = Epi-Jomon, O = Okhotsk, S = Satsumon, A- Ainu. Culture historical designations of relevance are shown on the curves. See Figure 2 for a more complete set of culture historical schemes.

For visual inspection, the 14 series in Figure 4 are arrayed in approximate geographical proximity and stacked as a whole set (Fig. 4a), as a Northwest Pacific stack (Hokkaido + Kuril Islands; Figure 4b), and as a Northeast Pacific stack (“Alaska Series,” from the Arctic coast to Kodiak; Figure 4c). Each trend plot includes a 95% confidence envelope, generated by Monte Carlo simulations, that better represents the changing effect of sample density through time in each series (see Supplemental Materials).

Visual inspection of Figure 4 allows comparison of trends between adjacent series. Notably, each series starts at a different time depth due to the uneven antiquity of archaeological evidence in each region (here, all series are truncated at 10 ka), and while some differences in the early part of the time series may be archaeologically meaningful, they render many TFD comparisons problematic in the early Holocene (and for some comparisons, prior to 2500 yr). As a result, we begin the visual analyses in the middle Holocene when archaeological data becomes available from the Bering and Arctic coasts of Alaska (and Siberia) (Fitzhugh, Reference Fitzhugh, Friesen and Mason2016; Tremayne and Brown, Reference Tremayne and Brown2017).

Before considering any subset of curves, we note that all curves can be decomposed into at least two scales of periodicity, which we will call long-period (>1000 yr) and short-period (<1000 yr) variability). Long-period variability is generally characterized by larger fluxes and enduring modes or regimes, while the transitions between them can be abrupt, often within 100–200 years. Short period variability is manifest as fine-scale fluctuations that reverse within a few hundred years. These short-period “sawtooth” fluxes have relatively low amplitudes and, for the most part, vary without transgressing the 95% confidence interval windows. This short-period variability cannot be distinguished from a null model of neutral growth and so is potentially meaningless in depicting population change. Even shorter transitions may have occurred that are obscured by the kernel density smoothing function used. Notable short- and long-period positive and negative trends are depicted in Figures 4b and 4c with colored shading (see Fig. 4 caption for explanation).

TIMI Method

For each time step, the TIMI procedure calculates the Moran I value (Anselin Reference Anselin1995) for each pairwise comparison of regional TFDs, where certain filtering criteria are applied based on the hypothesis under examination (in our case: adjacency, ecological similarity, or macro-regional unity overall and for coastal data). We compute running series of global Moran I values to evaluate each of the hypothesis. This is done by filtering pairwise relationships by means of a “1” (“include”) when the relationship is one of membership in the hypothesized in-group (adjacent/sharing a border; sharing a coastal or interior ecotype; sharing a macro-region/Alaska vs. Hokkaido+Kurils) and “0” (“exclude”) when not. For example, in testing for adjacency, we set TIMI to run paired correlation analyses only on those regions that share a border (coded as “1”) relative to the average correlation between all regions. The TIMI analyses allow us to evaluate changes in the autocorrelative relationships between different time series as they change through time. This makes it possible to evaluate the initial research questions/hypotheses temporally, to reveal previously unconsidered patterning in the data, and to further address questions about sample sufficiency and the evolution of human-environmental dynamics across the North Pacific.

Running Moran I plots are calculated on the mean adjusted growth rates (MAGR or $\bar{{\rm r}}$), which is a transformation of the TFD plots reflecting the changing slopes of the TFDs in Fig. 4 (the first derivative), not the ‘raw’ TFD values themselves. The Moran I statistic estimates, for each time step, the degree to which pairs of MAGR plots are doing the same thing (increasing, decreasing, staying constant) averaged across the within-group set (i.e., the set coded for either adjacency, identical ecosystem class, or membership in the same macro-region) relative to averaged pairwise comparisons of all regions. On the TIMI plots presented in Figure 5 (a, c, e, f), any index value above the mean expected value (E[l_t], solid red lines) indicates a greater than average positive autocorrelation in favor of the hypothesis under evaluation. Confidence intervals (CIs) of 95% (α=0.05) are calculated based on repeated permutations of the MAGR values randomly reassigned to dummy series (Supplement). With the confidence interval set at 95% and a relatively small number of regions included in the analyses, our TIMI results rarely exceed the “fail-to-reject” threshold. Any random error noise in the underlying TFDs would also tend to dampen TIMI values below the confidence interval levels, making it especially hard to achieve significant measures of autocorrelation. Given the typically smaller sample sizes earlier in each time series, we might expect less significant autocorrelations earlier in time. If so, then autocorrelation closer to the present may be an indication that the hypothesis is correct. Acknowledging those issues, we interpret the results with attention to both significance and the intervals spent substantially above or below the mean expected value (E[l_t]).

Figure 5. (a, c, e, f) Time Iterative Moran Index (TIMI) plots of global Moran I statistic (I) calculated on pairwise MAGR datasets to measure degrees of spatial autocorrelation as they change through time relative to hypothesized relationships. (b and d) Trends in Mean Adjusted Growth Rate ($\bar{r}$) plotted for each regional time series and color coded (in graph b) for coast (blue) vs. interior (red) and (in graph d) for macro-regional association (blue for Alaska; red for Hokkaido + Kurils). (a) Moran's I adjacency analysis (“Do adjacent regions tend to be more similar than non-adjacent regions?” Conventional spatial autocorrelation); (c) general eco-region analyses (“Are coastal vs. Interior data sets more autocorrelated within their sets than between them?”); (e) Shared macro-region analysis for all ecoregions, with stars denoting peak positive autocorrelations in the last 1100 years (“Do regions within Alaska and regions within Hokkaido autocorrelated with other regions in the same macro-region?”); (f) Shared macro-region analysis for coastal and near-coastal data sets (“Do TFD series for coastal and near coastal regions within the same macro-region autocorrelated relative to other regions and all regions in the other microregion?”). Gray overlays on panels a, c, e and f denote Moran I values that exceed the formal 95% “fail to reject” threshold (α=0.05) in a positive or negative direction. See text and Supplement for discussion of the TIMI approach to spatial autocorrelation of time series and the calculations applied.

Figure 5b and 5d show the mean adjusted growth rates ($\bar{r}$) for each time series. The trends are noisy, which is expected based on our analyses of the stacked TFD series in Figure 4. Figure 5b is color-coded to identify shared eco-regions (Coastal zones = blue; Interior regions = red). Figure 5d shows the same plots with each trend color-coded to indicate Alaska (blue) and Hokkaido + Kurils (red). These MAGR plots are included to allow visual inspection of the data underlying the TIMI analyses.

Next, we discuss these results in the context of the original hypotheses.

TIMI Results

Our application of the TIMI calculation for this hypothesis outputs above-expected values when growth rate trends in adjacent regions correlate more than with non-adjacent regions (hypothesis confirmation). As a result, a positive TIMI tendency reveals a trend in support of the first corollary (C1) of Hypothesis 1, that of shared growth or decline due to net fertility relative to mortality and/or in-migration relative to out-migration from non-adjacent regions. Neutral or negative pulling values would indicate more uncorrelated or inversely correlated trends in neighboring series, including net population movement from one adjacent region to another (C2).

In Figure 5a, the Moran I statistics cross the confidence threshold in the positive direction (indicating similar trends) four times, briefly around 2700 yr, 2000 yr, 1100 to 1000 yr, and again for a longer interval from ca. 800 to 450 yr. Given the sampling issues already discussed and the likelihood that the underlying TFD patterns are formed by multiple causes operating at different spatial scales, it is not surprising that the TIMI analyses rarely reach the 95% significance threshold. As a result, we are comfortable attending to persistent positive and negative excursions from the expected value (E[l_t]). Positive TIMI autocorrelations of adjacent regions are notable for most of the last 3000 years. Earlier patterns are variable or even predominantly negatively autocorrelated. In marked contrast to intervals before and after, the TIMI plot trends negative from 1900 to 1400 yr. This was an interval of considerable population flux from one region to another in Hokkaido. It was also a time of relative constancy but low level and uncorrelated flux in the Alaskan series (Fig. 4).

In broad terms, we interpret the results of the adjacency analysis to support the hypothesis that neighboring regions are more likely than not to share similar population trends, particularly in the most recent millennia, when data sets are most robust. Differences in regional scales between Alaska and Hokkaido likely complicate this analysis, something we return to in the Discussion.

This hypothesis gets at the idea that hunter-fisher-gatherer populations are, to a significant degree, dependent on the resources they can harvest within their regions. It assumes that those resources are more likely to respond in parallel to similar environmental changes, including changes that alter the net availability of food (holding harvesting methods constant), and that the productivity of food harvests is driven by bottom-up ecological factors. These assumptions underly both Hypotheses 2 and 3, but we approach the issue first without regard to differentiation in how our ecozones in different locations might experience environmental drivers differently. If human populations associated with North Pacific Rim ecosystems followed the same trends from Hokkaido to Alaska members in pan-subarctic/Arctic ecosystem, this hypothesis should be supported. Accordingly, our Eco-Similarity analysis lumps all coastal and near-coastal regions together and all Interior regions together and computes the TIMI autocorrelation within each group compared to the total set of regions. Figure 5c shows the results of the Eco-Similarity test. In this case, there is little if any autocorrelation, or even patterning, in the data—a characteristic seen also in the color-coded MAGR plot in Figure 5b. We can safely reject this hypothesis.

Hypothesis 3 embraces the spatially structured expectation of the LIONPac model. Figure 5d shows the same MAGR plot as in Figure 5c but colored for macro-region affiliation (blue for Alaska, red for Hokkaido). Despite the noisiness, there are clear correspondences between subsets of time series through various intervals. From the colors, it is particularly apparent that Alaskan series tend to track together predominantly from 2500 yr to the present, in distinction from the Hokkaido+Kuril series, which is both more chaotic and somewhat inverted relative to the Alaskan series over the last 800 years.

Figures 5e and 5f each evaluates a variation of Hypothesis 3. Figure 5e measures the extent of autocorrelation between regions within one or the other of Alaska and Hokkaido+Kurils macro-regions, irrespective of eco-zone (i.e., testing H3c). High positive values, such as the three dominant peaks between 1400 and 150 yr (starred), indicate positive correspondence within those macro-regions. Other intervals of high if not significant (at α=0.05) positive correlation are seen from 4000 to 3400 yr and, to a somewhat lesser extent, 2800 to 1900 yr. Because these analyses run both ways (looking for greater than average correspondence within both macro-regions), the greater discordance in the Northwest Pacific data sets are likely responsible for the low autocorrelation elsewhere in the series. It is notable that the interval from 800-450 yr is the most significant interval in all of the TIMI tests (Figs. 5a, 5e and 5f) except for Figure 5c (Eco-Similarity) and it is the most significant of any of the TIMI analyses in this test.

Figure 5f is an attempt to directly test the LIONPac model by looking at the coastal set of samples by macro-region (Hypothesis H3a). The result shows more accentuated positive and negative autocorrelation. While significant positive autocorrelation is indicated again for the interval of 650-450 yr in Figure 5e, it is less significant than for the total macro-region analysis (Fig. 5f). Based on these two tests, we would have to reject the narrowly defined LIONPac model. Our effort to quantify the degree of autocorrelations using TIMI forces us to conclude that the archaeological population proxies are not predominantly patterned as expected if only populations in coastal regions were affected by an east-west oscillating (marine) ecosystem dynamic. This does not mean that the LIONPac model is necessarily wrong, or that there are not some other forces driving regional population variability. It is highly likely that the rejection of the marine specified hypothesis (H3a) is in part a function of the greater heterogeneity of the Hokkaido TFDs compared to those of Alaska. We explore the factors that could be behind this below and in the Discussion.

This hypothesis formalizes the proposition that the inverse trends initially identified between southern Alaska and the Kuril Islands (Fig. 3d; Fitzhugh et al. Reference Fitzhugh, Brown, Misarti, Krupnik and Crowell2020) are, in fact, artifacts of the expansion of commodities trade networks, epidemic diseases, and ultimately colonial incursions, first in maritime Northeast Asia and later in Alaska. The increase in autocorrelation seen in Figures 5e and 5f after 1400 yr, but especially after 800 yr, and the more strongly aligned positive growth rates seen in the Hokkaido plots (red in Fig. 5d) from 900-500 yr might be taken as evidence for population growth connected to beneficial consequences of access to trade with Japan, Sakhalin, and the mainland by Satsumon, Okhotsk, and ultimately Ainu communities. The expansion the Okhotsk may also have been driven by the development of more productive sea mammal hunting techniques, potentially opening a new niche for Okhotsk populations to thrive and grow.

This explanation does less to help us make sense of the autocorrelated Alaskan series through time. Colonial impacts arrived in Alaska between 200 and 100 yr (1750-1850 CE), at the very tail end of radiocarbon discrimination. It is well-documented that Alaska Native populations were heavily impacted by introduced diseases and associated famines in waves during the 19^th and early 20^th centuries (Fortuine, Reference Fortuine1989), while Ainu communities were ravaged by conquest, resettlement and coerced labor starting in the 17^th century, but especially in the late 19^th and early 20^th centuries (Walker, Reference Walker2001). These colonial impacts doubtless contribute to the declining population trends on both sides of the North Pacific in the last couple of centuries, but they do not explain the earlier onset of these trends, especially in Alaska starting ca. 400 years ago.

One limitation in evaluating the colonial impacts hypothesis (H4) from the macroregional ecology hypotheses (H3) is the possibility that the TIMI analyses are sensitive to decreasing sample robustness and consequently more noise in TFD signals as we move back in time. While further investigation of the TIMI method will be needed to clarify how sample robustness relates to TIMI outputs, for now we recognize that the presence of strong autocorrelation in the macroregional data sets primarily in the last 1100 years (Figs. 5e and 5f), could support both Hypotheses 3 and 4, and both could be in operation together. For the time being, we can only cleanly reject Hypothesis 2. There is no support for the idea that coastal populations around the North Pacific Rim were locked into similar trends of growth and decline in distinction to terrestrial populations on both sides of the Pacific/Bering Sea.

What does it all mean?

Quaternary paleoecologists are well aware that past environmental changes have often exceeded the bounds of recent variability, and world-class climate modelers seeking to forecast future impacts under alternative scenarios of anthropogenic inputs regularly include paleoclimate proxy data to evaluate the performance of model assumptions. Even so, it is often challenging to compile robust paleo-records to track the dynamics of past climates and environments. Given the vital importance of fisheries to contemporary fishing communities and to the global food supply, it is as important to gain better understanding of the histories of marine ecosystems as it is to reconstruct past climates (which of course are related problems). This effort is still in its early stages. Arguably, it is even more important to understand how people have interacted with marine ecosystems and how those systems have in turn affected human ways of life and societies. People have been living around and by means of the greater North Pacific Rim ecosystem for 5000 to 10,000 years or more, and yet we have very limited understanding of how those communities managed their lives and economies in the face of significant climate and ecosystem shifts.

In this paper, we examined the proposition that a North Pacific Rim long-interval oscillating ecosystems dynamic has been operating at millennial scales through the Holocene. Lacking the data to consider the proposition directly, we considered it through the lens of archaeological paleodemography proxies. Preliminary analysis (Fitzhugh et al., Reference Fitzhugh, Brown, Misarti, Krupnik and Crowell2020) had shown that a selection of four proxy population reconstructions from southern Alaska and the insular Northwest Pacific (Kuril Islands) fluctuated as if facing century to millennial scale oscillating cycles of food availability and insecurity. We call the hypothetical socio-ecological model the long-interval oscillation in the North Pacific or LIONPAC model. The model envisions neighboring populations experiencing synchronous shortages in some intervals and surpluses in others, in opposite phase on either side of the North Pacific. One plausible mechanism to explain these patterns is the presence of a multi-century to millennial scale ecological regime dynamic driven by persistent intervals of strong and low Aleutian Low dynamics.

Building on the original analysis, here we have brought together fourteen archaeological demographic proxy models from across Alaska, Hokkaido, and the Kuril Islands to test the integrity of the original interpretations and to facilitate broader analyses. To our knowledge this represents a novel exploration of a large set of archaeological population proxies, taking advantage of the greater interpretability of neighboring distributions with partially shared archaeological, cultural, demographic, and environmental histories. We conducted qualitative visual comparisons and introduced the Time Iterative Moran I (TIMI) spatial autocorrelation method to compare archaeological temporal frequency distribution (TFD) trends quantitatively. Our population models are all built from cleaned archaeological radiocarbon data sets assembled into occupation-index-based temporal frequency distributions or TFDs.

Results indicate considerable population dynamism around the North Pacific. Adjacent and neighboring trends show strong associations in time and trend that reinforce the assumption that (for the most part) these proxy models represent meaningful variation in the intensity of archaeological deposition, usually as a reflection of changing population densities, reflecting net fertility, mortality, in-migration, and out-migration. The combined analyses culminated in the test of several hypotheses constructed around expectations for demographic relationships between neighboring regions (adjacency), coastal vs. interior locations (as a proxy for targeted ecological adaptation and susceptibility to ecologically differentiated fluctuations in resource availability), and membership in one or another microregion (Northwest or Northeast Pacific). Evaluating the predictions of these hypotheses, we show that adjacent regions were most likely to share demographic patterns and that membership in one or the other macro-regions contributed strongly to the likelihood of shared demographic trends, especially for Alaska. Interestingly but not surprising in hindsight, the ecosystem similarity test failed (Fig. 5b,5c), revealing that ecosystem dynamics, and hence human food security, is not strictly determined by whether people lived on the coast or interior.

Differences in spatial autocorrelation of Alaska vs. Hokkaido proxies suggests either that Hokkaido populations were far more spatially dynamic than those of Alaska, or that the different sized spatial units characterizing these two regions, respectively, open windows into different scales of demographic behavior. As such smaller regions, such as those within Hokkaido and the Kurils, may be characterized by more inter-regional movement and population shifts; whereas behavior at these scales operates ‘below the radar’ at the larger scales reflected in the Alaskan regions, which are thereby more reflective of more macro-scale (and interestingly more correlated) trends. Even so, Alaskan archaeological history is characterized by several major migrations, and these can also be detected in the population proxy modes across the Alaska data sets.

While the LIONPac model remains hypothetical at this point, the aggregate (macro-regional) comparison (Fig. 6d) suggests that even with greater heterogeneity in details, the Northwest Pacific data (Hokkaido+Kurils) encompass significant net growth and decline in generally opposing phase to that of the Alaskan series. This comes closest to supporting the hypothesis that populations were under the influence of some coherent millennial scale dynamic. The more regularly oscillating Alaskan pattern, in particular, comes closest to supporting something like the LIONPac model. This is especially true when viewed in consort with available paleoclimate and paleoecological proxy series that carry the signal of millennial scale fluctuations in Aleutian Low strength and position (Fig. 6).

Where do we go from here?

The paleodemographic proxy data presented in this paper contain a wealth of information that should be scrutinized, challenged, tested and interpreted by researchers interested in understanding the dynamics of past populations around the North Pacific. We should seek to fill the gap in distribution around the Chukchi and Bering Sea by building a more comprehensive set of radiocarbon data from sites around eastern Siberia and the Sea of Okhotsk beyond Japan.

The TIMI autocorrelation analysis is experimental and needs refinement to be most useful for comparison of contemporaneous trends and time series. We think the method holds potential for quantitative analysis in comparative time series analyses where visual evaluation is challenging. This could prove important in Quaternary paleoecology, paleoclimatology, archaeology or any similar field that produces and compares stacked time series to identify correlated features and infer causal relationships.

Finally, we believe that interdisciplinary perspectives such as those presented here are essential for improving our understanding of critical socio-ecological dynamics that have fallen out of knowledge but are especially critical in managing intensively exploited ecosystems under conditions of growing climate change and futures that are less predictable than the recent past. Interpolating this ancient knowledge from archaeological assemblages and environmental proxies and linking it to our best available understanding of human-environmental processes around the North Pacific is needed for anticipating potential future challenges and supporting the ongoing resilience of resident human communities and the ecosystems they depend on.