Island Colonization and Environmental Sustainability in the Postglacial Mediterranean

: Island environments present challenges to human colonization, but we have a poor understanding of how environmental difference drives heterogeneous patterns of insular settlement. In this paper, we assess which environmental and geographic variables positively or negatively affect the long-term sustainability of human settlement on islands. Using the postglacial Mediterranean basin as a case study, we assess the impact of area, isolation index, species richness, and net primary productivity (NPP) on patterns of island occupation for both hunter-gatherer and agropastoral populations. We ﬁnd that models involving area most effectively accounts for sustainability in hunter-gatherer island settlement. The agropastoral data are noisier, perhaps due to culturally speciﬁc factors responsible for the distribution of the data; nonetheless, we show that area and NPP exert profound inﬂuence over sustainability of agropastoral island settlement. We conclude by suggesting that this relates to the capacity of these variables to impact demographic robusticity directly.


Sustainability and Island Settlement
Islands present difficulties for human settlement-not only in terms of the physical challenges that reaching and colonizing them pose, but also because of biotic and physiographic processes peculiar to insular landmasses [1][2][3][4]. Nonetheless, humans have successfully, if differentially, established permanent settlement on many islands around the planet over the last 50,000 years, and possibly earlier [5][6][7][8]. Using the Mediterranean basin as a case study, we seek to understand which, if any, environmental factors have rendered human settlement on islands sustainable over millennial scales. In particular, we try to assess how the transition from hunter-gatherer (food-acquiring) to agropastoral (food-producing) subsistence behavior altered the parameters for long-term sustainable settlement on islands (i.e., any landmasses not connected to the mainland at the time of their earliest occupation).
The focus of archaeology on human-environment dynamics over the long term means that it is ideally situated to address questions of this type and scale. To do so, a relatively rare set of co-occurring circumstances are necessary. First, it is only possible to explore how environmental differences drive variance in the sustainability of insular settlement when such settlement is supported by broadly similar subsistence regimes distributed over broadly dissimilar island environments. Second, those environmental variables under investigation must be quantifiable in a consistent way [9]. Third, addressing the question of long-term island sustainability requires settlement histories with differing time depth, such that we can statistically evaluate the role of environmental variables in affecting the antiquity and persistence of island settlement. The postglacial Mediterranean basin, with its long record of settlement and its many and diverse insular landmasses, provides Sustainability 2021, 13, 3383 2 of 20 an excellent context in which to assess the roles of certain environmental variables in promoting or impeding sustainable settlement on islands.
Given the demonstrated ability of humans to reach even remote Mediterranean islands since at least the late Pleistocene, we take differences in the sustained settlement of those islands as reflecting their respective environmental attractions and limitations. We amass data on patterns of island settlement by both hunter-gatherer and agropastoral populations. We then assess the statistical relationship between colonization date for all islands with a surface area >10 km 2 and a series of variables likely to affect sustainability of settlement, including area, isolation index, net primary productivity and species richness, and (latterly) longitude. We begin by critically examining the complex relationship between the very long-term perspective adopted by archaeology and the constituent components of the sustainability concept.

Conceptualizing Sustainability over Archaeological Timescales
As the study of human behavior and its mutability over time, archaeology poses conceptual challenges for the concept of sustainability. As Fisher shows, this is because of the historically contingent [10] (p. 394) aspect of sustainable modes of subsistence. Human behavior and the biophysical feedbacks with which it intersects are intrinsically dynamic over centennial to millennial scales. Even if modes of subsistence develop that do not deplete or reorder the ecosystems in which they are situated, both these modes and their ecosystemic context represent moving targets. Thus, a reasonably sustainable mode of subsistence over meso-scales (for example, highly mobile hunting oriented towards migratory ungulates during glacial periods in western Eurasia) can become inherently unsustainable (for example, during interglacials as expanding boreal forest and associated fauna succeed tundra). When viewed over the very long term, it seems likely that no mode of subsistence is truly sustainable. Comparably, as Fisher [10] also stresses, we should not make the culturally uniformitarian mistake of assuming that humans recognize or value what we now describe as sustainable subsistence practices. Current human society is unique in impacting the majority of Earth's biosystems [11,12], but there are examples in deeper time of human behavior driving unidirectional change in the biosphere and beyond, suggesting that preindustrial lifeways were not intrinsically sustainable, e.g., [13][14][15][16][17].
With that said, certain types of human economic and subsistence organization can clearly persist over very long periods of time. Mid-to low-latitude mixed hunting and gathering in Africa and Eurasia, for example, is demonstrably sustainable over multimillennial scales. Beyond sustainability, the robustness of such modes of subsistence hints at their resilience-that is, their ability to adapt to changing socio-ecological conditions. Marston [18] shows the utility of resilience theory for the sustainability concept, and especially how it refocuses inquiry on adaptive cycles and explicitly allows us to address cause and effect. Critically, resilience theory emphasizes the importance of scales of analysis. Thus, over large spatial and temporal scales, we can conceptualize modes of subsistence as either resilient or transient. Without wishing to parse this terminology more minutely, and recognizing Fisher's [10] arguments about the utility of the sustainability concept over very long spans of time, we view resilient and sustainable economic organization as essentially identical: longevity of settlement that, premised on modes of subsistence, can persist over millennial scales without causing environmental degradation to the extent that it becomes unviable. Here, we do not address certain other aspects of sustainability (e.g., [19,20]); we simply examine durable settlement over millennial scales founded on known types of subsistence.

