Temporal, spatial and gender-based dietary differences in middle period San Pedro de Atacama, Chile: A model-based approach

To explore the possible emergence and lived consequences of social inequality in the Atacama, we analyzed a large set (n = 288) of incredibly well preserved and contextualized human skeletons from the broad Middle Period (AD 500–1000) of the San Pedro de Atacama (Chile) oases. In this work, we explore model-based paleodietary reconstruction of the results of stable isotope analysis of human bone collagen and hydroxyapatite. The results of this modeling are used to explore local phenomena, the nature of the Middle Period, and the interaction between local situations and the larger world in which the oases were enmeshed by identifying the temporal, spatial, and biocultural correlates and dimensions of dietary difference. Our analyses revealed that: 1) over the 600-year period represented by our sample, there were significant changes in consumption patterns that may evince broad diachronic changes in the structure of Atacameño society, and 2) at/near 600 calAD, there was a possible episode of social discontinuity that manifested in significant changes in consumption practices. Additionally, while there were some differences in the level of internal dietary variability among the ayllus, once time was fully considered, none of the ayllus stood out for having a more (or less) clearly internally differentiated cuisine. Finally, sex does not appear to have been a particularly salient driver of observed dietary differences here. While we do not see any de facto evidence for complete dietary differentiation (as there is always overlap in consumption among individuals, ayllus, and time periods, and as isotopic analysis is not capable of pinpointing different foods items or preparations), there are broad aspects of dietary composition changing over time that are potentially linked to status, and foreignness. Ultimately, these stand as the clearest example of what has been termed “gastro-politics,” potentially tied to the emergence of social inequality in the San Pedro oases.


Introduction
To explore the possible emergence and lived consequences of social inequality in the Atacama, we recently completed intensive bioarchaeological and biogeochemical analysis of a large set (over 600 individuals) of incredibly well preserved and contextualized human skeletons recovered from a series of cemeteries dated to the Middle Period of the San Pedro oases. In the present work, we explore one subset of our data in detail, the model-based paleodietary reconstruction accomplished through stable isotope analysis of human bone collagen and hydroxyapatite. We use this to explore local phenomena as well as the nature of the Middle Period, stressing the importance of the indigenous responses, as well as the interaction between local situations and the larger world in which the oases were enmeshed. Given that the distribution of, and access to, varied kinds or quantities of food is a common mechanism and/or consequence of intra-societal differentiation [9,10], knowledge of how patterns of consumption may have differed over time, among oases, and between various subsets of the population under study (e.g., males and females) offers unique insights into past processes of inequality. Combining these data with a rigorous program of radiocarbon dating and detailed analysis of grave contents, among other sources of information, has the potential to reveal how resource access may have varied along multiple axes of difference in the oases. Importantly, we argue that this integrative approach comes with a precision unobtainable by other archaeological means and allows us to explore in far greater detail the complex shifts occurring in local prehistory.

