Received 1 December 2015; accepted 1 February 2016; published 4 February 2016
The Puna de Jujuy is a characteristic region of Northwest Argentina; it is a high plateau with altitudes varying between 3500 and 4500 masl, cold and dry climate and wide temperature variations during day and night (Paoli et al., 2002) . The archaeological information available shows the presence of hunter-gatherers since the Early and Middle Holocene, but the evidence is most apparent in the Late Holocene (Aschero, 2011; Yacobaccio, 2011) . The emergence of farmer-pastoral societies constitutes the basis for the development of farmer-potter societies, which reach their peak in the period of Regional Development or Late period (Albeck & Ruiz, 2003; Albeck, 2007) .
The structure of a population may be determined by historical, political, social, economic, biological and geographical factors, all of which interact in complex ways with the processes of genetic drift, natural selection and gene flow, determining the temporal and spatial differentiation of subpopulations. We can approach the analysis of this structure by applying direct and indirect models, also known as bound and free models (Howells, 1973; Relethford, 1980; Relethford & Lees, 1982; Relethford, 2007) . These types of studies have been performed using gene frequencies, surnames, pedigrees and quantitative traits (Morton et al., 1968; Relethford et al., 1981; Varela, 1997) .
Several authors have made contributions to the study and estimation of different parameters in order to understand the structure of the population from quantitative traits; among these authors: Morton et al., 1968, Morton et al., 1971, Morton, 1973, 1975; Relethford, 1980, 1988, 1991a, 1991b, 1994; Relethford et al., 1981, Relethford & Lees, 1982; Relethford & Crawford, 1995 .
Following the explanation above, the aim of this paper is to study the genetic structure of the late population of the Puna de Jujuy using cranial metric traits. This involves knowing the existing biological variation within and among different subpopulations and the evolutionary factors that explain this variation.
One of the first contributions to establish the microevolutionary processes in the South Central Andes region consists of studying the temporal variation of five prehistoric groups from Northern coastal of Chile through cranial metric traits and the calculation of Mahalanobis’ D2 distance (Rothhammer et al., 1982) . The authors use a model that predicts an exponential increase of kinship over time. The model is applied to both metric and non-metric traits in populations from the Valley of Azapa and Pisagua, Chile (Cocilovo, 1995; Cocilovo & Rothhammer, 1996a, 1996b, 1999; Cocilovo & Varela, 1998) . For the prehispanic population of San Pedro de Atacama, the authors study the population structure using free and bound models in order to estimate kinship based phenotypic and chronological distance (Varela, 1997; Varela & Cocilovo, 2000) .
The method suggested by Relethford & Blangero (1990) for quantitative traits proposes a balance between gene flow and genetic drift, and allows us to estimate Wright’s FST as an indicator of genetic divergence between groups. The key assumptions of this model are to consider characters to be selectively neutral, that there is a similar mutation rate for all populations, that migration is constant and that the number of migrants is similar between any pair of populations. The deviations in relation to the model may be accounted for by variations in the effective size or by differences in the gene flow rate (Relethford & Blangero, 1990; Relethford, 1994, 1996; Relethford & Harpending, 1994; Relethford et al., 1997) . The model is applied to prehispanic populations from San Pedro de Atacama (Varela & Cocilovo, 2009, 2011) and Arica (Cocilovo et al., 2001; Varela & Cocilovo, 2002; Varela et al., 2004; Varela et al., 2006) in Northern Chile, with the aim of studying the temporal and spatial genetic divergence between the different periods of cultural development.
For Northwest Argentina, a research was carried out in the Calchaquí Valley which proposes the existence of a biological unit in the Late period, as a result of a common biological history and a shared sociopolitical and cultural organization (Baffi & Cocilovo, 1989-1990) . In the Quebrada de Humahuaca, it was possible to prove the existence of a biological heterogeneity higher than expected, which reflects the intervention of different entities in the configuration of groups that inhabited the area (Bordach & Cocilovo, 1991; Cocilovo et al., 1999; Varela et al., 1999; Cocilovo et al., 2001) . Finally, in the eastern sector of the Puna de Jujuy, what is proposed is an emerging structuring of the late population and relations with the Quebrada de Humahuaca and San Pedro de Atacama (Mendonça et al., 1990-1991) .
2. Materials and Methods
The skulls sample of the Puna de Jujuy is stored at the Museum of La Plata (University of La Plata) and in the Ethnographic Museum (University of Buenos Aires). The material belongs to different archaeological sites: Doncellas, Agua Caliente, Casabindo, Río Negro, Sorcuyo and Queta (see Figure 1) and their radiocarbon dating places the population of the Puna within the Late period, 1000 - 1450 A.D (Fuchs & Varela, 2013; Fuchs, 2014) .
The sex of the individuals was determined by observing the morphological traits of the skull (Acsádi & Nemeskéri, 1970; Bass, 1981; Buisktra & Ubelaker, 1994; White & Folkens, 2005) . The contribution of ancient DNA allowed us to certify that the sex of the individuals was the same as that determined by morphological traits (Postillone et al., 2013) . The age determination was performed according to Milner et al. (2000) , Meindl & Lovejoy (1985) , Molnar (1971) , Ubelaker (1984) and White & Folkens (2005) . The artificial cranial deformation was determined according to Dembo & Imbelloni (1938) , and the following types were considered: erect tabular, oblique tabular, erect circular, oblique circular and non-deformed. The observation of the metric traits was performed according to Bass (1981) , Buikstra & Ubelaker (1994) and Comas (1966) . In this study, 27 cranial metric variables were used (Table 1).
The sample included 302 adults of both sexes, with different types of artificial deformation, and from six sites (Agua Caliente, Casabindo, Doncellas, Queta, Río Negro and Sorcuyo) (see Table 2). In order to study the structure of the local population, the data were corrected for sex and deformation, using the general linear model (Seber, 1984) , i.e., we use the residuals of the regression between the dependent variables (cranial metric characters) and independent categorical variables (sex and artificial deformation). Later, the discriminant analysis was performed and Mahalanobis’ D2 distances between groups were calculated. Finally, a dendrogram was made with these distances using the Ward’s method for the six localities.
Furthermore, the genetic divergence index was calculated (FST) applying the model proposed for quantitative traits by Relethford & Blangero (1990) . The model can be expressed as: E (VPi) = VPw (1 − rii)/(1 − FST), where E (VPi) is the expected average phenotypic variance of the population i, VPi is the pooled average within-group phenotypic variance, rii is the genetic distance of population i to the centroid and FST is the average genetic distance to the centroid. For this study, we used the Rmet 5.0 software to calculate the value of FST, assuming equal relative population size and heritabilities of 0.55 and 0.78. Heritability (h2) of 0.55 was used since it is the value most used in the specific literature and an h2 of 0.78 which represents the estimated value for the population of the Puna de Jujuy based on the calculation of repeatability (Fuchs et al., 2014) .
Figure 1. Sites in the Puna de Jujuy, Argentine.
Table 1. Craniometrical variables.
Table 2. Sample used to analyze the population structure.
The upper hemimatrix of Table 3 displays Mahalanobis distances and the lower hemimatrix shows the probability values for the comparison between pairs of localities. Figure 2 shows the distribution of the groups based on the first two canonical variables which account for 67% of the variability between groups. The first canonical variable separates mainly Agua Caliente and Casabindo, while the second canonical variable separates Queta and Sorcuyo.
Figure 3 shows the results of cluster analysis where Agua Caliente, Doncellas and Río Negro make up a group, Casabindo is associated later, then Sorcuyo, and finally Queta.
The results of the fixation index for h2 = 0.55 shows FST biased of 0.047, FST unbiased of 0.025 and standard error of 0.005, whereas for h2 = 0.78, values were 0.034 (FST biased), 0.013 (FST unbiased) and standard error of 0.005.
Table 4 shows the results of applying the model of balance between genetic drift and migration proposed by Relethford & Blangero (1990) : unbiased genetic distance of population i to the centroid (r(ii)), standard error (se), observed variance (VObs) and expected variance (VEsp) and the difference between both variances for heritability of 0.55 and 0.78. The sites of Queta and Sorcuyo showed less observed variance than expected, Agua Caliente
Figure 2. Locality distribution according to the centroids of the first two can- onical variables.
Figure 3. Cluster diagram among sites of the Puna de Jujuy.
Table 3. Mahalanobis’ D2 distance and values of probability for the Puna de Jujuy.
Wilks’ Lambda: 0.49313 approx. F (135, 1336) = 1.5261, p < 0.0002
showed greater observed variance than expected, and the rest of the sites presented similar values between observed and expected. Figure 4 shows the connection between the phenotypic variance and the distance of each group to the population centroid.
4. Discussion and Conclusion
The calculated distances suggest that the differentiation between localities is moderate, 60% of them are not statistically significant (Table 3). FST calculated with values of h2 = 0.78 indicates that 1% of the total genetic
Table 4. Results for the model of balance drift/migration.
variation explains the difference between localities, and the remaining 99% constitutes the existing variability within localities. These results are compared to FST calculated for other local populations, for example, the Quitor Cemetery in San Pedro de Atacama, where there is a study on the temporal genetic divergence between the Formative, Tiwanaku and Regional Development periods (Varela & Cocilovo, 2009) . For this population, results show FST values of 0.027 (with h2 = 0.8), 0.046 (with h2 = 0.55) and 0.018 (with h2 = 1). Estimated values for h2 = 0.8 and 1 in San Pedro de Atacama are very similar to the values of the Puna de Jujuy, where FST values are 0.025 (with h2 = 0.55) and 0.013 (with h2 = 0.78). Later, the phenotypic differentiation between the Early, Middle and Late periods of the ayllus of San Pedro de Atacama is analyzed. Results show a moderate process of genetic divergence between periods, with an average FST value of 0.0481. The FST for the Early, Middle and Late periods is 0.114, 0.005 and 0.026 respectively (Varela & Cocilovo, 2011) . This last FST value for the Late period is very similar to the value obtained in the population of the Puna de Jujuy.
In the prehistoric groups of the coast and the Azapa Valley (Northern Chile), evidence shows that gene flow is lower on the coast than in the valley, minimal estimated FST for the total population (coast and valley) being 0.02, 0.01 for the coastal populations, and 0.006 for valley populations (Varela & Cocilovo, 2002) . FST values for the total population and for the population of the coast of Azapa are similar to the values obtained for the Puna de Jujuy.
The results of applying the drift/migration model (Relethford & Blangero, 1990) show that the sites of Queta and Sorcuyo present observed variance values lower than expected for those localities; this can be due to a gene flow that is lower than the average. In Agua Caliente, a higher variance than expected is observed, which may reflect a higher migratory effect than the expected average. The rest of the sites show values close to the expected by the model of balance between gene flow and genetic drift.
The site of Agua Caliente (Alfaro, 1983) , because of its greater gene flow, will show less kinship within-group and minor isolation with other sites of sub-region. This site is located on the same line than Doncellas but on the opposite edge of the river, it is a conglomerate or semi-conglomerate undefended, with rectangular and quadrangular enclosures forming simple and compound units.
By contrast, the lower gene flow of Queta and Sorcuyo is related to increase within-group kinship and isolation with other groups in the sub-region. Queta has been cited by archaeologists who carry out research or visit the site, but the only available information is Boman (1908) and Madrazo & Ottonello (1966) , who define it as a conglomerate village without defense. In a later contribution, Alfaro (1988) walks around the site and can only find a partially-opened cist and three utilitarian urns with the remains of infants.
The archaeological site of Sorcuyo (Casanova, 1938) or Pueblo Viejo de Tucute (Albeck, 1999; Albeck & Ruiz, 2003) is a large semiconglomerate village located in two hills on both sides of the Tucute stream, and appears to be inhabited as from late 900 A.D. until early 1500 A.D. (Albeck & Zaburlín, 2008) . Some years ago, Albeck (2007) proposed the entry of Highland tradition groups, which settled in the Casabindo area. The estimated date of this event is previous to the year 1300 A.D., when the village of Pueblo Viejo appears to be settled permanently. These two groups may have caused greater isolation as a consequence of the decrease in genetic exchange when compared to the rest of the population, or of practicing a reproductive system with more frequent mating within-group than between groups.
The distances calculated between the six locations show a moderate phenotypic differentiation between them. The results of the model application drift/migration (Relethford & Blangero, 1990) , reveals a differential distribution of morphological variability between sites. Agua Caliente, Casabindo, Doncellas and Río Negro show higher morphological variability, product of a higher gene flow, whereas Queta and Sorcuyo show lower variability. This differential distribution of the morphological variation could be explained by the entry of groups possibly of altiplano origin, or other regions in the Late Period. The greater isolation of these two groups may be the result of a decrease in the genetic exchange regarding the rest of the population, or by practicing a reproductive system that promotes mating within the group.
We thank the institutions that financed this research: Consejo Nacional de Investigaciones Científicas, Secretaría de Ciencia y Técnica-Universidad Nacional de Río Cuarto and Agencia Nacional de Promoción Científica y Técnica. We also want to thank the following people: Myriam Tarragó, Claudia Aranda, Karina Zucala, Héctor Pucciarelli, Mariano Del Papa, Julián Valetti and Maia Upstein.