Sustainability and the Mediterranean
Since the Last Glacial Maximum~20,000 years ago, two modes of subsistence have dominated in the Mediterranean basin. Until about 9000 years ago, hunting and gathering was the prevalent mode of economic organization, as in much of Africa and Eurasia. Initially focused on targeting high-ranked species moving over extensive ranges, by the terminal Pleistocene and into the early Holocene hunter-gatherer strategies had shifted towards more limited seasonal mobility involving the exploitation of more varied taxa [21,22]. This subsistence behavior supported a small population around the Mediterranean littoral. Subsequently, a form of agropastoralism, ultimately derived from Southwest Asia, emerged as the predominant mode of subsistence. The proliferation of this form of economic organization was rapid, essentially occurring in an east-to-west cline from approximately 9000 to 7000 BP so that by 7000 BP much of the Mediterranean littoral (as well as parts of their peninsular interiors) was subject to agropastoral exploitation [23,24]. Broadly premised on the mixed utilization of domesticated ungulates, cereals, pulses, and latterly fruits in heavily managed and anthropogenically impacted landscapes, this agroeconomic arrangement has proved durable over the long term. Notably, it has supported high population densities and early forms of state-level societies in a relatively liminal environment [25]. Moreover, the longevity of this agroeconomic mode is substantial enough that the physiographic organization in the Mediterranean reflects a high degree of anthropogenic disturbance (e.g., coastal/riverine progradation episodes reflecting land-use change and corresponding sediment transport), and the current biota reflect strong anthropogenically influenced selection (e.g., the prevalence of sclerophyllous vegetation under selective pressure from ungulate herbivory).
Which of these modes of subsistence was more sustainable in an insular context, and what made it so? Did both hunter-gatherers and agropastoralists express clear preference for certain types of island environments, one that has broader implications for human insular ecodynamics? The Mediterranean confers several advantages when it comes to addressing how environmental variability drives differential outcomes with regard to sustainability, and these derive from: (a) relatively clear spatial-temporal structure in hunter-gatherer and agropastoral island settlement throughout the basin; and (b) the basin's physiographic organization. Hunter-gatherer groups, and after 9000-7000 BP farming groups, established themselves on the continental littoral broadly, but the same is not true of the Mediterranean islands. Some (including Sicily, Sardinia, Corsica, Cyprus, and Crete) were clearly targets for both hunter-gatherer and agropastoral settlements, but other islands appear to have been ignored by early Holocene groups and during subsequent processes of Neolithization during the early and middle Holocene. In stark contrast to other western Mediterranean islands, for example, the Balearic Islands have no pre-farming archaeological record and the arrival of the farming lifestyle there postdates its arrival on the neighboring Iberian coast by three millennia [26,27].
As a result, some Mediterranean islands have deep and enduring histories of settlement; others have much shorter histories, sometimes marked by apparent periodic abandonment and resettlement [28] (pp. 180-208). Why is this? In other contexts, such as the insular Pacific, size, distance, and degree of isolation-geometric variables that drive biogeographic patterning more broadly-have been shown to structure insular settlement, e.g., [2,29]. Indeed, island area has, since Cherry's early work [30][31][32], consistently been viewed as a major factor in explaining variation in the data, both for agropastoral settlement [28] and, more recently, in preceding hunter-gatherer exploitation [33]. There is substantial evidence, however, that Mediterranean hunter-gatherers had the technological capacity to access relatively small and remote islands when the benefit existed to do so, e.g., [34]. Early farmers moving westward were evidently also accomplished sea-goers [23]. As a result, it is probably not behavioral or geographic constraint which drove this heterogenous temporal patterning in human settlement. Humans have been able to move deep into the insular Mediterranean since the last glaciation, but have either chosen not to do so, or to do so in a highly discontinuous manner. Why?
We assume that potential colonists preferred environments in which their mode of subsistence was more intrinsically viable and sustainable, and avoided those that were hostile or challenging from a subsistence perspective. Given the general absence of scalar constraints on dispersal limiting modern human populations, antiquity of island colonization is therefore an acceptable proxy for intrinsic environmental attractiveness. Thus, we propose that the spatio-temporal structure in the settlement of the Mediterranean islands by hunter-gatherers and farmers should reflect meaningful environmental differences that promote or inhibit long-term sustainability for each economic mode. Might we therefore predict a simple relationship between island size and capacity to promote hunter-gatherer and agropastoral sustainability? Possibly, in that gross area can permit raw population loads greater than smaller areas, other factors being equal; but as Cherry and Leppard [33] stress, the problem is more complicated. Area in and of itself cannot be assumed to be the variable accounting for patterning in the data, in that it can co-vary with other environmental variables which might reasonably be supposed to impinge on settlement dynamics. Area, for example, can correlate with elevation and thereby hydrological patterns (via orographic precipitation), and it can have non-trivial impacts on biotic richness, both in terms of basic biodiversity but also in taxon-specific demographic loads, all relevant factors in accounting for settlement histories. We need to evaluate the importance of these factors individually. Moving beyond simple geographic properties such as area, Mediterranean islands also differ in other dimensions (e.g., latitude, geological and pedological diversity, degree of remoteness, and primary terrestrial and marine productivity) which are less correlated or intrinsically less sensitive to one another but may also affect, albeit to a lesser degree, suitability for agropastoral or hunter-gatherer settlement.
We are interested in quantitatively evaluating the impact of environmental and geographic variability on sustainability in insular human populations over the long term. As a result, we need to separate these variables and assess their impact on the temporal structure of the data independently and quantitatively. As many of these characteristics are quantifiable, we should in principle be able statistically to test the relationship between these variables and settlement histories. We do this with a focus on those environmental variables which research suggests are most likely to explain patterning in the data efficiently: area, degree of isolation, and primary biotic information.