Geographical and archaeological background
The Atacama Desert comprises over 100,000 km 2 of northern Chile and southern Peru between ca. 18˚and 30˚South latitude. The Pacific Ocean and the Andes, which flank the Atacama on the west and east, respectively, both contribute to the desert's hyperarid conditions. Among the world's deserts, the Atacama receives a classification of "Ea23" [11], denoting extreme aridity, a lack of seasonal precipitation, and temperatures between 10˚and 30˚C. Instrumental data indicate 1-6 mm/yr of rainfall at weather stations across the region [12], multi-year studies are consistent with measurable rain only on a decadal scale [13], and the region's aridity seems to have been a persistent feature for millennia [14], if not millions of years [15,16]. Indeed, a recent review of paleoclimatic conditions in the region suggests that aridity equal to, if not exceeding, that seen in the present prevailed throughout the Holocene [17]. The exception to this hyperaridity was the existence of a somewhat wetter time between 1000-2000 years ago, encompassing the period examined here. As discussed below, this may have facilitated the spread of agriculture and the permanent settlements that characterize the Middle Period.
Within this arid landscape, the oases of San Pedro de Atacama offer access to stable conditions for the long-term settlement of humans and have sustained agropastoral societies for millennia. The oases are located at the northern edge of the Atacama salt flat, in the alluvial delta formed by the San Pedro and Vilama Rivers, before they drain into the subsurface water table that forms the salt flat. This flatter region concentrates humidity that allows the development of natural oases which have been managed and even expanded by humans in the past. The oases have been continuously occupied for over 4000 years as populations moved down from the surrounding Puna (the cold arid plateaus of the high Andes; [18][19][20][21][22]). With the adoption of horticulture as well as arboriculture and silvopastoralism in the oases, populations aggregated into larger settlements, taking advantage of the rivers and arable land [20,23,24]. From the Formative Period (~1200 BC) onwards, there is strong evidence of aggregated settlements, large and small habitation sites, as well as substantial cemeteries that document the long-term occupation of the oases (e.g., [20,[24][25][26]). Settlement pattern analyses have revealed a substantive concentration of sites, both cemetery and habitation, in the central and northern oases during and just before the Middle Period (AD 500-1000) [20,27]. Agüero [27] argues that this demonstrates the considerable growth and ultimate stability of local villages during the Middle Period. This time also witnessed the widespread integration of the Atacama oases into the interregional networks of the southern Andes. Nu 0 ñez [21] has argued that this time also sees the earliest evidence for the production of surplus, a phenomenon that facilitates the growing caravan trade and is a key factor in the rise of persistent inequality. This is bolstered by a surge in population as well as incipient craft specialization, evidenced in new metallurgy practices, the growth of agriculture, and individuals who focused their energies on the caravan trade [21, [28][29][30][31][32].
In the San Pedro de Atacama oases, this rise in complexity and affluence appears to have manifested in different fashions across the oases and ayllus. The ayllu is the traditional Andean form of kin-based social organization that can also denote territorial boundaries (ascriptive descent groups [33][34][35]). Interestingly, Goldstein [36] has argued for the ayllu as an ethnic community that, while based on kin, is bounded territorially by huacas or ancestral shrines. For example, in 1642, an idol ("ydolo") or huaca was identified as belonging to Sequitur (today's ayllu of Sequitor); it consisted of a black stone that was worshipped by the ayllu's residents and was hidden in a type of vault. Similarly, in the ayllu of Contituques Cantal Acapana (today's Conde Duque or San Pedro) was an idol called Tocol, which was buried in the courtyard of the ayllu's church, having been previously burned in the public square [35]. Locally, ayllus have typically been interpreted as naturally differentiated oases (e.g., Solor, Tulor, or Coyo) as well as the divisions of the major oasis into internal sectors marked by territorial features and contemporary boundaries (Fig 1).

Food and intra-societal differences
Anthropologists have expended countless pages on the ways in which dietary differences can mirror, reinforce, or even establish social inequalities. Arnold [37] asserts that "food was, and continues to be, power in a most basic, tangible and inescapable form." The collusion of the basic sustaining power of food, its necessity for life, and the demands of systems of social inequality, means that "variation in what people eat reflects substantive variation in status and power and characterizes societies that are internally stratified into rich and poor, sick and healthy, developed and underdeveloped, overfed and undernourished" [38]. In this way, the consumption of, or access to, particular foodstuffs has come to mark membership in empowered or subjugated intra-societal groups, and/or as an indicator of internal social schisms.
In spite of the inherent tendency of food to homogenize difference and bring people closer together [9], many stratified societies have developed means by which this innate homogenizing effect of food is overridden through a process of semantic inversion and serves as a tool that, in Appadurai's words, can "serve to regulate rank, reify roles, and signify privileges," and which can, "sustain relations characterized by rank, distance, or segmentation," [9]. Goody [10] maintained that differential control of, and access to, food can play a central role in the creation and maintenance of multiple intra-and inter-societal hierarchies. While food may have a central social role in even the least stratified society, increasing social stratification offers the opportunity (and perhaps establishes the demand) for even more intense uses of food in the negotiation of relations of power [39]. In such societies, the power of food as a symbol or code for status differentiation can be fully exploited as a means of establishing and furthering inequality. In, perhaps, the purest expression of the furtherance of social inequality through food, elite individuals will sometimes go so far as to effectively monopolize certain foodstuffs, especially if that food is already an established symbol of their position in society, for example, the supposedly exclusive consumption of freshwater fish and wildfowl by the nobility of 15 th century England [40]. This monopolization gives rise to an internally differentiated cuisine, a food system in which not all people have equal access to all foodstuffs, as a consequence of differences in social status [10]. The specifics of the consumption of restricted foods: who can, and cannot, have access, when, how much, and under what circumstances, are the means by which elite individuals can attempt to build and maintain their higher status position within a given society. It is this use of food, as de-facto ammunition in internal conflicts over cultural, economic, or political resources, that Appadurai aptly labels "gastro-politics" [9].
That not all societies stratify their food is undoubtedly a reflection of the fact that there are many possible paths by which a society may become stratified [2], and that it is not obligatory that dietary practices follow broader patterns of social stratification, especially in societies without formal stratification. Nevertheless, it is quite often the case that stratified societies do possess some form of differentiated cuisine and the development of this differentiation seems often to be coeval with the development of early stratified societies [41].

Research questions
Building on this exploration of social differentiation and diet, our research into the paleodietary practices of the Middle Period inhabitants of San Pedro de Atacama was guided by a series of questions that sought to identify the temporal, spatial, and biocultural correlates and dimensions of dietary difference.
Time. As the individuals included in this analysis represent nearly a millennium of the occupation of the San Pedro oases, our initial interrogation of the dietary data was diachronic in nature. Specifically, using the corpus of individuals for whom we possess both isotopic and radiometric data (n = 167), we sought to determine if there were any coherent temporal trends (correlations) in the consumption of any of the modeled food groupings that might have resulted from diachronic changes in resource availability/access. Such changes in access or consumption may, in turn, be connected to the establishment, institutionalization, or intensification of social inequality in the San Pedro oases over time.
Furthermore, in the temporal realm, we sought to test the dietary dimensions of a hypothesis advanced in Pestle and colleagues [42], which proposes a disjuncture between what we term the early/incipient (pre-600 calAD) and late/established (post-600 calAD) Middle Period. As above, the large corpus of dated individuals allows for a specific examination of whether there existed dietary correlates of this proposed inflection point in San Pedro's history.
Space. In our previous examination of the chronology of the six ayllus under consideration here [42], we noted the apparent independence of the use-life patterns of different ayllus. Based on this, we contend that the dynamics of human activity in each ayllu would appear to have been particularly sensitive to local socio-cultural and environmental dynamics, a conclusion that reinforces previous archaeological and bioarchaeological findings of significant differences among ayllus during the Middle Period [43][44][45][46]. Here, we attempt to quantify any differences in dietary practices among the ayllus through comparison of modeled dietary contributions. This line of inquiry flows from the notion that differences in access to different foodstuffs among ayllus could be a consequence of different positioning vis-à-vis the interregional exchange systems of the period, or indeed of differences in positions in the status hierarchy of the ayllus. As the different ayllus sometimes belong to different portions of the Middle Period, we augmented our comparisons of the "raw" modeled foodstuff proportions with a second comparison of the residuals of modeled contributions after regression against calibrated radiocarbon dates (derived from the Inter-Ayllu date model of [42]). These time-corrected dietary differences should allow us greater insights into the mechanics of life in the various ayllus after the variable of time is fully accounted for.
Sex. Finally, working from the complementary theories that: 1) differences in diet are often associated with, or pattern on, sex and understandings of gender, and 2) that inequality (broadly construed) is experienced differently by men and women, we next examined the relationship between biological sex (as determined osteologically) and diet. First, at the coarsest level, we compared modeled diet between all males and females in the study sample (grouping all individuals from all ayllus together). Next, in light of the aforementioned differences among the ayllus, we compared male and female diet from each of the four ayllus for which we possessed sufficiently large samples of sexed individuals. This comparison was made on both raw modeled foodstuff estimates and with the date-regressed residuals of foodstuff values to control for time (see above). Finally, in an attempt to integrate considerations of both social/ spatial (ayllu-based) and sex-based dimensions of dietary difference simultaneously, we compared the modeled diets of males and females from each of the four best-represented ayllus with each other. As above, this comparison was made on raw modeled foodstuff estimates and with the date-regressed residuals of foodstuff values.
The three vectors of difference explored here are not meant to be an exhaustive list of the possible ways that inequality may have been partitioned in the San Pedro societies. However, they are the ones where the information available in the bioarchaeological record allow us to explore the differentiation on this spatial scale with a high degree of confidence. For instance, another important vector of social differentiation in many human societies is age, but this is an aspect that we are unable to explore with the current available skeletal record, since 1) the study sample included less than 6% subadult individuals, and 2) most of the individuals included in our study were represented by cranial remains alone, rendering precise age determinations too unreliable to be useful in statistical models.
Of relevance to the three aforementioned vectors of difference, in the present work, we examine diet as a proxy for underlying social, cultural, and political dynamics, remaining conscious that the observed changes in consumption are likely reflections of changes within society. Put differently, when confronted with diachronic changes in consumption, we are aware that time itself is not changing diet, but rather that something integral to the arrangement of Atacameño society is changing with time, which, in turn, is affecting diet. While such proxy evidence can be challenging given that it is not a direct measure of the phenomenon of interest (social change), the accumulated decades of anthropological work on the linkages between consumption and social structure add support to the validity of our investigation along these three lines.

Isotopic basics
Stable isotope analysis remains the best means available to archaeologists for the estimation of the diet of ancient individuals. Indeed, nearly fifty years of the application of this method in archaeology has built a corpus of experimentally validated methods that permit accurate estimation of broad aspects of human paleodiet (see [47,48] for general reviews of the method). Assuming sufficient preservation of target biomolecules (bone collagen and hydroxyapatite) and accurate knowledge of the local foodweb and fractionation (the offsets between, for instance, consumed foods and consumer tissues), stable isotope analysis can be used to estimate broad aspects of individual paleodiet (C 3 vs. C 4 plants, marine vs. terrestrial protein) with a high degree of accuracy. The basics of this method have remained effectively unaltered since the earliest archaeological applications of the technique [49,50], and the real power of the method remains in the scalability of the data it generates, which permits inference of everything from individual diets to the consumption practices of broad subsets of a past society. Given the incredible preservation of skeletal biomolecules in the Atacama Desert, and the relative simplicity of the local foodweb (which lacks marine or freshwater foods), stable isotope analysis is particularly well suited for paleodietary reconstruction in the region.
Classical applications of stable isotope analysis to archaeological questions of diet have tended to proceed as follows. Collagen (and sometimes hydroxyapatite) are extracted from human bone samples, and mass spectrometry is used to generate estimates of stable isotope ratios (expressed using the delta (δ) notation, a ratio of heavier to lighter isotope relative to an international standard) of carbon ( 13 C/ 12 C) and nitrogen ( 15 N/ 14 N) for the extracted biomolecules. These consumer values are then compared, typically on an ad hoc basis, with measured or stipulated values of potential food groups as to generate bounding (more than/less than), ordering, or rough proportional estimates of past consumption. Only rarely are the effects of fractionation, disparate macronutrient composition of different foodstuffs, and the effects of differential routing of different food fractions to different consumer tissues accounted for systematically and sufficiently.
Fortunately, recent decades have witnessed the development of powerful new tools for the interpretation of isotopic data, in particular mixing models, which permits probabilistic quantification of source contribution. Such models greatly improved the explanatory power of isotopic estimation of paleodiet because these approaches use a complex suite of user-stipulated inputs to model bounded source contributions to a given consumer's diet, rather than relying on basic visual comparisons of plotted source and consumer data, or making simplistic assumptions (for example, that hydroxyapatite isotope values reflect plant diet). While the earliest mixing models [51,52] could only solve for a limited number (n+1) of food sources, later developments [53][54][55] expanded calculations to non-determined systems, and provided means for dealing with uncertainty and incorporating priors based on a Bayesian approach.
Most recently, models such as FRUITS [56] have been built to accommodate the input of consumer, foodweb, macronutrient composition, and routing parameters, as well as user-specified priors, to determine, probabilistically, bounded source contributions while accounting for uncertainty in all input data. These recent Bayesian approaches "offer a powerful means to interpret data because they can incorporate prior information, integrate across sources of uncertainty and explicitly compare the strength of support for competing models or parameter values" [53]. The substantial potential of these methods for addressing questions of significant paleodietary interest has consequently resulted in their increased use in the South-Central Andes (e.g., [57][58][59][60][61][62]).

Skeletal collections
The skeletal sample included in this study represent individuals excavated from 13 cemeteries dated to the Middle Period. The sites were excavated over several decades by different researchers, including the amateur archaeologist Father Gustavo LePaige, who recovered most of the skeletal collections in the region. A detailed review of the history of research and the preservation status of the collections can be found in [24,63]. The chronological contextualization of the cemeteries has been detailed in [42].
The accumulated body of archaeological work on the region suggests that the cemeteries were not precincts reserved for specific classes of individuals, as there is not systematic bias in sex representation in the sites [63], and as individuals of all ages are found in the burial spaces [28,64]. Therefore, they are assumed to be representative of local societies. However, the different archaeological recovery practices employed over time have resulted in a skeletal collection that is very biased towards adult individuals (which is reflected in our data), and most of the collections are represented only by human skulls, as LePaige, in particular, did not save postcranial remains from his excavations [63]. Moreover, during the military government in Chile, the skeletal collections were separated from their funerary context, resulting in the loss of individual context information for most of the collections excavated by LePaige. Therefore, detailed archaeological information and complete skeletons are available only in cemeteries excavated more recently, such as Solcor 3, Coyo 3, Quitor 6 Tardío, and Casa Parroquial. For these reasons, the analysis of vectors of inequality within cemeteries is not possible across all of our samples.
Demographic characteristics of the individuals included in this study (age at death and sex) were estimated based on the analysis of skulls and pelvis, following traditional bioarchaeological methods of analysis [65]. However, as mentioned before, for most of the samples only skulls are available, which makes the estimates of age at death unreliable [66,67]. For this reason, the age distribution of the local populations cannot be estimated properly, and this demographic aspect was not explored in the current study.
The cemeteries were grouped inside ayllus, as they represent long-standing geo-political units in the Atacama oases (see above), far predating the colonial period. Cemeteries are included in each Ayllu based on their geographic location within their boundaries (Fig 1), given that these boundaries have been stable since at least the beginning of colonial period [68], if not far earlier. In total the cemeteries represent the human presence of six ayllus, as detailed in Table 1 and Fig 1. Sampling. Sampling for isotopic analysis took place in the Museo Gustavo Le Paige in San Pedro de Atacama, Chile. We selected samples from a total of 288 individuals representing six ayllus (Conde Duque, Coyo, Larache, Quitor, Solcor, and Tchecar). The selection criteria for the individuals attempted to generate representative samples from each the cemeteries and the sexes, while at the same time minimizing the damage to complete bones. For this reason, whenever possible, samples were collected from fragmented skeletal remains. The sample composition is heavily biased towards adults, as subadults are not common in the collections excavated by LePaige. While chemical preservation reduces the sample under consideration to 257 individuals (89.2%), this still represents an enormous study sample by isotopic standards. In all cases,~1.0 g (target weight) samples of dense cortical bone were removed from available skeletal elements using a diamond cut-off wheel (Dremel #545) mounted on a handheld Dremel Rotary Tool. The rotational speed of the cutoff wheel was kept at a minimum effective speed as to avoid causing any thermal degradation of bone collagen, and wheels were sterilized between individuals to avoid cross-contamination. The element and location of sampling varied from ayllu to ayllu, and on an individual level, based on completeness and the state of preservation. All sampling was performed in consultation with museum curatorial and conservation staff. Following extraction, samples were bagged in pre-labeled sterile sample bags. Permits for sample export were granted by the Chilean Consejo de Monumentos Nacionales (Permits 3682/12, 3219/15, 4276/16, and 5084/17).
Extraction and analysis. The extraction of collagen and hydroxyapatite from human bone samples was performed at the Archaeological Stable Isotope Laboratory at the University of Miami. Initially, each sample was ground by hand using a sterilized ceramic mortar and pestle and separated, using geological screens, into different size fractions, which were stored in sterilized scintillation vials until extraction.
Collagen extraction followed the protocol of Longin [69], as modified by Pestle [70]. Briefly, 0.5 g of the 0.5-1.0 mm size fraction was weighed and placed in a 50 ml centrifuge tube. Samples first were demineralized in 30 mL of 0.2 M HCL on a constantly spinning rotator for 24 h, at which time those samples that had achieved neutral buoyancy were then rinsed to neutral

PLOS ONE
pH through a repeated process of centrifugation, decanting, and the addition of distilled water. Samples that had not yet been completely demineralized after 24 h had their acid changed and were returned to the rotator for a further 24 h period. Humic treatment consisted of adding 30 mL of 0.0625 M NaOH to each sample for a period of 20 h. After this time elapsed, the samples were again rinsed to neutral pH by the repeated application of distilled water, centrifugation, and decanting. The remaining "collagen" was then gelatinized in 10 −3 M HCL at 90˚C, and filtered using single-use 40 μm Millipore Steriflip1 vacuum filters, condensed, frozen, and freeze dried. Start and end weights were recorded (to the nearest 0.1 mg) and were used to calculate collagen yield (wt%) for each sample. Hydroxyapatite extraction followed a protocol first established in Lee-Thorp [71] and Krueger [72] and modified by Pestle [70]. Approximately 0.1 g (weighed to the nearest 0.1 mg) of the 0.125-0.25 mm fraction was placed in a 50 mL centrifuge tube. After weighing, each sample underwent a 24 h oxidation of organics using 30 mL of 50% bleach (~2.5% NaOCl). The bleach treatment was then repeated, with fresh bleach solution, for an additional 24 h period, for a total of 48 h of treatment. Samples then were rinsed to neutral pH through a repeated process of centrifugation, decanting, and addition of distilled water. Next, labile carbonates were removed by the addition of 30 mL of 0.1 M acetic acid to each centrifuge tube for a total of four hours, with a 5 min vacuum treatment after 2 h. After the acid treatment, each sample was rinsed to neutral pH before being placed in a 50˚C oven overnight to dry the resulting product. Start and end weights were recorded for all hydroxyapatite samples (weighed to the nearest 0.1 mg) and used to calculate the hydroxyapatite yield (wt%).
Collagen and hydroxyapatite isotopic analysis was performed in the Marine Geology and Geophysics Stable Isotope Laboratory at the Rosenstiel School of Marine and Atmospheric Science, University of Miami. Collagen samples were packed into tin capsules and analyzed using a PDZ Europa ANCA-GSL elemental analyzer interfaced to a PDZ Europa 20-20 isotope ratio mass spectrometer (IRMS). This analytical process yields information on elemental carbon and nitrogen composition as well as relative abundance of the stable isotopes of carbon and nitrogen, data which are used in the generation of δ 13 C co and δ 15 N co . Hydroxyapatite samples were analyzed using a Kiel-IV Carbonate Device coupled to a Thermo Finnigan DeltaPlus IRMA, providing δ 13 C ap values. Collagen results were calibrated using acetanilide and glycine, and an OCC (optically clear calcite) standard calibrated to NBS-19 was used for hydroxyapatite. Standards were analyzed in every sample set at the beginning and end of the run, as well as in-between the analyzed samples to ensure instrumental stability. Check samples of all three standards were also run as unknowns in every run to verify measurement accuracy [73]. Precision (as determined by replicate analysis of samples included in the present study) averaged 0.1‰ for δ 13 C co , 0.2‰ for δ 15 N co , and 0.1‰ for δ 13 C ap .
As all samples considered herein were extracted and analyzed using the same protocols and instruments, the issue of inter-laboratory differences in resulting stable isotope signatures [74][75][76] are not of concern.
Foodweb sampling. Paleodietary reconstruction requires knowledge of the isotopic composition of both consumers and potential foodstuffs. Thus, in addition to human bone samples, we also built a representative sample of the plants and animals that make up the local foodweb. To make the foodweb appropriate for the reconstruction of ancient dietary practices in San Pedro de Atacama we targeted, inasmuch as possible, archaeological flora and fauna from sites in the region, which would better represent the local foodweb available to the human consumers of interest here. As necessary, however, the archaeological foodweb was augmented through the collection of modern floral and faunal samples, again restricted to the immediate environs of San Pedro de Atacama. Ultimately, this foodweb comprised sixty-two samples representing four distinct groupings: beans, C 3 plants, C 4 /CAM plants, and terrestrial fauna. Unpublished values for six camelid samples, analyzed and graciously provided by Francisca Santana Sagredo were included in the average values for the terrestrial faunal grouping, but are not detailed individually here. As with human bone samples, permits for foodweb sample export were granted by the Chilean Consejo de Monumentos Nacionales (see permit numbers above).
Modeling. The Bayesian model known as FRUITS (Food Reconstruction Using Isotopic Transferred Signals; [56]) was used to quantify individual dietary composition. This multisource mixture modeling technique is one of several developed with the goal of better bounding estimates of food source contribution and, as such, offers clear advantages over both classical approaches to paleodietary reconstruction and other model-based approaches. While such an approach offers notable advantages (in that it permits statistical testing of multiple hypotheses, does not rely on ad hoc explanations of observed dietary differences, and fully accounts for the many sources of error in paleodietary reconstruction) it is far from a panacea. Indeed, such models require user stipulation of a large number of parameters, many of which can be difficult to obtain/estimate. Model outputs are only as good as the data upon which any model iteration is provided built, and while such models generate mathematically viable solutions, any model is invariably specific to the assumptions informing it, while their reliability to broader scenarios are frequently not straightforward or easy to assess. As such, their use requires detailed knowledge of both model mechanics and the many parameters needed for their generation. Therefore it is central to the defense of any model that the specifics parameters used are detailed, and the following paragraphs present the parameters used for the present application.
Only well-preserved consumer samples (human bone; collagen yield >0.5 wt%, carbon yield >4.5 wt%, nitrogen yield >0.9 wt%, atomic C/N ratio between 2.9-3.6) were included in FRUITS modelling. The consumer data for the model were based on the isotope ratios generated by IRMS. To account for fractionation, we first determined the consumer-foodstuff offset (and error) for δ 13 C co using the method of [77], which takes into account both the measured δ 13 C co and the spacing value (Δ 13 C ap-co ). Fractionation of δ 13 C ap was stipulated as 10.1±0.4‰, following Fernandes et al. [78]. Finally, for δ 15 N co we employed a trophic fractionation value of 3.6±1.2‰, as recommended by several experimental studies of omnivorous animals [79][80][81][82][83][84].
Foodweb isotope values comprised the edible portions of the sixty-two samples discussed above. Any modern data included in this reference sample had δ 13 C values corrected by +1.5‰ to account for recent fossil fuel burning [85]. Furthermore, the δ 13 C value of bone samples were adjusted by -2.0‰ to account for bone collagen-edible tissue offset [86,87]. Macronutrient (protein, carbohydrate, lipid) composition of each food group was determined by reference to the USDA National Nutrient Database for Standard Reference [88]. Elemental composition (particularly %C) of each foodstuff/macronutrient group was based on formulae provided in Morrison and colleagues [89], and differences in digestibility were accounted for following Hopkins [90].
To account for differential routing, all nitrogen in bone collagen was stipulated as coming from dietary protein, the carbon in hydroxyapatite was stipulated as reflecting all dietary carbon, and the carbon composition of bone collagen was set as reflecting a 3:1 ratio of dietary protein to energy [78]. Carbon isotope offsets between measured bulk food isotope values and the isotopic values of a foodstuff's fats (bulk δ 13 C-6‰) and carbohydrates (bulk δ 13 C+0.5‰) were based on data from Tieszen [91]. The carbon isotope signature of a measured bulk foodstuff's protein was determined using a mass-balance equation, such that a proportional/ weighted average of the δ 13 C of protein and energy (fats and carbohydrates) would equal the measured δ 13 C bulk value (corrected for the concentration of carbon in each macronutrient and foodstuff-appropriate macronutrient concentration).
For purposes of this analysis, we divided the available foodstuffs into four groups reflecting taxonomy and photosynthetic pathway (beans, C 3 plants, C 4 /CAM plants [which, in spite of photosynthetic differences, are inseparable isotopically], and terrestrial fauna). Some of the local C 3 plants, in particular leguminous trees of the family Fabaceae (which includes key local taxa such as algarrobo [Prosopis spp.] and chañar [Geoffroea decorticans]) but also the Cyperaceae, can exhibit nitrogen-fixing behaviors resulting in very low δ 15 N values. Given, however, that these N-fixing behaviors would appear to be related to the age of the tree(s) in question (at least for Fabaceae), with more mature trees yielding much higher δ 15 N values [92,93], and as we see a huge range in δ 15 N values in our samples of these taxa (from 1.1-15.1‰), we have chosen to include such taxa with the other C 3 plants, leaving beans proper (frijoles/porotos, with high protein content and very low δ 15 N values) in a distinct category.
Consumption of protein was limited to between 10 and 50% of protein as energy (using the FRUITS a priori data option), reflecting the lower and upper limit of possible human protein intake [94]. All FRUITS simulations were performed using 10,000 iterations, as recommended by the program's developers.
Given differences in foodweb composition, we note that the results of the modeling presented here are different from, and not directly comparable to, previous dietary modeling we have presented for the region (e.g., [60,61]).
Dating. Temporal control was obtained through a rigorous program of direct radiocarbon dating of bone collagen. Selection and collection of samples for dating were performed concurrently with osteological assessment of remains and the collection of samples for stable isotope analyses. AMS dating for the majority of these samples was performed at the National Science Foundation-Accelerator Mass Spectrometry Laboratory at the University of Arizona (n = 145), with the remaining twenty-nine samples dated by Beta Analytic of Miami, Florida, USA. All told, our sample included 167 individuals (65%) for whom we possessed both isotope values and direct AMS dates. The temporal structure of the dates used in our statistical considerations below are based on Bayesian modeling of available dates conducted at two scales (all ayllus and inter-ayllu), as fully detailed in [42].
Statistical analysis. All statistical analyses were performed in R [95], and, unless otherwise noted, nonparametric methods were employed. An alpha of 0.05 was employed for all analyses.
Correlations were assessed using Spearman's rank correlation coefficient (ρ), the nonparametric equivalent of Pearson's r, which has the added capability of assessing the strength of non-linear bivariate correlations. Bootstrapped Spearman's correlation analysis was performed using an R function (UncertaintyCorrelations) to address uncertainty in the calculation of correlations. The function was written by one of the authors (MH) and the script is provided in S1 File. The use of a bootstrapping approach allows for full accounting of uncertainty in both dates and modeled foodstuff contributions, a step that is not taken often enough in such diachronic analyses. Two-sample comparisons of variances were made using Levene's test for the equality of variance, and we tested paired differences using the Wilcoxon signed-rank method, the nonparametric version of the more familiar Student's t-test. Comparisons of multiple group means was accomplished using the Kruskal-Wallis one-way test of variance instead of ANOVA, with post-hoc pairwise comparisons made using the Wilcoxon signed-rank test.
In several instances, it was judged appropriate to remove the possibly confounding effect of temporal differences among individuals and groups such that dietary effects of differences in ayllu and/or sex could be assessed independently of time. In those cases, we regressed (linearly) the variables of interest (dietary contributions by food group) against median modeled calibrated dates and captured the residuals (by food group) for subsequent analysis. Both the original food group contributions and these residual values were visualized using Metric Multidimensional Scaling [96] (using the function cmdscale from R(stats)), which permitted us to envision the similarities among groups as "distances" in Cartesian space [97].

Results
The results of our analyses are presented in five parts. We begin by considering the quality of preservation of the samples under analysis. Next, we present a summary of the human (consumer) isotopic values, which is followed by a rendering of the isoscape of the local foodweb. Subsequently, we discuss the workings of the model by means of analysis of the relationship of model inputs and outputs. Finally, we turn to the primary focus of this analysis, the results of the FRUITS dietary modeling, with considerations of differences and trends on diachronic, all ayllu, inter-ayllu, and intra-ayllu scales. For the reasons stipulated above, we believe that these modeled results are significantly more useful and meaningful than raw isotope (delta) values.

Quality control
Of the 288 human bones samples analyzed, 257 (89.2%) had sample quality values above accepted cut-offs or within acceptable ranges (collagen yield >0.5 wt%, carbon yield >4.5 wt %, nitrogen yield >0.9 wt%, atomic C/N ratio between 2.9-3.6) commonly employed in archaeological studies [98][99][100]. It should be noted that for two samples (I-32 and I-62), the elemental analyzer failed to generate elemental yields, and that, lacking those data, atomic C:N ratios also could not be calculated. Based on the high collagen yields of these individuals, they were nonetheless included in later analysis and discussion. Individual data are to be found in Table 1. The summary of the sample quality data is presented in Table 2 and only the 257 samples meeting these stringent quality standards are considered/discussed henceforth. On the whole, the hyperarid conditions of the Atacama lend themselves to excellent preservation of buried organic materials (human skeletons included), and the integrity of the samples in this case is consistent with this. As has been previously observed [60,101], the only notable exception to this pristine preservation are those individuals from cemeteries located very close to water sources (e.g., Quitor 8), which seem to have been subject to frequent inundation.
Examination of the relationship between sample age and the state of preservation produced heartening results, as the negative consequences of differential preservation would seem to be minimal. As seen in Fig 2, for those 167 individuals for whom we possess both radiocarbon dates and isotopic data, there were no statistically significant correlations between the modeled median calibrated radiocarbon dates and collagen yield, apatite yield, carbon or nitrogen concentration, or atomic C:N ratio. Furthermore, as seen in Fig 3, there are no significant correlations among any of the sample quality variables (collagen yield, apatite yield, wt% C, wt% N, and atomic C:N) and the primary isotopic values of the samples (δ 13 C co , δ 15 N co , and δ 13 C ap ).
While there are weak, but still significant, correlations between collagen yield and Δ 13 C apco , on the one hand (ρ = -0.14, p = 0.03, n = 257), and Δ 13 C ap-co and atomic C:N on the other (ρ = -0.15, p = 0.02, n = 256), as Δ 13 C ap-co is a value derived from two primary isotope values

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile (δ 13 C co and δ 13 C ap ), these correlations are more likely a consequence of systematic relationships between δ 13 C co and δ 13 C ap than any real effect of preservation. The ultimate implication of these findings is that the diachronic trends in isotopic and modeled foodstuff values detailed below are not the consequence of differential preservation due to sample age, or the effects of differences in preservation on isotope values, but rather stem from bona fide differences in practices of consumption. The remaining correlations presented in these matrices are as expected in terms of directionality and strength. For instance, concerning Fig 2, collagen and apatite yields are significantly negatively correlated (ρ = -0.6, p<0.01, n = 257) and wt% C and wt% N are significantly positively correlated (ρ = 0.94, p<0.01, n = 255). Similarly, as seen in Fig 3, the strong positive correlations between δ 13 C ap and δ 13 C co (ρ = 0.88, p<0.01, n = 257) or δ 13 C co and δ 15 N co (ρ = 0.5, p<0.01, n = 257) are in line with a priori understandings of isotope systematics and the nature of the local foodweb (see below). As such, none of these correlations present any cause for concern.

Isotopes
While the primary focus of the present work is the FRUITS modeling of dietary makeup, rather than raw isotopic values, we briefly discuss those raw isotopic data below. All individual data for the three primary (δ 13 C co , δ 15 N co , and δ 13 C ap ) and one derived (Δ 13 C ap-co ) isotope systems are presented in Table 1, and summary data for the 257 well-preserved individuals

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile included in the present study are presented in Table 3. The fact that all of the primary isotope systems show ranges in excess of 8‰ speaks to substantial variability in dietary practices across time and space, variability that we contend patterns with a number of temporal, spatial, and bio-cultural differences.
As shown in Fig 3, δ 13 C co , δ 15 N co , and δ 13 C ap are all significantly positively correlated with one-another, as are δ 13 C ap and Δ 13 C ap-co . Given the positively (but not significantly) correlated δ 13 C and δ 15 N values of the overall foodweb (see below), and the fact that both collagen and

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile hydroxyapatite derive much of their carbon from shared sources, none of these relationships are particularly surprising.

Foodweb
Summary stable isotope data for the sixty-two samples included in the local foodweb are presented in Table 4 and Fig 4, with individual values for fifty-six of the samples found in Table 5

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile (Continued ) (see note above in Methods and Materials). Macronutrient concentrations and tissue-specific isotope data for these food groups are found in Table 6. Importantly, we note that the δ 13 C and δ 15 N values of the Atacama foodweb overall are not significantly correlated with one another (ρ = 0.19, p = 0.14, n = 62). There is exclusivity in the δ 13 C signatures of the two energy-rich plant food groups (C 3 plants averaging -23.7±2.0‰ for δ 13 C versus -11.2±1.6‰ among C 4 /CAM plants), and in the δ 15 N values of the two foods likely to have accounted for the majority of dietary protein (terrestrial animals, which have an average δ 15 N of 9.5±1.8‰, and beans, which average 0.7 ±3.0‰). Such exclusivity, even if it occurs in only one isotope system, is highly useful in attempts to determine the proportional contribution of these food groups to consumer diets. And while the overlapping isotopic values of certain groupings (e.g., beans and C 3 plants, which have almost completely overlapping δ 13 C values, with some greater separation in δ 15 N) may appear troubling for resolving foodstuff contributions, the inclusion of macronutrient composition in the FRUITS model nonetheless permits source differentiation. Indeed, in this example, the roughly 1:3 ratio of protein to energy in beans, which is effectively the inverse of that of the animal flesh (which has an approximately 3:1 protein to energy ratio) renders their respective dietary contributions resolvable in spite of their isotopic similarities.

Model mechanics
Finally, before presenting the results of our modeling, we briefly discuss the relationship between the isotope data used as model inputs and the modeled outputs, in an attempt to better understand and convey the functioning of these iterations of the model. As seen in Fig 5, if we consider isotope-modeled food relationships on a foodstuff-by-foodstuff basis, several evident patterns emerge. In terms of the consumption of carbohydrate/ energy-rich plant foods (C 3 and C 4 /CAM plants), there are strong correlations with the two primary carbon isotope values (δ 13 C co and δ 13 C ap ). Specifically, modeled C 3 plant consumption is significantly negatively correlated with both δ 13 C co and δ 13 C ap , whereas modeled C 4 / CAM consumption is significantly positively correlated with both of those isotopic values. Weaker, but still significant, relationships exist for both of these plant groupings and δ 15 N co (which is negatively correlated with consumption of C 3 plants and positively correlated with the proportion of C 4 /CAM plants), reflecting the far higher δ 15 N values of the C 4 /CAM grouping. Δ 13 C ap-co follows similar patterns of relationships with the modeled consumption of these two plant food groups, being negatively correlated with modeled consumption of C 3 plants and positively correlated with the proportion of C 4 /CAM plants.
Turning to the more protein rich food groups, beans and terrestrial animal meat, there is a rather straightforward set of relationships that would seem to govern the model outputs. Given their 13 C and 15 N depleted signatures, it is not surprising that modeled bean consumption is strongly negatively correlated with all three of the primary isotopic values. In contrast, all three of these values are moderately, but still significantly, positively correlated with the modeled consumption of terrestrial meat, a finding that is in line with what we know of the isotope values of that food group.
In the briefest possible terms, the model is performing as one would expect a priori, given the isotopic makeup of the locally available foods. Deviations from 1:1 correlations between isotopic inputs and modeled foodstuff contributions can be attributed to the effects of differential routing and the distinct macronutrient composition of the different food groups employed. As mentioned previously, one of the strengths of this model-based approach is that all these effects are fully and honestly accounted for.

Modeling
In light of the research questions outlined above, we consider the results of the FRUITS modeling along three main potential vectors of difference (time, geography/space, and sex).
Overall. As seen in Fig 6, the average modeled contribution of the four food groupings for all 257 individuals included in the present study hover around 25% each. In spite of this apparent similarity, these values nonetheless differ significantly, both in the aggregate (Kruskal-Wallis χ 2 = 18.2, df = 3, p<0.01), and with pairwise comparisons showing that C 3 plant consumption differs significantly from all three other groups. Moreover, these average values obscure an enormous amount of individual-level variability in the consumption of each food (3 to 9 times, depending on the food group). Indeed, for beans, mean modeled individual contributions ranged from 9.8-46.8%, for C 3 plants 8.2-46.6%, for C 4 /CAM plants 6.5-56.2%, and for terrestrial meat 9.7-35.6%. It is this broad individual level variation that we ultimately return to in our discussions below. Fig 7, based on a correlation analysis of all 167 individuals for whom we possess AMS dates and isotope values, there are statistically significant temporal trends in the consumption of all four food groups. For C 3 plants and beans, there are significant decreases over time, whereas for C 4 /CAM plants and terrestrial meat, modeled consumption increases significantly over time. While we used the median dates of the modeled calibrated ranges for these analyses, given the nearly 1:1 correlation between unmodeled and modeled dates (Fig 8), little changes if the unmodeled median values are used instead.

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile Given that there are associated uncertainties in the estimates of both radiocarbon dates and modeled foodstuff contribution values, we also performed a bootstrapped analysis of the observed correlations. As seen in Table 7, at their respective 95% confidence limits, the observed positive correlation between radiocarbon dates and C 4 /CAM plants still holds (as the 2.5% and 97.5% ρ values are still positive), as does the negative correlation between date and beans. In the case of C 3 plants, the previously observed negative correlation with radiocarbon dates is not fully supported, as, at the upper end of the 95% confidence interval, the bootstrapped correlation values extend just beyond/above zero (ρ = 0.01). The same is true with terrestrial meat, as the previously observed positive temporal correlation extends just slightly into negative values (ρ = -0.02) at the lower limits of the 95% confidence interval. However, given that the overwhelming majority of the bootstrapped values for both C 3 plants and terrestrial meat support the initial/original correlations, we are inclined to subscribe to their reality, although we must acknowledge that such a conclusion is not supported at a 95% confidence level.
Comparisons of modeled diet were next made between those individuals who had modeled median calibrated AMS dates that preceded and postdated 600 calAD. As seen in Fig 9, this analysis found significant differences (using Wilcoxon rank sum tests) for all four of the modeled food groupings. Pre-600 calAD individuals consumed an average of 33.4±4.7% C 3 plants
Inter-ayllu. Comparison of modeled diet on an inter-ayllu basis found significant differences at an aggregate level for all four modeled food categories, although pairwise comparisons found a somewhat inconsistent pattern of differences, as displayed in Table 8 above the diagonal.
In summary, For C 3 plants, the Kruskal-Wallis test indicated significant overall differences (χ 2 = 11.6, df = 5, and p = 0.04), and significant pairwise differences (Wilcoxon rank sum test) were found between five of fifteen pairings. In the case of C 4 /CAM plants, aggregate differences among the six ayllus were significant (χ 2 = 15.1, df = 5, and p = 0.01), and significant pairwise differences were found for seven of fifteen pairs. For beans, there was a larger aggregate difference in consumption (χ 2 = 16.6, df = 5, and p<0.01), and there were accompanying

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile significant pairwise differences among five of fifteen pairings. Finally, for terrestrial meat, we found the largest aggregate difference (χ 2 = 22.3, df = 5, and p<0.01), and significant differences for seven of fifteen pairs.
A graphical summation of these differences is presented in Fig 10, in which we present the results of a multidimensional scaling analysis of all four food groups by ayllu. The relative similarity of Tchecar and Solcor is striking, as is the degree of difference between Conde Duque and Larache.

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile A final consideration of the "raw" modeled foodstuff values relates to the patterns of internal variance in consumption within the different ayllus. As seen in Table 9, there were significant differences in variance (assessed using Levene's test for equality of variance) among the ayllus for three of the four modeled food groups: C 3 plants (F = 3.6, p<0.01), C 4 /CAM plants  As the different ayllus belong to different portions of the Middle Period, we augmented this comparison of the "raw" modeled foodstuff proportions with a second level of analysis considering the residuals of modeled contributions after regression against Inter-Ayllu modeled calibrated radiocarbon dates [42]. Even after accounting for time, some significant differences in modeled consumption remained, as seen in Fig 11 and with pairwise comparisons shown in in Table 8 below the diagonal.

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile For C 3 plant residuals, the Kruskal-Wallis test indicated significant differences overall (χ 2 = 20.7, df = 5, and p<0.01), and significant pairwise differences (Wilcoxon rank sum test) were found for five of fifteen pairings. In the case of C 4 /CAM plant residuals, aggregate differences among the six ayllus were significant (χ 2 = 25.7, df = 5, and p<0.01), and significant pairwise differences were found for six pairs. For bean residuals, there was only a borderline significant aggregate difference (χ 2 = 11.0, df = 5, and p = 0.05), and there were only three accompanying significant pairwise differences. Finally, for terrestrial meat, there was no significant overall difference (χ 2 = 5.3, df = 5, and p = 0.38), and only one significant pairwise difference. Having thus controlled for time, it would appear that inter-ayllu differences persist, at least for plant consumption and, potentially, for beans as well, although other differences (e.g., Conde Duque's distinctiveness) moderate once time is accounted for.
These differences are presented in graphical form in Fig 12, which shows the results of a multidimensional scaling analysis of all four food group residuals by ayllu. While the pattern, after accounting for time, differs somewhat from the raw value MDS plot by ayllu presented above, the overall structure is not markedly different. The relative proximity of Tchecar and Solcor remains a striking feature, as does the continued dissimilarity of Conde Duque versus Larache.

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile A comparison of the variances of the residuals among the six ayllus under consideration found no significant differences for any of the four modeled food groupings once differences in time were fully accounted for by regression. These results are presented in Table 10. Despite the lack of significant differences, the general pattern observed in the raw data holds true, as the individuals from Larache and Quitor show the greatest degree of individual-level variability in consumption across the food groupings, whereas Conde Duque and Tchecar show more internal consistency in dietary makeup. Having controlled for time, however, these differences are not significant.
Sex. The first level of comparison of sex-based differences in modeled consumption was made at the level of all sexed individuals (n = 179) across all ayllus. As the temporal structure of sexed and dated individuals (n = 112) does not differ significantly (Fig 13), this seemed an appropriate first comparison. As seen in Fig 14, in spite of a working presumption that female and male diets might differ significantly, statistically significant difference was only observed

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile for one of the four modeled food groups (beans). Indeed, for C 3 plants, C 4 /CAM plants, and terrestrial meat, there was substantial overlap in male and female consumption and insignificant sex-based differences in the modeled outcomes. On the other hand, females appear to have consumed significantly more beans (26.3±6.2%) than did their male counterparts (23.2 ±7.3%). This difference was not exclusive, however, as is evident in the Fig 14, which shows the substantial overlap between the proportion of beans consumed by females and males. Acknowledging that substantial variability in sex-based difference might have been obscured at the previous scale of analysis by temporal differences among the ayllus, the next level of analysis was conducted on an intra-ayllu scale, with male and female modeled diet from each ayllu compared with one-another. A sufficient number of sexed individuals were available from only four ayllus: Coyo, Larache, Quitor, and Solcor (Table 11). Conde Duque was excluded from this analysis as it comprised just seven individuals for whom we could determine sex, while Tchecar, which had fifteen sexed individuals, only included five males. We chose to include Larache in this analysis because, while the sample size was smaller-thanpreferred (n = 18), that number did include an equal number of females and males, and

PLOS ONE
Temporal, spatial, and gender-based dietary differences in middle period San Pedro de Atacama, Chile Larache is a cemetery with a number of intriguing contextual elements, upon which many arguments about the Middle Period in San Pedro have been built [102][103][104].
For three of these four ayllus, there were no significant differences in modeled consumption between females and males. Among the Coyo individuals (n = 54), the p-values for comparisons (Wilcoxon rank sum tests) of male and female diet ranged from 0.22-0.66. For Quitor
The final iteration of this sex-based analysis combined ayllu and sex, forming twelve groups, as to determine the effects of sex on diet across ayllus. Ayllu/sex pairs with significant differences are noted in Table 12 above the diagonal. For C 3 plants, the Kruskal-Wallis test found an overall significant difference among these twelve groups (χ 2 = 21.2, df = 11, and p = 0.03), and significant pairwise differences (Wilcoxon rank sum test) for nine of the sixtysix pairwise comparisons, all involving males from Larache. In the case of C 4 /CAM plants, aggregate differences among the twelve groups were significant (χ 2 = 20.3, df = 11, and p = 0.04), and pairwise comparison found significant differences for the same nine pairings as the C 3 plant analysis (Larache males being significantly different from the same eight other subgroups).
For beans, the comparison of the twelve groups again found an aggregate significant difference (χ 2 = 23.2, df = 11, and p = 0.02), and thirteen significant pairwise differences of the possible sixty-six, as noted in Table 12 above the diagonal. Finally, for terrestrial meat, we again found a significant aggregate difference (χ 2 = 22.3, df = 11, and p = 0.02), and significant differences for ten of the sixty-six pairs. A graphical summation of these differences is presented in Fig 15, in which we present the results of a multidimensional scaling analysis of the twelve ayllu-sex groups that accounts for all four modeled food groups. As presaged by the analysis presented above, there are some differences between the Quitor females and Conde Duque males, with the Larache males appearing most distinctive of all.
As before, given that the different ayllus belong to different portions of the Middle Period, we augmented this comparison of the "raw" modeled foodstuff proportions with a second level of analysis considering the residuals of modeled contributions for the ayllu-sex groups after regression against Inter-ayllu modeled calibrated radiocarbon dates. It should be noted that, for this level of analysis, the total number of groups decreased to eleven, as there were no dated males from Tchecar. In spite of accounting for time, some significant differences remained at this level of the analysis. Ayllu/sex pairs with significant differences in residual values are noted in Table 12 below the diagonal.
For C 3 plant residuals, the Kruskal-Wallis test did not find an overall significant difference among the eleven ayllu-sex groups (χ 2 = 17.8, df = 10, and p = 0.06). Nonetheless, there were significant pairwise differences for seven of fifty-five pairs. Similarly, for terrestrial meat, the aggregate differences among the eleven groups was not significant (χ 2 = 6.8, df = 10, and p = 0.74), and there was only one pairwise difference. For beans, however, the comparison of the ayllu-sex groups did find an aggregate significant difference (χ 2 = 22.6, df = 10, and p = 0.01), and twelve significant pairwise differences of a possible fifty-five. Likewise, for C 4 / CAM plants, a significant aggregate difference was again found (χ 2 = 23.6, df = 10, and p<0.01), and significant differences were found for fourteen of the fifty-five possible pairs, as noted below the diagonal in Table 12.
A graphical summation of these differences is presented in Fig 16, in which we present the results of a multidimensional scaling analysis of the eleven ayllu-sex groups that accounts for the residuals of all four modeled food groups. As above, the Larache males stand out as the most notably distinct subgroup, followed by the Conde Duque males. Additionally, it is noteworthy, as seen in Fig 17, that the females from the various ayllus form a rough "cluster" and the males from the ayllus exhibit greater variety. This suggests a greater homogeneity of timecontrolled diet among females and a broader range of dietary practices among males.

Discussion
In this work, we have used diet as a proxy for underlying social dynamics, acknowledging that any observed changes in consumption likely flow from, or reflect changes in, the structures of the societies under study (and, conversely, that food could have been used to affect certain social changes). In our consideration of the three research axes structuring this work (time, space, and sex), we highlight how the gastro-political lens of the study of paleodiet suggests new insights into the social dynamics of the Middle Period societies of San Pedro de Atacama.
Beginning with the considerations of time, our analysis revealed that: 1) over the 600-year period represented by our sample, there were significant changes in consumption patterns that may evince broad diachronic changes in the structure of Atacameño society, and 2) at/near 600 calAD, there probably was an episode of social discontinuity that manifested in significant changes in consumption practices.
Changes over time have been at the core of the archaeological discussion of the occupation of the San Pedro de Atacama oases [19,20,24,25,68,[105][106][107][108][109][110][111][112], given the long-term dynamics of local human occupation and their connections to neighboring regions. Despite this focus on changing social and biological aspects of life in the past, there has been a strong tendency to define local prehistory through discrete (non-overlapping) periods. In that context, much of the previous discussion about diet and health among past Atacameños has focused on comparisons between such periods [43,44,113,114]. In recent years, however, there has been an important shift away from typological views of the past, with archaeological [115][116][117][118][119][120] and

PLOS ONE
bioarchaeological work [42,46,60,61,[120][121][122], which has highlighted both synchronic variability and gradual historical processes that were products of intrinsic social developments, like population growth, emerging social complexity, and local cultural developments, as well as the exchange of people, resources, and ideas from other regions.
This new lens on the local past presents a more detailed and nuanced view of local prehistory, allocating agency to local actors, and situating the oases as important nodes inside the social network (or meshwork) that organized and influenced the South-Central Andes (rather than just a peripheral society within the sphere of influence of larger polities). While there is a long-term engagement in the scholarly literature with the important networks and polities of the Middle Period, our data reveals the import of local traditions and resources in shaping the human experience in the San Pedro de Atacama oases. In other words, this type of work creates a space where the discussion of the past can be seen as the result of the interplay between local and regional. The results we observe in diet over time in the Middle Period fit well within this discussion, as they show a gradual change in diet in the oases, stressing the variation that existed internal to this period.
As pertains to the observed long-term diachronic changes, the increased consumption of C 4 /CAM plants and terrestrial meat over time could both be argued to reflect a broad increase in the status and/or wealth of individuals across the oases. It is noteworthy that these classes of food have been associated with certain categories of power/privilege. While likely holding a variety of economic, ritual, political, military, and social functions throughout Andean prehistory, in later periods (in particular the Inka and colonial ages), chicha (a maize beer) is imbued with particular social and political power [123][124][125][126]. Interestingly, some of the higher C 4 values in our analyses come from Larache, where some of these individuals are interred with gold keros, the name given to a type of waisted drinking vessel. Keros are found in both the Middle Period Tiwanaku polity and the later Inca and are frequently understood as chicha drinking vessels [127]. While we cannot say with certainty that the observed increase in C 4 /CAM consumption is the result of the wider availability or distribution of chicha-indeed it could simply represent an increase in cultivation of C 4 plants/maize in the oases over the period of study-the possibility cannot be discounted, especially given the importance that chicha has as a ritual beverage across the Andes.
Similarly, there is a broad anthropological literature (from non-human primates [128] to complex human societies [129]) that associates greater meat consumption with higher status. The increasing role that both of these classes of (potentially) high-status food played over the centuries in Atacameño diet could appear to confirm the long-held narrative that describes the Middle Period as one of unparalleled prosperity [19,43,44,130] from which all or many of the inhabitants of the oases benefitted. The change in diet may reflect the impact of the growing networks of interregional exchange (e.g., [43,[120][121][122]131]), bringing not only wealth, but different dietary practices to the oases. As we have previously argued [60], incorporation into these exchange networks likely served to accelerate or accentuate intra-societal dietary variability.
This notion of significant change over the course of the Middle Period itself has and should continue to be explored via the material culture and settlement patterns of the area. In that vein, the observed dietary disjunction at/near 600 calAD would appear to have correlates in other categories of archaeological material. Both Stovel [132] and Gallardo and colleagues [133] have identified significant coeval shifts in patterns of ceramic production and decoration at/near this time that should be further explored in other forms of material culture. These parallel findings illuminate new foci for future research; it would be interesting to examine, for instance, whether 600 calAD represents an inflection point in the intensity, level, or mechanism of foreign presence in, or interactions with, the people of the San Pedro de Atacama oases.
Turning to socio-spatial considerations, our analysis of dietary makeup by ayllu revealed several interesting phenomena. First, while dietary variability was greater at Larache and Quitor and lower at Conde Duque, Coyo, and Tchecar, differences in internal variability among the ayllus decreased to insignificance once temporal differences were accounted for. Put differently, while there were some differences in the level of internal dietary variability among the ayllus, once time was fully considered, none of the ayllus stood out for having a more (or less) clearly internally differentiated cuisine. This could be a consequence of limitations of diet variability at any given time in what is, after all, an ecologically constrained environment. However, the lack of internal dietary variability (a truly internally differentiated cuisine) could also be a feature of emergent, rather than established systems of inequality. Thus, a true internally differentiated internal cuisine might only be expected in locales/cultures with established social hierarchies, such as at Moundville [41], Cahokia [134], Classic Period Oaxaca [135], or during the height of the Inca Empire [136]. Or, as Goody [10] noted to be possible, it could be the case that internally differentiated cuisine failed to appear in spite of the existence of bona fide institutionalized inequality, the former not always requiring the latter.
In contrast, at the level of both raw foodstuff contribution values and date-regressed residuals, the inhabitants of Larache (and to a lesser degree Conde Duque) stand apart from the individuals from the other four ayllus. Larache's distinctiveness in terms of dietary practice is, perhaps, attributable to the greater presence of non-local individuals among those people interred there. Previous strontium isotope studies [104,137], have identified Larache as an ayllu with a high percentage of non-local persons, with five of eighteen analyzed individuals possessing strontium signatures indicative of birth/youth spent outside of the San Pedro oases. That four of the five presumptively "foreign" individuals interred at Larache were male, and that the diet of the Larache males is the most distinctive of any ayllu-combination, supports this notion that it is the place of origin and/or continuing maintenance of non-local dietary traditions among a subset of Larache individuals that is driving the observed ayllu-based difference(s).
Larache's distinctiveness has long been noted, and it was once considered the site of a colony of elite Tiwanaku settlers in the Atacama [102,138]. While this speaks to an earlier perspective that gave primacy to the role of Tiwanaku in shaping the Atacama oases, our bioarchaeological data more clearly support Larache as "the burial place for a diverse yet culturally integrated and potentially elite and well-connected segment of the Atacameño population" [104], and not a space and people inherently tied to Tiwanaku. Given the interesting results seen in the modeling data, it is perhaps worth considering Martínez's [139] idea of interdigitated populations instead of assuming a pattern of ethnic enclaving. Using data from the Revisita de Atacama of 1683, Martínez details the interdigitated relationships between Atacama, Lípez, Tucumán, and Chicha regions, noting the significant movement between, and occasionally the permanent settlement of, members of these different groups in each of the other areas [35,139,140]. In some cases, the Revisita records over 40% of a population being from outside a given region. While the patterns we note in both strontium and carbon and nitrogen data from San Pedro de Atacama are not so pronounced, the data from Larache do add credibility to the idea that there was an interdigitation of the Atacama oases with neighboring areas long before the colonial incursion. Similarly, it supports the argument that important cultural interactions were not limited to a monolithic engagement with the Tiwanaku polity, but rather with numerous groups scattered across the highlands. If we were to consider the idea of interdigitation more systematically, we might be able to explore hypotheses grounded in dietary and morphological variability and diversity. This cursory consideration of our paleodietary data bear this out with evidence of variety internal to local diet, and more specifically, a varied diet at Larache.
In reference to our third structuring question, sex does not appear to have been a particularly salient driver of observed dietary differences among the Middle Period inhabitants of San Pedro de Atacama. This finding contradicts earlier studies that suggested important differences in access to resources according to sex during the Middle Period [43,44]. The lack of observed sex-based differences is surprising given that reconstructions of human movement in this period are tied to interregional exchange driven by llama caravans which have been understood, based on ethnographic and ethnohistoric evidence, to be primarily led by male kin [141][142][143], which would lead one to predict that the diets of males should show greater variability (see below). The lack of definitive sex-based differences raises the possibility that diet on these excursions was closely related to diet in the oases, at least as visible isotopically. This is especially true given that most exchange seems to have been focused on other highland communities and not on the more dietarily variable coast.
Importantly, our data suggest that diet-or at least access to particular foodstuffs as we cannot explore differences in quantities of foods-was not necessarily inherently differentiated by sex, notwithstanding our understanding of the breadth of Andean prehistory, which suggests the importance of duality in Andean cosmology. Our data suggest that males and females were not counterpointed in terms of their access to dietary resources in these Atacama populations despite some established sex differences in other practices and material culture distributions [121,144,145]. We would, however, note that our results do not entirely preclude the possibility of sex-linked dietary differences on a broader scale, in that females show greater homogeneity of dietary residuals across different ayllus, with their male counterparts exhibiting greater variety and consequently a possible range of different dietary practices (see Fig 17). As concerns this, our analyses are noteworthy, in that only Larache has significant male-female differences in modeled dietary contributions (observed in three of four food groups). Given the foreign origins of many of that ayllu's males, it is not clear whether the observed modeled differences are a consequence of dietary practices tied to sex or they speak more clearly to distinctions based on place of origin.

Conclusion
In closing, we return to Appadurai's [9] notion of gastro-politics, and ask whether any of these results could be taken as revealing the dietary dimensions of the construction or maintenance of meaningful institutional social difference within the Middle Period societies of San Pedro de Atacama. While we do not see any de facto evidence for complete dietary differentiation (as there is always overlap in consumption between/among individuals, ayllus, and time periods, and as isotopic analysis is not capable of pinpointing specifically different foods or preparations), there are broad, potentially status-linked, changes in diet over time, with foreignness (and perhaps male-ness) playing a role in determining dietary composition.
The three levels of analysis presented here contribute directly to the recent archaeological and bioarchaeological debates about social organization in San Pedro de Atacama populations during the Middle Period. On a larger level, the analyses across ayllu show a clear diachronic trend in consumption (more C 4 /CAM plants and more terrestrial meat with time), suggesting a rather uniform increase in status-linked diet, if the literature regarding the sumptuary dimensions of these foods are to be believed. This pattern fits well with the traditional view of the Middle Period in San Pedro de Atacama as a time of increasing social complexity [19,29,68], highlighting the regional commonalities among the different ayllus. This temporal shift in diet militates against the existence of a form of zero-sum gastro-politics, in that all diets are improving, irrespective of the inequalities that exist within local societies.
The more focused levels of analysis presented here, however, demonstrate the complexities of trying to define the impact that the emergent local social complexity had on dietary composition. Across all the ayllus, the results do not support a clear pattern of sexual differentiation of diet composition, with a suggestion that male diet may have been more diverse, contrasting with some previous studies of specific cemeteries [43,60,61]. Nevertheless, there are exceptions to this on an ayllu by ayllu basis, which adds support to recent discussions about the importance of local identity between ayllus and larger regional heterogeneities [46,117,120,121].
Along this line, in the males of Larache we see potential evidence for the intersection of a sex and origin-based difference in diet given that males, who were more often non-local in origin, consumed more C 4 /CAM plants and terrestrial meat. Again, if the literature on these classes of foods is to be believed, the distinctiveness of the Larache male diet could be seen as a testament to their higher status/access to prestige resources. That these were non-local individuals, who presumably had more direct access to the long-distance exchange networks for which the Middle Period is known, speaks to their differentiation being a possible manifestation of a network based strategy (a la Blanton et al. [2]), whereby status gained through network connections beyond the local conferred status upon aspirant leaders, who signified/reinforced their positioning through the consumption of different kinds/combinations of food. The richness of the material assemblage from a number of Larache burials (which includes gold keros and other rarities), provides a confirmation of the uniqueness and prestige of this subset of individuals. This may stand as the clearest example of the kinds of "gastro-politics" we set out to identify.
Supporting information S1 File. R Script for calculation of bivariate correlations when either/both variables have associated uncertainties. (R)