Methodological Assumptions
To test the assumptions reviewed above, we created statistical model sets of Mediterranean islands, their associated environmental and geographic attributes, and dates associated with their occupation histories (see Supplementary Material for RMarkdown file and associated files). This work was conducted using the R programming language and statistical software [35]. We utilized the global GADM (Database of Global Administrative Areas [36]) to acquire a complete and high-resolution dataset of Mediterranean islands, numbering 2019 in total. This dataset was processed and subset prior to analysis, by which procedure we selected only those islands 10 km 2 or greater in area. This threshold was a largely arbitrary choice, balancing the need for a sufficiently large dataset against the likelihood that appropriately high-resolution dates would be available and published. We also took into consideration the fact that the majority of islands smaller than 10 km 2 in area are so small that settlement of any kind has probably always been precluded. We recognize that larger islands under this size (i.e., those between 10 km 2 and 1 km 2 ) appear to experience highly stochastic and discontinuous settlement [37], and this threshold thereby skews our analysis towards islands that are less likely to experience discontinuous settlement; moreover, many islands that today fall below the 10 km 2 threshold would have been somewhat larger (via glacial melt-driven eustatic sea-level change) in the late Pleistocene and earlier Holocene. Yet, if sustainability of settlement is sensitive to area, this should emerge during the analysis, allowing us to make secure predictions about likely dynamics on islands smaller than 10 km 2 in area.
With this 10 km 2 threshold, we obtained a dataset of 161 islands, or 8% of the total 2019 islands acquired via GADM. Of these 161, 16 (9%) had published dates for huntergatherer occupation and 106 (66%) had published dates for earliest agropastoral occupation ( Figure 1). Future iterations of this research are free to apply different thresholds, but our island dataset covers much of the expanse of the Mediterranean and presents a represen-Sustainability 2021, 13, 3383 5 of 20 tative sample of Mediterranean island ecologies. For each of these islands, we calculated area, isolation index, mean terrestrial net primary productivity (NPP), and mean species richness values. We selected these variables based on the relative ease of their calculation given the size of our study area and their relevance for general environmental processes, including the organization of trophic systems and their likely capacity to affect human subsistence behavior. eustatic sea-level change) in the late Pleistocene and earlier Holocene. Yet, if sustainability of settlement is sensitive to area, this should emerge during the analysis, allowing us to make secure predictions about likely dynamics on islands smaller than 10 km 2 in area.
With this 10 km 2 threshold, we obtained a dataset of 161 islands, or 8% of the total 2019 islands acquired via GADM. Of these 161, 16 (9%) had published dates for huntergatherer occupation and 106 (66%) had published dates for earliest agropastoral occupation ( Figure 1). Future iterations of this research are free to apply different thresholds, but our island dataset covers much of the expanse of the Mediterranean and presents a representative sample of Mediterranean island ecologies. For each of these islands, we calculated area, isolation index, mean terrestrial net primary productivity (NPP), and mean species richness values. We selected these variables based on the relative ease of their calculation given the size of our study area and their relevance for general environmental processes, including the organization of trophic systems and their likely capacity to affect human subsistence behavior.

Area
Other factors being equal, island size directly relates to demographic and biotic dynamics in two primary ways related to human settlement. First, greater area can in principle support larger populations of a given taxon, humans included. For a number of reasons-including the high energetic needs and extended juvenility of our species, as well as the cross-cultural prevalence of incest avoidance-smaller populations generally

Area
Other factors being equal, island size directly relates to demographic and biotic dynamics in two primary ways related to human settlement. First, greater area can in principle support larger populations of a given taxon, humans included. For a number of reasons-including the high energetic needs and extended juvenility of our species, as well as the cross-cultural prevalence of incest avoidance-smaller populations generally exhibit a fragility which increases allometrically over larger populations. This implies the existence of scalar thresholds below which human populations on islands must necessarily maintain an exogamous relationship with an off-island metapopulation or risk exposure to exaggerated stochastic effects [38]. Second, island area can covary grossly with scale (and thereby complexity and ultimately resilience) of trophic systems, and with biodiversity via species-area relationships (SARs) [39][40][41][42]. Large trophic systems are intrinsically likely to promote sustainability, especially in food-acquiring subsistence economies; moreover, biodiversity may affect sustainability for reasons we discuss below.

Isolation
For the reasons outlined above, we suspect that "isolation" is not likely to be a primary predictor of patterns of colonization and we suppose colonization date to be a reasonable proxy for the intrinsic potential of an environment to promote sustainability. In other contexts, however, the physical geometry of the organization of landmasses clearly has patterning effects [2,3], and as a result we consider it here for its possible implications in terms of date of settlement and resulting degree of connectedness to mainland resources and populations. Isolation is difficult to define in quantitative terms, since straight-line distance from the nearest mainland or island can gloss relevant, intervening geographic organization (such as strings of stepping-stone islands), and the place from which a source population was derived can only rarely be known with certainty. As such, the isolation index is defined here as the Euclidean distance to the nearest island or landmass of equivalent or larger size, modified from Dahl [43].

Productivity
NPP is a commonly employed index of biological productivity based on the rate at which plants store photosynthetic energy as biomass. Island NPP values therefore correspond to the vigor of flora with consequent implications for its human and animal consumers. This has clear relevance for the sustainability of hunter-gatherers situated in the higher tiers of trophic systems, and for agropastoralism premised on exploiting cereals and pulses. More broadly, we note statistical support for a relationship between human population and high-productivity environments, albeit modified by the relationship between such environments and pathogen loads [44,45]. The NPP data used here were collected in 2019 by the Moderate Resolution Imaging Spectroradiometer (MODIS) on NASA's Terra satellite and subsequently made available at 500 m resolution in the MOD17A3HGF Version 6 product [46]; this product was downloaded through Google Earth Engine.

Biodiversity
Lastly, we employ species richness in this study as a measure of biodiversity. Species richness is calculated as the sum of all species present in an area, without consideration of species evenness, abundance, or the number of individuals. Our use of this metric derives from the assumption that greater numbers of species should provide greater opportunities for risk management and expanded diet breadth during times of stress, e.g., [45]. Beyond this, biodiversity is a strong predictor of ecosystem functioning [47], with concomitant impacts on sustainability. To measure this phenomenon, we downloaded global data on terrestrial (mammal, bird, and amphibian) primary species ranges from www.biodiversity.org, which are based on data from the International Union for Conservation of Nature [48]. Species richness was calculated using a 10 × 10 km grid, within which the total number of species present was summed.

Methodological Considerations
Clearly, these data reflect the modern environmental organization of the subset islands, not their organization in past states. While we recognize the limitations of these modern datasets for measuring environmental conditions in the past, particularly in the late Pleistocene/early and middle Holocene, most Mediterranean-wide environmental differences should be consistently distributed between the past and present [49] (p. 12124). That is, even if absolute values diverge, relative difference is likely to be fairly constant over millennial scales, at least when considering productivity values. Any effects on NPP that may be the result of the long-term island occupation that we are evaluating can also be assumed to be consistent across all islands with documented occupation. As our analysis only makes use of those islands for the inter-island variability we aim to measure, modern NPP data are a useful means of appropriately differentiating island ecologies.
Species richness is more problematic. In common with many island environments, human settlement in the insular Mediterranean has substantially reconfigured its biotic organization [32], with an 89% extinction rate of mammalian endemics at the Pleistocene-Holocene boundary [50] (pp. 194-197, Figure 9). Nonetheless, the Mediterranean remains highly biodiverse. It is also conceptually challenging to extract biodiversity values from palaeoenvironmental data, and this includes zooarchaeological and palaeoethnobotanical data: while clearly useful for determining taxonomic presence, intervening taphonomic and sampling effects nonetheless generate highly lacunose assemblages that are ineffective in reconstructing the total range of palaeodiversity. Human-mediated Holocene extinctions across the insular Mediterranean have exhibited parallel structures, however, and the same types of endemics have been extirpated and replaced with broadly homogenous, anthropically assembled ecologies.
As a consequence, we anticipated that modern species-richness data are appropriate for our purposes here. That there is no statistical relationship between area and species richness is consequently problematic, especially when other data support SARs for the insular Mediterranean, e.g., [51,52]. We speculate that the bias in our data towards terrestrial chordates and away from other taxa is skewing the relationship. At present, however, there are no extant pan-Mediterranean datasets with enough taxonomic breadth to resolve this issue. Thus, we are able to evaluate only the impact of terrestrial vertebrate biodiversity on patterning in the colonization data. We also were unable to include information pertinent to the productivity of marine ecosystems. This is partly a function of data availability, but also a function of assigning discrete values to often mobile resources, and furthermore categorically associating them with one island or another. Considering the emphasis placed on terrestrial, as opposed to maritime, food sources by initial Mediterranean farming populations, however, we do not think this exclusion has necessarily limited our analysis to any great extent.
We gathered data on hunter-gatherer and agropastoral settlement from several recent synthetic sources, most notably Dawson [28], Cherry and Leppard [33], and Leppard [23]. While this distinction between modes of subsistence is less stark in other global theaters, in the Mediterranean there is a demonstrable structural change from majority food-acquiring to majority food-producing economies during the Neolithization process [23]; thus, we argue that our hunter-gatherer versus agropastoral distinction is empirically coherent. We also note that "colonization" might best be conceptualized as a spectrum of behaviors [28] (pp. , rather than simple presence/absence, a factor that impinges on our discussion below. The upper and lower temporal bounds of our study are the Last Glacial Maximum and an arbitrary date of 1 BC/AD 1, as a point by which the vast majority of islands for which data exist had witnessed some form of agropastoral activity. Accordingly, we here avoid the controversy surrounding pre-LGM island colonization by our species and pre-Upper Palaeolithic island colonization by other hominin species [53,54]. Our analysis of both hunter-gatherer and agropastoral settlement took different forms. As regards the latter, every island in our dataset has almost certainly witnessed some form of agropastoral settlement over the last eight millennia, cf., [37], but many islands lack published archaeological dates or data. Regarding hunter-gatherer settlement, the number of islands with relevant archaeological evidence (of whatever scale or duration) is substantially smaller. This very probably reflects the fact that, for reasons we wish to evaluate, hunter-gatherer settlement is constrained to a subset of islands. Finally, we emphasize that-as we are interested in diachronic patterning-we do not discuss in detail the individual sites implicated in our analysis; the reader is directed to the relevant sources [23,28,33].
Ideally, we would assess not only colonization date as a proxy for environmental variability driving sustainability, but also permanence of occupation as meaningful to understanding environmental factors which promote sustainability. Cyclical patterns of settlement and abandonment (over whatever temporal scales) might be taken to indicate that some unknown environmental factor is intrinsically limiting for sustained settlement. There is evidence for such abandonment among agropastoral populations in the Mediterranean [28] (pp. 180-208). However, except for very small islands that have been intensively investigated, or for extended (i.e., millennial) breaks in cultural sequences, it is extremely challenging to distinguish actual abandonment from potential gaps in the archaeological record. We can fairly confidently assert that very large islands with extensive archaeological records (e.g., Sardinia) have never been entirely depopulated since initial Neolithization, and that very small yet thoroughly investigated islands (e.g., Antikythera, Palagruža) almost certainly experienced near complete abandonment episodes. Between these poles, however, it is very difficult to assert categorically that an island experienced definitive habitual abandonment. The subset of islands for which we can demonstrate abandonment or a lack of it is too small for statistical analysis, so this avenue of investigation could not be pursued.
The quality of the archaeological data itself, and especially chronostratigraphic data, is problematic. Radiocarbon assays in the insular Mediterranean in general are very patchily distributed. For example, Katsianis et al. [55] amassed a database of 3159 radiocarbon dates from 353 sites in Greece, but few dates came from the Greek islands. A majority of these only exhibit dates from single sites, and no radiocarbon dates at all have been reported for Andros, Tinos, Ikaria, Karpathos, Serifos, Astipalaia, and Kalimnos. Dates on material in a relevant stratigraphic relationship with-for example-early and middle Holocene insular Mediterranean archaeological contexts are even more sporadically distributed, with investigation not unreasonably but nonetheless often biased towards certain key sites and stratigraphic sequences (often on larger islands), e.g., [56]. Beyond problematic distribution, too frequently radiocarbon dates from the Mediterranean are bedeviled by a series of methodological limitations including inappropriate paucity within discrete deposits, incomplete publication (e.g., 48% of published Greek dates lack δ 13 C values [55]), lack of marine reservoir effect correction, and excessive interpretive weight placed on single dates with large standard deviations.
Since for many islands within our sample no relevant absolute dates exist, presumed dates are based on the morphological or stylistic affinities that cultural assemblages have with comparable assemblages situated in established dating frameworks. Often, we can have confidence that dating via morphological or stylistic parallels is robust. For example, because of detailed research on western Mediterranean impressed ceramics [57,58], the presence of this type of material in otherwise culturally agropastoral deposits is highly diagnostic temporally. In other situations, the temporal boundaries of the use of a given ware or type are not usefully constrained. The practical consequences of this are several. In some cases, based on immediate context, we might suspect seasonal hunter-gatherer exploitation, but it is not evident (e.g., in the coastal islands of the Dalmatian archipelago [59]). More concerningly, because we are particularly interested in temporal structure of Neolithization, a further result is that our analysis lacks granularity as colonization horizons can rarely be narrowed down to centennial scales. More often, we are working with a margin of error of perhaps half a millennium. For example, when a colonization horizon is reported as "fifth millennium BC", we arbitrarily select a date which is likely to be least inaccurate (e.g., 6500 BP), but we have no means to evaluate inaccuracy in absolute terms. An obvious solution would be to incorporate only those islands with radiocarbon dates, but this reduces the sample to an unacceptable size for statistical analysis. For now, beyond exhorting Mediterranean island archaeologists to acquire high-quality radiocarbon dates as a best practice, we offer no solution. Despite these limitations, our analyses show that general patterns and relationships are evident and aligned with our expectations. From this, we infer that the majority of the dates employed here are largely accurate, but also stress the potential for these models to improve further with more refined chronologies.
For the purposes of our model then, and our evaluation of how environmental and geographical factors impact the sustainability of island occupation, we take the antiquity of island habitation to be indicative of its permanence and longevity. That is, our model assumes that islands remain occupied continually from the earliest documented date of settlement. As our dates of earliest habitation are in recorded as BP, these dates thus also indicate duration of occupation. Those occupied first are therefore those occupied Sustainability 2021, 13, 3383 9 of 20 longest, providing a general, albeit coarse, means by which to evaluate the effects of various variables on island suitability.
Finally, we note here the problematic issue of palaeogeographic reconstruction and geographic changes since the LGM. Mediterranean coastlines were very different at the LGM from their present configuration, e.g., [60]. Much of our analysis here covers the last ten thousand years, but even during the Holocene coastal geography has been slightly yet significantly different, e.g., [33]. This has evidently altered the area and NPP of some of the islands under consideration and potentially affected the composition of their biotas. It has also potentially submerged sites relevant to understanding overall patterning in settlement and colonization. These factors constrain the overall analysis when considering Upper Palaeolithic and Mesolithic data, although by the initial Neolithization of the Mediterranean sea levels more closely approximated current levels.

Analysis
Having calculated environmental attributes for our subset of Mediterranean islands and collected their associated dates of settlement, we began by conducting several simple exploratory tests to examine the distribution of these variables relative to islands with documented occupation and those without. We conducted these tests, in part, to further investigate potential biases in our dataset, but also to begin to identify differences in environmental variables impacting island settlement. T-tests were used to determine if there were significant differences in the environmental variables of islands with and without documented dates for both hunter-gatherer (HG) and agropastoral (AGR) occupations.
We next sought to determine which of these environmental variables, either independently or in some combination with each other, best modeled the chronological variance in island occupation dates by hunter-gatherer and agropastoral populations. That is, although we have identified the variables reviewed above as potentially having significant influences on when islands were occupied (or not), a multiple linear regression model that uses all these variables is not necessarily one that best models the timing of island occupation. A model that uses all possible variables may fit the data better than a model with fewer variables, but in doing so may limit the predictive power of the model on other datasets. Here, the principle of parsimony is essential, by which models with the least number of variables for a maximum data fit are preferred over models with more variables [61] (pp. [31][32][33][34][35]. Through a process of model selection, we then investigated which environmental variables are most essential for modeling island occupation and identified a model that best fits our data and allows us to evaluate the relationship between island environments and occupation dates. For model selection, we generated sets of multiple linear regression models for each subsistence strategy that included all combinations of our predictor variables. The most inclusive linear regression model, for example, is that which uses area, isolation index, NPP, and species richness to predict island occupation dates, while the most exclusive uses only one of those environmental variables (e.g., area). We identified the best model of each set using Akaike's Information Criterion (AIC). AIC values are calculated based on the fit of each model and the number of variables used. The model with the lowest AIC value is that which underfits or overfits the data least, relative to other models in the set [61]. More specifically, we use second-order AIC (notated as AIC c ) due to our sample size, with the main difference being an additional bias-correction term [61] (p. 66). An additional metric we report for the best fitting model is the Akaike weight (AIC c Wt), which indicates what proportion of the total predictive power of the model sets each individual model holds. This statistic is useful for providing additional information about how well the best-fitting model predicts the data relative to other models, or in other words the probability that that model is the best of the set. Importantly, though AIC values help us identify the best available model, they provide no information about the predictive strength of that model (R 2 ) nor its statistical significance (p-value). For that, we report R 2 values from the best performing models, which indicate how much of the variance in occupation date is predicted by those models. Due to issues of multiple comparison that arise from model selection, p-values from individual models are misleading and are not considered here [61] (p. 84). We then examined the parameters and R 2 value of those models with the lowest AIC to make inferences about environmental effects on island occupation histories.

Results
Examining solely environmental/geographic differences between islands with attested early occupations and those for which there is little data or no evidence, we can make several observations. For HG occupation (Figure 2), we find that islands with documented occupation (n = 16) have a much wider range of areas than islands with no documented occupation. This is true as well for the isolation index of occupied islands. Island NPP and species richness values, on the other hand, were not significantly different between occupied islands and those with occupation only presumed or unlikely. Of these four variables, island area was the only variable with a significantly different distribution (p < 0.05) between islands with occupation and those without. Given the small sample size of islands with documented HG occupation dates, however, we must be cautious in drawing any conclusions from these results.
For AGR occupation (Figure 3), area, isolation index, and species richness all had significantly different distribution between islands with and without documented occupation. The distribution of island area and isolation index values for occupied islands was again far more variable than for islands without. Interestingly, species richness was significantly lower on islands with agropastoral occupation than those without, perhaps reflecting modern-day legacies of long-term occupation and biotic homogenization, the bias in the data towards terrestrial chordates notwithstanding. As our model only analyzes islands with documented evidence of occupation, which as described above we assume to have uniformly impacted species richness, our subsequent analysis is unaffected by these modern signatures.
The results of our linear regression models for HG island occupation indicate that the best performing model is the model with area (β = 0.16) and NPP (β = −2915.61) as independent variables (AIC c = 279.31; AIC c Wt = 0.45). The adjusted R 2 value for this linear model is 0.60, indicating a moderately strong correlation between these two environmental variables and the relative earliness of HG island occupation. This indicates that area and NPP together are the best predictors of relative antiquity of HG occupation, at least among the variables we examined. Curiously, the coefficient for NPP is negative. We suspect that, again, this may be an artifact of sample size, but we discuss the implications for late Pleistocene/early Holocene insular sustainability below.
We obtained similar results for agropastoral island occupation, with the best performing model (AIC c = 1858.11, AIC c Wt = 0.35) again being the one with area (β = 0.13) and NPP (β = 1219.93) as independent variables. In contrast to HG island occupation, however, the adjusted R 2 value for this linear model is 0.10, indicating a very poor fit between these variables and the earliness of agropastoral occupation. We also note that the coefficient for NPP in the best performing AGR model is positive, whereas it is negative for the HG model.
For both types of subsistence organization, the second-best performing model for predicting the earliness of island occupation was the one with only area as an independent variable (HG: AIC c = 280.98; AIC c Wt = 0.20; ∆AIC c = 1.67; AGR: AIC c = 1859.4; AIC c Wt = 0.19; ∆AIC c = 1.29). Given the comparatively poor fit of the best-fitting model for predicting agropastoral island occupation, we decided to investigate the otherwise well documented effect of east-west directionality in Neolithization [23,24] on agropastoral dates. We subsequently incorporated longitude as an additional spatial parameter and expanded the model set to include those additional combinations. As we suspected, the model with area (β = 0.14), NPP (β = 1725. 16), and longitude (β = 46.55) then became the best performing model (AIC c = 1855.66, AIC c Wt = 0.34), though without any notable gain in predictive power for that linear model (R 2 = 0.13). For AGR occupation (Figure 3), area, isolation index, and species richness all had significantly different distribution between islands with and without documented occupation. The distribution of island area and isolation index values for occupied islands was again far more variable than for islands without. Interestingly, species richness was significantly lower on islands with agropastoral occupation than those without, perhaps reflecting modern-day legacies of long-term occupation and biotic homogenization, the bias in the data towards terrestrial chordates notwithstanding. As our model only analyzes islands with documented evidence of occupation, which as described above we assume to have uniformly impacted species richness, our subsequent analysis is unaffected by these modern signatures. The results of our linear regression models for HG island occupation indicate that the best performing model is the model with area (β = 0.16) and NPP (β = −2915.61) as independent variables (AICc = 279.31; AICcWt = 0.45). The adjusted R 2 value for this linear model is 0.60, indicating a moderately strong correlation between these two environmental variables and the relative earliness of HG island occupation. This indicates that area and NPP together are the best predictors of relative antiquity of HG occupation, at least among the variables we examined. Curiously, the coefficient for NPP is negative. We suspect that, again, this may be an artifact of sample size, but we discuss the implications for late Pleistocene/early Holocene insular sustainability below. To assess which factors were inhibiting the performance of the best-fitting model, or whether another major environmental variable we did not consider was impinging on the data, we investigated the residuals from this best-performing model of Neolithic island occupation. These residuals reflect the amount of variability in the data remaining unexplained by the predictor variables. Mapping these residuals (Figure 4) helps illustrate which islands have agropastoral occupation dates far earlier (blue) or later (red) than our model otherwise predicts. The presence of a pronounced spatial structure to these residuals suggests that there are additional variables that remain unaccounted for in our model.

Discussion
Our analysis allows us to be moderately confident in asserting that various combinations of environmental and geographic variables have meaningful implications for long-term settlement sustainability, but also that other factors complicate the distribution of the data and contribute to patterning. We first turn to the residuals in the best performing AGR model, and then address more broadly the issue of area and NPP as these variables pertain to sustainability.

Residuals in the Agropastoral Model
In examining the distribution of residuals from the area + NPP + longitude model, we found that the most evident spatial structure to the residuals is the strong west-east gradient of values. We argue that this pattern is attributable to the behavioral specifics of the process of Mediterranean Neolithization. Cyprus and the Balearics are demonstrable outliers, situated at the eastern and western extremities of the basin. Cyprus is best considered as a component of core Pre-Pottery Neolithic settlement processes between ~12,000 and 9000 BP, a period during which the agropastoral lifestyle was essentially stalled at 32 degrees east [23]. The Balearic colonization horizon remains strangely late, for reasons that may relate to local calculi of intrinsic environmental attractiveness [27]. In both cases, specific cultural factors drive the data away from modeled expectations.
The second group of residuals we discuss can broadly be categorized as central Mediterranean small or coastal islands, ranging from the northern Adriatic to the Libyan Sea. The southernmost of these islands (i.e., in the Sicilian Strait and points south) have Neolithic settlement horizons that are much more recent than the model predicts; those to

Discussion
Our analysis allows us to be moderately confident in asserting that various combinations of environmental and geographic variables have meaningful implications for long-term settlement sustainability, but also that other factors complicate the distribution of the data and contribute to patterning. We first turn to the residuals in the best performing AGR model, and then address more broadly the issue of area and NPP as these variables pertain to sustainability.

Residuals in the Agropastoral Model
In examining the distribution of residuals from the area + NPP + longitude model, we found that the most evident spatial structure to the residuals is the strong west-east gradient of values. We argue that this pattern is attributable to the behavioral specifics of the process of Mediterranean Neolithization. Cyprus and the Balearics are demonstrable outliers, situated at the eastern and western extremities of the basin. Cyprus is best considered as a component of core Pre-Pottery Neolithic settlement processes betweeñ 12,000 and 9000 BP, a period during which the agropastoral lifestyle was essentially stalled at 32 degrees east [23]. The Balearic colonization horizon remains strangely late, for reasons that may relate to local calculi of intrinsic environmental attractiveness [27]. In both cases, specific cultural factors drive the data away from modeled expectations.
The second group of residuals we discuss can broadly be categorized as central Mediterranean small or coastal islands, ranging from the northern Adriatic to the Libyan Sea. The southernmost of these islands (i.e., in the Sicilian Strait and points south) have Neolithic settlement horizons that are much more recent than the model predicts; those to the north (i.e., in the Tyrrhenian, Ionian, and Adriatic Seas) are much older than the model predicts. This divergence away from the trend we would again interpret according to the specifics of the Neolithization process. Accumulating evidence strongly indicates that the east-to-west spread of agropastoralism was mediated by coastal travel, almost certainly involving maritime transport technology. From recent assessment of the radiocarbon dates [23], it seems likely that the Neolithic entered the southern Adriatic from the southern extremity of the Balkan Peninsula. It was then transported westwards across the Strait of Otranto and simultaneously northwards along the Dalmatian coast, with these processes essentially indistinguishable in radiocarbon terms. Its subsequent spread through the Gulf of Taranto to Sicily and then northwards was almost certainly mediated via coastal voyaging with the possible use of small islands (Ustica; the Tuscan Archipelago) as waystations between large colonist populations. It may behoove us accordingly to make a distinction between different types of island settlement: more large-scale (in the case of, for example, Corsica or Sicily) and more expedient in terms of maintaining evident long-distance connectivity (e.g., [62,63]) in the immediate aftermath of Neolithization. Conversely, the available evidence suggests that the Neolithic was not carried across the Sicilian Strait to Cape Bon in North Africa, but rather that agropastoralism arrived in the Maghreb via the Strait of Gibraltar and the Alboran Sea [64]. This accounts for the outlier status of the coastal North African islands.
Finally, we turn to the Aegean, where the majority of islands are a reasonably good fit for the model, which is skewed towards Aegean dynamics by virtue of how overrepresented islands of that region are in our dataset (84 of 161, or 52%). Nonetheless, two exceptions are worthy of comment. Thera is precociously early according to the model; but this may reflect the massive change in gross area the island experienced following the eruption in the early fourth millennium BP, or its key location connecting the Cyclades to Crete, also settled very early. We reiterate the paucity of radiocarbon dates for the smaller Greek islands, and the reliance on morpho-typological dates with large margins of error. Large age ranges may drive residual status in the southern Sporades and Ikaria, for example. Beyond that, some vague yet potentially interesting patterns emerge, including a slightly later-than-expected horizon for the southern Cyclades and the islands of the "Western String" from Kea to Melos, while Naxos, Paros, and Antiparos experienced earlier than anticipated horizons.

Area and Long-Term Sustainability
For both the HG and AGR datasets, the best-performing models included area as an independent variable. As noted earlier, we have known that area was important as an underlying factor in the patterning of agropastoral island colonization since the work of Cherry [30][31][32]. It has also been argued to be similarly important with regard to huntergatherer settlement [33]. However, the statistical absence of isolation indices in accounting for the distribution of the data suggest that this is not related to the impact that area has on target/distance ratios. That is, it is not the case that larger areas make islands more discoverable from an island biogeographic perspective, but rather that larger areas are preferable from a settlement perspective, for both hunter-gatherers and agropastoralists. Why does area exercise such influence on sustainable settlement? Cherry and Leppard [33] emphasized the elevated trophic level and high calorific demands of humans practicing hunter-gatherer subsistence. In non-human species, the populations of taxa that occupy apical trophic positions are intrinsically constrained by area, to the extent that area constrains population size in their prey taxa, e.g., [40]. They consequently suggested that only very large Mediterranean islands-in the absence of the kinds of rich marine resources which facilitate insular hunter-gatherer settlement at high latitudes-would be able to support hunter-gatherer groups on a permanent basis. Where small islands do provide evidence for hunter-gatherer settlement, they are usually so close to larger landmasses that the insular evidence likely represents exploitation of only part of the total range of these groups. Area allows for larger and more complex trophic systems. That area also explains relative antiquity within the HG data lends further support to an intrinsic constraining link between scale of trophic systems and HG settlement. This is because the late Pleistocene to early Holocene transition in general witnessed expanding breadth in hunter-gatherer subsistence [21,22], including increasing reliance on taxa lower in trophic structures. Expanding dietary breadth should decrease the area necessary to support a given population of hunter-gatherers, other factors being equal; and this is reflected in our data, as earlier settlement of larger islands is followed by later settlement of smaller (though still relatively large) islands.
Our results lend strong credence to area as a major factor affecting long term sustainability for insular hunter-gatherers. The negative coefficient for NPP in the best performing model is, consequently, counter-intuitive, considering the importance of this variable elsewhere (44)(45). Other global contexts with evidence for long-term sustainability in HG settlement most often include very large islands, low-latitude islands with very high productivity indices, or both, e.g., [8]. The general exception to this is HG settlement of islands at mid-or high latitudes with high marine productivity; massive marine trophic systems founded either on seasonal plankton blooms or specific conditions (e.g., the Humboldt current; kelp ecologies) support large populations of marine mammals and/or fish occupying apical positions, in turn sustaining human subsistence in very low-productivity terrestrial environments. We suppose then that our NPP results are either a function of small sample size, or reflect that NPP varies so insubstantially between the contexts in question (being distributed primarily longitudinally and not latitudinally) that the corresponding effect on trophic scale and complexity is entirely negated by differences in area. In any case, that area and NPP exercise profound constraint on sustainability in HG insular settlement, to the extent that we suppose a dynamic threshold below which such settlement is inherently unviable, is a proposition that can be falsified by further testing, and it will be desirable in future work to ascertain statistically the degree of patterning in global HG insular settlement.
Area and NPP (with a positive coefficient) as determinants of sustainability in agropastoral settlement is a more complex proposition. As we have seen, culturally specific processes (and in particular the specifics of the spread of the agropastoral lifestyle throughout the basin) underpin the considerable variance in the data. Nonetheless, despite the noise introduced by contextual and cultural factors, island area and NPP still seem to exercise influence over antiquity of settlement. Considering that cereal/pulse and ungulate agropastoralism essentially represents an artificial compression of trophic structure, massively expanding the calorific potential of a landscape for human populations, area should not in principle continue to exercise such substantial constraint over agropastoral settlement. This is supported by cross-cultural evidence for agropastoral populations subsisting over the long term on some very small islands, for example in the Pacific or Caribbean, e.g., [65,66]. The role of area and NPP accordingly requires explanation. We suggest that area and NPP play a role in determining sustainability in agropastoral settlement for related yet slightly different reasons than those pertaining to HG settlement. The structuring principle we invoke for HG dynamics is essentially that population dynamics are intrinsically sensitive to productivity and thereby area; in the HG case, we are dealing with population dynamics among primary producers and low-ranked heterotrophs imposing limits on apical taxa. In the agropastoral context, food-production frees human populations from the absolute limits operating under food-acquisition strategies [67], but agropastoralists nonetheless remain subject to non-taxon-specific factors related to area, including intrinsic vulnerability to stochastic effects in small populations. Put simply, if deleterious inbreeding is to be avoided (e.g., [68,69]), human populations require reliable access to robust meta-populations. Palaeogenomic, bioarchaeological, and archaeological data suggest small and highly dispersed populations during Mediterranean Neolithization [23]. Overall fitness in these populations would be increased by retaining access to larger meta-populations and decreased by lowered access. Maritime movement, however, was a comparatively high-risk behavior in the Mediterranean. Better, then, to have access to a wider demographic pool on the same island, and not beyond it. This would favor settlement of islands capable of supporting bigger populations, with these large and stable populations proving sustainable over millennial scales. Assessing degree of persistence of populations is challenging, but island-specific palaeogenomic data hint at long-term meta-stability in populations-for example, an essential lack of turnover on Sardinia from Neolithization to the Iron Age [70,71].

Conclusions
Hunting-and-gathering and, subsequently, agropastoral subsistence persisted (and the latter still persists) over millennial scales within the Mediterranean. Populations practicing these subsistence strategies had both the knowledge and the means to access and, in principle, occupy any of the islands within the basin, islands which vary enormously in their environmental organization. They did so in a highly discontinuous fashion. We suggest that this discontinuity can be understood to reflect the capacity of various types of environmental and geographic organization to affect long-term sustainability of settlement. Accordingly, we analyzed which of five variables, in which combination, most effectively account for temporal and spatial patterning in both HG and AGR settlement. We show that area and NPP perform best when explaining the distribution of the HG data; and that area, NPP, and longitude perform best when explaining the distribution of the AGR data, although this relationship has a very low R 2 value of 0.13. We can in part explain the poor fit of this model by considering contextual effects in the residual data.
We conclude that area affects long-term sustainability in both HG and AGR populations, with NPP also affecting the latter. We argue that area is important not because it promotes discovery during colonization processes, but because of its subsidiary effects. Specifically, it permits the operation of larger trophic systems, and we suspect this impinges directly on the viability of insular HG settlement. For agropastoral settlement, in which smaller areas can support much greater demographic loads, we propose it is the relationship between area and population size that renders settlement on larger islands more sustainable, thereby biasing the data towards larger islands earlier in the Holocene. In essence, we are describing a form of Allee effect, by which increases in population and co-occupation lead to increases in individual fitness up to a certain threshold [72]. On a given island, however, demographic growth should eventually increase pressure such that the Allee benefits of population size are offset. It would be of interest to establish whether permanent settlement of smaller islands corresponds to Ideal Free Distribution-type predictions as larger islands become less valued niches (cf., [73][74][75]).
Considering how the transition from food-acquiring to food-producing subsistence behavior altered the parameters for sustainable island settlement, it is clear that the transition to food production removes threshold limitations which otherwise prohibit HG settlement on smaller islands. Supposing that one subsistence mode is intrinsically more or less "sustainable" is essentially meaningless, however, in that area and primary productivity have implications for both subsistence strategies. From a resilience perspective, we can suggest that smaller population sizes render hunter-gatherer groups less resilient, but that large agropastoral demographic loads (already potentially predisposed towards boom-andbust dynamics [76]) in highly sensitive insular contexts also risk low resilience because they drive demographic inelasticity. To that extent, we predict that islands most clearly resistant to HG settlement were also those that exhibit the best evidence for periodical AGR abandonment, probably a consequence of boom-and-bust dynamics intersecting with changing environmental parameters. In future, we aim to assess periodic abandonment in more detail.
In this analysis, sustainability in island settlement essentially becomes a function of robustness in food acquisition and production strategies and in demography. This is, in part, a result of the latitudinal position of the Mediterranean, its proximity to the core of Southwest Asian domestication processes, the average size of its islands, and their relative proximity to the continental African and Eurasian landmasses. There are clearly broader implications, however, especially for insular contexts which parallel the Mediterranean, either in terms of climate [25] or geographic organization [74]. More generally, we suggest that the multiple potential subsidiary effects of area-as promoting discoverability, population size, and resulting environmental dynamics that affect sustainability-need assessing separately in order to build larger, non-specific models of colonization. To provide a non-Mediterranean example, we need to be able to differentiate whether Honshu in the Japanese archipelago was settled comparatively early because it is easy to locate; can drive Allee dynamics; has orography-promoted robust hydrological systems; or a mixture of all three factors. Current theory would struggle to choose between these scenarios, so the type of analysis we conduct here should prove productive over the longer term.
Our analysis of the ancient patterns of human settlement and sustainability on Mediterranean islands also has implications for the future sustainability of human settlements, not only regionally but globally. Climate change and attendant anthropogenic processes threaten Mediterranean communities-not least via rising sea levels, biodiversity loss, and desertification of the southern fringes of the basin [77][78][79]. Our foregoing analysis is germane to understanding the likely human impacts of these processes in two ways. First, it indicates which communities are more exposed to converging environmental pressures-specifically those on smaller islands, which are intrinsically more sensitive to ecodemographic crunches and the loss of ecosystem services associated with biodiversity loss. Second, and in a more positive vein, it reinforces lessons drawn in other global contexts-that curation of natural systems, while an end in itself, has clear benefits from the perspective of long-term sustainability, boosting resilience in ecosystem services [80,81]. The heterogenous, highly biodiverse Mediterranean has proven capable of supporting contextually large-scale human populations for much of the Holocene, and all 161 islands in this analysis remain occupied to the present day. Substantive and immediate mitigation of the causes and effects of climate change and associated processes could, in principle, continue to preserve this capacity over millennial scales.
Supplementary Materials: The following are available online at https://www.mdpi.com/2071-105 0/13/6/3383/s1, Supplemental materials for this paper include a detailed RMarkdown document outlining all steps of analysis and figure creation (with the exception of Figures S1 and S4, which were created in QGIS 3.14), as well as additional data tables and plots showing model variables and AIC and model results. Alongside the RMarkdown document are additional files necessary to replicate this study, including tables containing chronological information for island occupations and raster layers used to calculate NPP and species richness.

Data Availability Statement:
The data presented in this study are openly available at https://osf. io/98gxh/.