Ancient and modern mitogenomes from Central Argentina: new insights into population continuity, temporal depth and migration in South America

Abstract

The inverted triangle shape of South America places Argentina territory as a geographical crossroads between the two principal peopling streams that followed either the Pacific or the Atlantic coasts, which could have then merged in Central Argentina (CA). Although the genetic diversity from this region is therefore crucial to decipher past population movements in South America, its characterization has been overlooked so far. We report 92 modern and 22 ancient mitogenomes spanning a temporal range of 5000 years, which were compared with a large set of previously reported data. Leveraging this dataset representative of the mitochondrial diversity of the subcontinent, we investigate the maternal history of CA populations within a wider geographical context. We describe a large number of novel clades within the mitochondrial DNA tree, thus providing new phylogenetic interpretations for South America. We also identify several local clades of great temporal depth with continuity until the present time, which stem directly from the founder haplotypes, suggesting that they originated in the region and expanded from there. Moreover, the presence of lineages characteristic of other South American regions reveals the existence of gene flow to CA. Finally, we report some lineages with discontinuous distribution across the Americas, which suggest the persistence of relic lineages likely linked to the first population arrivals. The present study represents to date the most exhaustive attempt to elaborate a Native American genetic map from modern and ancient complete mitochondrial genomes in Argentina and provides relevant information about the general process of settlement in South America.

Introduction

The peopling of the Americas has been extensively studied and discussed from a range of different disciplines including linguistics, history, archeology and population genetics, among others. For more than three decades, mitochondrial DNA (mtDNA) has been widely used as a powerful tool to address this issue (1–4). As a result, there is a growing body of evidence that describes the existing variability within South America of the Native American-specific haplogroups A2, B2, C1, D1 and D4h3a (5–10). While studies of the mtDNA control region (CR) are still the most prevalent in the literature, the development of high-throughput sequencing technologies allows the sequencing of complete mtDNA genomes in ancient and modern samples. The identification of locally differentiated sub-haplogroups from whole mitochondrial genome sequences increases the possibility to interpret population differentiation within the continent and to develop demographic models at a regional scale.

It has been suggested that South America was initially populated following two main routes along the Pacific and Atlantic coastlines (11) (and bibliography therein). On the one hand, the earliest archeological sites are located near the Pacific coast and date to the Late Pleistocene–Early Holocene transition as documented by the sites in Chile, Peru and Ecuador (11–16). There is also archeological evidence of an early occupation of the Andean highlands (17) that could be interpreted as a signal of an additional inland Andean route. However, the archeological record and ancient DNA (aDNA) studies show an interconnection between highlands and coastal groups (17,18), suggesting that the Western Andean occupation occurred through the Pacific route. On the other side, the earliest South American sites near the Atlantic are located in the middle courses of the rivers, indicating the entry into the inland following the river shores (19) through least cost routes (20,21). The rise of sea level after the last glaciation displaced the Atlantic coastline tens or hundreds of kilometers inland, raising the possibility that the earlier sites on the east coast are today under water. Evidence of this early human occupation in the Southern Cone of South America is registered in the Argentinean Pampas region, in the site Arroyo Seco 2, which is dated to 12 170 ¹⁴C ybp (14 060 cal ybp) (22). The inverted triangle shape of the South American continent places Argentina at the meeting point between these two principal peopling streams that followed the Atlantic and the Pacific routes.

This research addresses, for the first time, the analysis of modern and ancient mitogenomes from Central Argentina (CA), defined in this paper as the territory encompassing Córdoba province, the Sierras in San Luis province and the southwest of Santiago del Estero province (Fig. 1). Given its geographical position as a crossroads, CA has a crucial importance for understanding the demographic processes that took place during the peopling of the Southern Cone of South America. It may further have served as a major route of dispersal to the surrounding geographic areas and, consequently, maintained gene flow with those populations as suggested by archeological (23), morphological (21,24,25) and molecular evidence (6,26,27).

Figure 1

Locations of the samples analyzed in the present study. The grey circles represent modern populations, with the diameter being proportional to the number of individuals analyzed (also given as text). The black circles correspond to the approximated location of archeological samples, labeled with the sample ID. References for populations abbreviated names: Santiago del Estero province: Sumampa (SUM), Villa Atamisqui (VA); Córdoba province: Córdoba capital (CBA), Amboy (AMB), Cruz Alta (CA), La Carlota (LCA), Jovita (JOV), Villa de Soto (VSO), La Para (LPA), San Carlos Minas (SCM), Villa Dolores (VDO), Chancaní (CHA), Río Cuarto (RIV); San Luis province: La Toma (LTM), Santa Rosa del Conlara (SRC), Tilisarao (TIL).

Archeological research has revealed that the CA region has been inhabited since the Pleistocene–Holocene transition, at least 11 000 ybp (28,29). These records are characterized by the presence of fishtail projectile points, a lithic design that is widespread in the early occupations of the eastern and western slopes of South America (30–32). After a period of scarce archeological evidence, the signal of human occupation reappears by 8000 ybp (33), characterized by sites with a technology of lanceolate or foliate projectile points, known as ‘Ayampitín’, which has been identified as a signal of Andean influence (28,34,35). During the Middle Holocene (8000–5500 ybp), a series of changes appear in the archeological record: reduction in the size of lithic projectile points and acquisition of a triangular shape (36,37), decrease in mobility, intensification of production and diversification of residential site types in different regions, suggesting an increase in population (38,39). By 3000 ybp, these changes were adopted throughout the region, albeit with noticeable regional variations (40). For this period, foreign rocks have been identified as prime material for manufacturing tools (38), and this, combined with the appearance of clam shells in burials (41), indicates the development of long-distance interaction networks. Between 1500 and 1000 ybp, a process of diversification in the way of life arose among CA communities in terms of exploitation of resources, occupation of environments and population growth. Moreover, biological changes were observed from the analysis of craniofacial morphology (42), secular variation in height (43) and mitochondrial haplogroups of ancient populations (44). Between 1200 and 1000 ybp, the practice of small-scale cultivation was established (45) as a complement to foraging. During the last centuries before the Spanish conquest, there were signals of stress owing to population growth, resource scarcity or environmental pressure, all of which continued later with the arrival of the Spaniards (46–48). In the colonial period, many populations of CA settled along the so-called ‘Camino Real’, which was the main route to ‘Alto Perú’ (present-day Andean Bolivia). This major passage for people and goods prompted important socio-cultural and demographic exchanges. The CA populations along the ‘Camino Real’ occupied a strategic position for trade, both legal and illegal, between Buenos Aires, Chile and ‘Alto Peru’ (49).

The evolutionary history of present-day human populations of CA has been investigated using non-recombinant uniparental markers (mtDNA and Y chromosome) (26,27,50) and autosomal markers (51). These studies indicate that, after five centuries of colonization and cultural exchange, the modern populations of CA maintain most of their maternal Native American gene pool (~80%). Genetic distances based on mtDNA CR sequences assigned to the Native American haplogroups show that there is no significant genetic differentiation among CA populations, probably owing to their common origin and sustained gene flow (27).

Focusing on a larger geographical scale, haplogroup-level analyses suggested that present-day populations of CA have genetic affinities with populations from Patagonia (26), while studies of mtDNA CR sequences demonstrated similarities with the populations from West Central Argentina (WCA) (27,52,53). In contrast, morphological studies suggest biological relationships between CA and other populations from surrounding regions, which differ according to the variables and/or configurations of the samples used (54,55). Furthermore, according to the traditional model for the peopling of CA, supported by archeological evidence, the region should be closely linked to the Central Andean region (56–59). In addition to this model, there are different and opposing hypotheses regarding the origin of the CA populations: WCA (60,61), Patagonia (24,62) and Northeast Argentina (NEA) (20,21) have been suggested as the source areas. Nonetheless, the contrasting hypotheses concerning the origin and genetic affinities of CA populations have never been explored with high-resolution genetic markers such as whole mitochondrial genome sequences.

Several questions remain open about the demographic processes in CA within a more general context: From where did human populations originate? To what extent did this region maintain gene flow with populations from surrounding regions?

The aim of this study is to reconstruct the maternal history of the populations from CA through the description of the spatial and temporal patterns of molecular diversity within this region. Additionally, we provide relevant information about the general process of settlement of the Southern Cone of South America and describe new clades within the mitochondrial phylogenetic tree.

Results

Haplogroup characterization

To survey the mitogenome variation in CA, a total of 114 modern and ancient complete mitochondrial genomes (Fig. 1) were successfully sequenced with ~100% mitogenome coverage. The 92 modern and 22 archeological mitogenomes were sequenced to an average coverage depth of 1027× (120×–2615×) and 80.1× (14.3×–189×), respectively (Supplementary Material, Table S1). Damage profiles and estimates of mitochondrial contamination for the 22 archeological samples indicate authentic endogenous ancient DNA (Supplementary Material, Table S1b). We found no significant differences in the frequency distribution of haplogroups between modern and archeological samples (P = 0.229 ± 0.005) according to the exact test (63). From the total dataset, 19 individuals were assigned to haplogroup A2 (17%; 3 ancient and 16 modern), 25 to B2 (22%; 8 and 17), 25 to C1b (22%; 7 and 18), 5 to C1c (4%; 0 and 5), 10 to C1d (9%; 1 and 9) and 30 to D1 (26%, 3 and 27). The major haplogroups of Native American origin found in the modern mitogenomes show similar frequencies to a larger present-day sample set from CA (27), indicating that the sample analyzed here is largely representative (Supplementary Material, Table S2).

Population diversity statistics

Similar genetic diversity levels (Hd and π) were detected in the contemporary and ancient samples, although the former was slightly higher (Table 1). This result was consistent with the high diversity levels found with the mtDNA complete CR sequences in a separate study (27).

Table 1

Diversity parameters estimated from mitogenomes for 92 modern and 22 ancient samples of Central Argentina

Population	N	n	K	S	Hd ± SD	π ± SD	Tajima’s D	Fu’s F(s)
Modern mitogenomes	92	16 569	84	334	0.9980 ± 0.0020	0.0021 ± 0.0010	−1.4687 (P = 0.0411)	−24.0992 (P = 0.0007)
Ancient mitogenomes	22	16 564	20	140	0.9910 ± 0.0170	0.0020 ± 0.0010	−0.2014 (P= 0.4663)	−1.7345 (P= 0.1974)

N, sample size; n, number of observed nucleotide sites; K, number of different haplotypes; S, number of polymorphic (segregating) sites; Hd, gene diversity; π, nucleotide diversity; Tajima’s D, test of selective neutrality; Fu’s F(s), test of selective neutrality; SD, standard deviation. Values in italics indicate statistical significant difference (P<0.05).

Neutrality tests [Tajima’s D and Fu’s F(s)] did not reach significance for the ancient samples, whereas they did for the modern samples. Significant negative values were observed for both tests with modern sequences, indicating population growth, although the Tajima’s D was only marginally significant (Table 1). This is likely because Fu’s F(s) test is more sensitive to recent population expansions than Tajima’s D (64). To further explore this issue, we tested the hypothesis of neutrality in mtDNA variation from the modern set of samples by estimating both parameters but by separating the samples by haplogroup and partitioning the mitogenomes into CR, coding region and whole mitogenome components. From this analysis, we observed statistically significant values for most cases, mainly for the Fu’s F(s), further supporting population expansion. The exception was haplogroup D1, for which Tajima’s D test showed non-significant statistical values for the three partitions (Supplementary Material, Table S3).

We observed no significant genetic differentiation between the ancient and contemporary population samples: analysis of molecular variance’s (AMOVA’s) F_ST= 0.018 (P = 0.112) and Weir and Cockerham’s F_ST= 0.020 (P = 0.114). Since the ancient sample includes individuals from different dates and locations, and cannot be considered to be a real population, this result suggests that the population continuity in the region should be considered as only prospective.

Phylogenetic and phylogeographic diversity in CA

We compiled ~3700 sequences, including all available mitogenomes of Native American origin, both modern and ancient, to describe the phylogenetic relationships between the mitogenomes from CA and those from the rest of the Americas. Through the reconstruction of the maximum parsimony trees for each Pan-American haplogroup (A2, B2, C1b, C1c, C1d and D1; Supplementary Material, Figs S1–S6), we identified a large number (N = 85) of novel clades and subclades at different levels of branching. We also modified the definition of some clades (B2z, B2am, D1g5 and C1d1e), correcting the diagnostic motif previously established according to PhyloTree Build 17 (65) and other publications (see Supplementary section S1 and Supplementary Material, Figs S2 and S4–S6).

Each clade was classified as belonging to one of three categories—local, recent and relic—according to the phylogenetic and phylogeographic relationships identified within CA and between the other regions of the Americas. These are described in the three following subsections.

To facilitate the interpretation of the results, we define the geographic regions of South America to which we refer throughout this study. CA constitutes a transitional region between the eastern valleys of the Andean chains, here called as WCA, the Pampas-Patagonian lowlands in the east and south, and the plains of the Gran Chaco toward the north and northeast. The Salado River in Santiago del Estero province serves as the border between CA and the Gran Chaco, while there are no natural boundaries between CA and the rest of the surrounding areas. The Paraná and Paraguay rivers constitute the eastern geographical limit of the Gran Chaco with the subtropical forests of the Southern Cone, including the NEA region, while the 500 m.a.s.l. contour line to the west separates the Gran Chaco from the Andean region of Northwest Argentina (NWA).

Local clades

Local clades refer to those clusters that exhibit local differentiation and therefore encompass mainly CA samples. Phylogeographic and/or phylogenetic analyses revealed that some clades (A2bb, A2bc, B2ak, C1b16, C1d1b, D1g5a and D1j) had links with other regions of the Southern Cone, while for other local clades (A2ax, A2ay, A2ba, C1b₊₂₀₇, C1b17, C1c10, C1d1e and D1v) no such link could be identified. In what follows, we summarize the main results concerning these clades (additional details are provided in Supplementary section S2 and Supplementary Material, Figs S7 and S8). Phylogeographic analyses for local clades show that A2₊₁₅₀ (A2bb and A2bc), B2_+16 142 (B2ak1) and C1b₊₁₄₆ (C1b16) are specific to CA and WCA, although they are occasionally found in the Central Andes, Gran Chaco, Patagonia, Bolivia and Chile (Fig. 2). The remarkable correspondence between higher frequency and greater diversity is consistent with the hypothesis that these clades trace their origins to CA and WCA and subsequently spread from there.

Figure 2

Phylogenetic, phylogeographic and isofrequency analyses for CA local lineages. From left to right within (A–C), a summarized phylogenetic tree, the median joining network and the isoline map are represented. (A) Lineage A2₊₁₅₀ includes 77 samples, (B) lineage B2_+16 142 (B2ak1) includes 54 and (C) lineage C1b₊₁₄₆ (C1b16) includes 119. The phylogeographic analyses were performed using mtDNA sequences from positions 16 032–16 544 and 51–555. The black arrows show the nodal haplotypes. The dots in the isoline maps indicate the geographic location of the populations included in the analysis. Corresponding references and population grouping at the regional level are detailed in Supplementary Material, Tables S5 and S6.

The estimated times to the most recent common ancestor for the major haplogroups ranged from ~17 000 (A2) to ~12 000 (C1c) ybp, which are similar to those obtained in previous studies (66,67). Within them, the coalescence ages estimated for the local CA lineages varied enormously, with values ranging from ~11 000 ybp (D1j′w) to ~1200 ybp (D1v) (Fig. 3 and Supplementary Material, Table S4). The oldest dates are consistent with the archeological evidence that dates the peopling of CA back to at least 11 000 ybp (32,36) (and literature cited within).

Figure 3

Time to most recent common ancestor for founding haplogroups and CA local lineages. The coalescence ages are presented in Supplementary Material, Table S4. The TMRCA distributions were obtained across the sampled iterations from the BEAST procedure with the 937 mitogenomes and are represented with both gray violin plots and white boxplots. Braces: (1) internal branches for local clades, with estimated TMRCA consistent with a human presence since the Pleistocene–Holocene transition (~9000–11 000 ybp; upper shaded area), and (2) outermost nodes for local clades that occurred during the major pulse of population growth (~7200–4200 ybp; lower shaded area).

Recent clades

Recent clades refer to those clades that have low frequency in CA and are more frequent and diverse in other regions of South America. The presence of these clades in the region is most likely explained by recent migration processes because (i) none of the ancient CA mitogenomes fall into them and (ii) modern CA haplotypes belonging to them present few or no mutational steps from those described in other regions. In this study, we defined recent migrations as those demographic movements that have occurred since the 15th century, although there is no clear temporal break (68).

It is widely documented from the archeological and historical records that the Inca empire expansion first, and the European conquest later, caused demographic changes that affected the CA populations. The displacement of entire communities was a state policy of the Inca empire, employed to reaffirm their domination over newly annexed territories (69). Once the territory was occupied by the Spaniards, they established ‘Indian villages’, where Native people from different regions were grouped together (70). In addition to these coordinated forced relocations, many Indigenous groups voluntarily moved away from their original territories, avoiding Inca or European domination, and occupied unconquered neighboring regions such as the Gran Chaco lowlands (71). Recent migration processes also include all those movements of individuals or families that occurred in the 20th century, mainly to the Córdoba province, which is the second most populated province in Argentina and has many industries that attract workers from other regions (72).

In what follows, we summarize the main results concerning this category. The geographic distribution and genetic diversity observed for recent clades allow us to infer the sources of migrations into CA. First, the presence of A2bd, B2aj1a, B2z, D1f4 and D1u in CA is most likely explained by the recent migration from the Central Andes.

In addition, we identified a new clade called B2b2b, which has a wide distribution as it occurs at low frequency in modern populations from Amazonian Bolivia, NWA, Gran Chaco, Chile and CA, while it accounted for 7 out of 19 individuals from the archeological site of Pampa Grande, NWA (73). Since this clade occurred at very low frequency in relation to its wide geographical dispersion, it is very difficult to make a firm judgment about its specific origin. Within it, we define the sub-branch B2b2b1 that was found only in the CA and Gran Chaco mitogenomes, probably representing a rare and highly derived local lineage of B2b2b that likely spread along the Pacific route. Furthermore, clades A2az, B2b15, C1b18 and C1b19, which are marginally represented in the CA populations, also indicate a recent link between CA, Gran Chaco and NEA. In a wider geographical context, the clades B2e1 and the newly defined C1b15 and B2am show a genetic link between CA and Amazonia, Paraguay and NEA (Supplementary section S3).

Finally, it is not straightforward to interpret the presence of clades A2p and B2j in CA (Supplementary Material, Figs S1 and S2), as their known continuous distribution is restricted to Mexico, Central America, the Caribbean and Venezuela. Given their absence in the intermediate populations of Ecuador, Peru and Brazil, we propose that the CA sequences represent southern isolated occurrences that may be related to the late migrations in colonial times which were associated with movements of soldiers, merchants, priests, Crown officials, etc. with their families, servants and/or slaves. In Supplementary section S3, we detail other lineages included in the recent clades category.

Relic clades

CA relic clades are those that nest in the internal branches of previously described clades that span broad and diverse geographic regions but are separated by several mutational steps from them. A possible explanation for their presence in CA is that they are local survivors of clades that used to have wider geographic distributions. We identified two clades, C1b3a and C1c9a, that appeared at low frequency in CA but differed from the recent clades with discontinuous distribution by having well-differentiated haplotypes.

C1b3 has a disjointed geographical distribution in the United States (US) and South America (Supplementary Material, Fig. S3). The three US sequences belong to the contributors of Mexican or Hispanic origin (74–76), while the South American sequences come from admixed individuals from Argentina and Peru (77,78). A phylogeographical analysis of C1b_[phone] (C1b3a) (Supplementary Material, Fig. S9) shows a distribution in Argentina centered in CA, WCA and NWA, with sporadic cases in Gran Chaco and NEA, Northern Chile (79) and Lima, Peru (80).

The new clade C1c9 is differentiated into three sub-branches and has a broad geographical distribution in the American continent (Supplementary Material, Fig. S4). The first of the three sub-branches, C1c9a, is divided into two branches, one exclusive to CA and the other to the Gran Chaco, that represent the southernmost distribution of C1c9 lineages. The second clade, C1c9b, is restricted to Southern North America (75,81), while the third one is represented by a unique well-differentiated haplotype from Peru (77). The haplotype network for the southern branch (C1c_+16 169) does not exhibit any nodal haplotype but, instead, many unique haplotypes (Supplementary Material, Fig. S10A). The isofrequency map based on the C1c9a diagnostic polymorphism 16 169 shows a very low frequency for this clade, which is present in populations of the Gran Chaco, NWA and, to a lesser extent, in CA (Supplementary Material, Fig. S10B).

Single branches: allochthonous lineages or extinction?

We define single branches as those represented by one or two individuals who do not belong to any internal clade already defined in the phylogenetic trees. We observed a relatively low number of such single haplotypes that were not related to any others. This was the case for 7% of the total sample, including five modern and three ancient individuals (Supplementary Material, Figs S2–S5). Focusing on the SNPs in the CR, we searched for matches between CA and other South American populations and found none or only distantly related sequences (detailed information can be found in Supplementary section S4), so the presence of these single haplotypes in modern CA individuals is not consistent with a recent migration scenario. Therefore, we suggest that they are representatives of clades that underwent extinction (in the case of ancient samples), or that are in the way of being lost by genetic drift (in the case of modern samples), possibly affected by the reduction of the Indigenous population owing to epidemics, warfare and slavery related to European colonization (82).

Discussion

To date, most studies of complete mitochondrial genome sequences in the Americas have focused on specific haplogroups (7,10,67,83,84) and on understanding the human evolutionary history at broad geographical and temporal scales (10,66,67,85). Only a few publications have documented regional patterns of variation in order to address questions linked to a shorter time scale or interactions among closer populations (86–88).

In the present work, we focused on the sampling of a restricted geographical region, CA, but with a relatively large time transect. With this sampling strategy, we could define new clades through phylogenetic analyses of whole mtDNA genomes from the ancient and modern populations. We were also able to revise phylogenetic interpretations of previously characterized clades (see Supplementary section S1 and Supplementary Material, Figs S2–S5 and S8).

Based on our results, we observed, within each of the Pan-American haplogroups, a wide range of differentiation into derived branches, including 85 newly defined sub-haplogroups of which 15 are specific to CA. From the definition of these new clades and the analysis of mitochondrial diversity and geographic distribution, we formulated alternative hypotheses to explain the main evolutionary processes that led to the formation of the populations in CA.

Local clades as evidence of great temporal depth and population continuity in CA

Most of the local clades are defined here for the first time. They represent 64% of all the CA samples analyzed, and for most, the distribution is centered in CA and/or WCA, supporting the hypothesis of a close relationship between both regions (27,89). For several local clades, the coalescence ages were estimated to date between approximately 10 800 and 8500 ybp, revealing a great temporal depth (Fig. 3 and Supplementary Material, Table S4). These dates are consistent with the archeological evidence that indicates human presence in the region since the Pleistocene–Holocene transition (32,36).

Although some ancient samples analyzed here show private mutations that are absent in present-day populations, most of them (19/22) belong to clades that also include modern CA samples (Supplementary Material, Figs S1–S5). Moreover, most of those ancient haplotypes (17/22) nest inside branches of local clades having great temporal depth, suggesting long-term population continuity in the region. Based on radiocarbon dating, it is possible to establish a genetic link between modern people and individuals who lived at least 5000 ybp (the mean calibrated age of the oldest sample included in the clade B2ak1; sample #199, Supplementary Material, Table S1b). Further support for the continuity hypothesis includes both the frequency distribution of haplogroups and AMOVA’s F_ST results that show no significant differentiation between the modern and archeological samples. To discard a potential bias owing to the sample size difference between both groups, we also applied Weir and Cockerham’s F_ST index, and we report non-significant P-values estimated through random permutations.

To further investigate population continuity in CA, we constructed two mitogenomic trees, one of them using only the newly generated lineages (Fig. 4) and the other including these data together with an extensive comparative dataset (Supplementary Material, Fig. S11). The trees show that shared branches between ancient and modern haplotypes bifurcate along the stem of Native American clades, at times spanning over a wide range. Our results contrast markedly with those obtained by Llamas et al. (66), who sequenced 92 whole mitochondrial genomes from pre-Columbian South American inhabitants, mostly from Peruvian archeological sites. The phylogenetic analysis using their novel lineages and a comparative dataset of 206 modern Native American mitogenomes from both Indigenous and admixed populations showed that none of the ancient haplotypes shared a common ancestor with a modern haplotype more recently than ~9000 ybp, suggesting that ancient sequences represent unique branches from the stem of Native American haplogroups. Conversely, we find that the bifurcations between the ancient and modern samples in the mitogenomic trees span a large time range between 6500 and 500 ybp, much more recently than 9000 ybp (Fig. 4 and Supplementary Material, Fig. S11). The apparent contradiction between our results and those of Llamas et al. is likely attributable to different regional dynamics and/or the more comprehensive database used here, which provides a more representative spatial overlap between archeological sites and the present-day population locations. Analyzing the genetic diversity at a fine regional scale, leveraging modern and ancient samples that cover the same restricted geographical area proves to be essential to discriminate between extinction (irrecoverable loss of genetic diversity at some time) and the population continuity (i.e. stable genetic diversity across time) scenarios (90).

Haplogroup B2: extinction or recent gene flow?

Although haplogroup B2 represents only 17% of the Native American lineages in modern CA populations (27), it is much more frequent in archeological samples, reaching 42% in the mountain region of Córdoba province (44). The high diversity but low frequency observed in modern populations for B2 has been previously attributed to the extinction of ancient lineages (91). Under this evolutionary scenario, intermediate haplotypes would represent genetic ‘gaps’ that accentuate the difference between lineages in the phylogenetic tree (92). However, the phylogenetic analysis based on mitogenomes shows that 52% (13/25) of modern B2 samples nest in clades that are characteristic of other regions. Therefore, the high diversity found for the B2 haplogroup in modern CA populations is better explained by recent gene flow caused either by the expansion of the Inca empire, the impact of European conquest and/or by workforce arrivals from the surrounding areas during the 20th century.

Nonetheless, we also describe a new clade within B2, called B2ak, with important temporal depth (8945 ybp). This local clade exhibits great internal genetic diversity, and although it is highly represented among ancient samples (ranging between ~550 and ~5000 ybp), it is still found in present-day populations. Together, these patterns reveal that different demographic processes are involved in the evolutionary history of the B2 haplogroup. This example demonstrates the importance of considering the internal branches in mtDNA studies, beyond studies at the haplogroup level.

Although the evidence does not support the hypothesis of massive extinction of lineages owing to the impact of European colonization, it is indisputable that the European conquest led to a series of transformations that affected the entire ‘Indigenous world’ by which most of these societies were disrupted and eventually ‘disappeared’ (93). However, other groups survived and persisted, some of them reduced in ‘Indian villages’ within the system of Hispanic domination (46), and many other individuals, especially women, were incorporated into the colonial society.

Genetic relationships between neighboring regions

We performed phylogenetic and phylogeographic analyses with the modern and ancient mitogenomes, including an extensive database of mitogenomes and mtDNA CR sequences from the Americas, to analyze genetic links that reflect the interactions of CA populations with other regions in South America.

Relationships with WCA

Despite the lack of published mitogenomes from WCA, the phylogeographic analysis based on CR sequences shows that most of the CA local clades are also highly represented in the modern samples from WCA, indicating a possible common origin and/or sustained gene flow between those regions as supported by morphological (42,55), linguistic (94) and ethno-historical (60,61) evidence. Further analysis of WCA mitogenomes will be needed to understand the nature of the links between both regions.

Relationships with Central Andes and NWA

In total, 7% of the CA samples belong to Central Andes-specific clades. These are generally haplotypes that are separated by just a few mutational steps from those published for samples from Peru, Bolivia, Chile and NWA, indicating a recent link with these regions. The presence of Andean haplotypes probably reflects gene flow associated with the Inca expansion or post-conquest processes. It is worth mentioning that B2aj, a characteristic haplotype of the Andean Plateau populations (95) that spread successfully into NWA from Bolivia (10), was only found in a single modern individual from the Córdoba province. This finding, together with the fact that most of ancient B2 samples nest into local clades, allows us to review the proposal of Nores et al. (44) that the high frequency of B2 found in past populations after 1200 ybp from the Córdoba province was owing to significant contacts with the populations from the Andean region.

Relationships with Gran Chaco and NEA

As mentioned before, some clades marginally represented among the CA samples are also present in the Gran Chaco and NEA. In all cases, they share similar haplotypes, indicating a recent contact between these regions, with gene flow probably in both directions. For example, the presence of the lineage C1d1b_[phone] in archeological samples from Santiago del Estero (96) and Córdoba (sample #122, dated to ~2300 ybp) provinces as well as in modern inhabitants of Gran Chaco suggests gene flow between both regions (96). The archeological evidence suggests movements of Gran Chaco groups into the south of Santiago del Estero province by at least 1000 ybp (97,98). For the Mar Chiquita Lagoon (northeast of the Córdoba province), contact between these regions since at least 2000 ybp is supported by the ethnohistorical (99), bioanthropological (100) and archeological evidence (41,45,101,102). Ethnohistorical evidence further supports the relationship with Gran Chaco as it describe post-contact south-to-north movements of Indigenous people who took refuge in the Gran Chaco from Spanish conquerors (71), and as a result of 19^th century movements to occupy new territories for breeding (103).

Relationships with the Pampas

As only two mitogenomes have been published for the Pampean region of Argentina, it is therefore difficult to decipher the genetic links of CA with this region. On the one hand, these two sequences, reported from the Arroyo Seco 2 site and dated to more than 7000 ybp (66), are not related to CA sequences. On the other hand, HVRI sequences from an archeological site in the Pampas–Patagonia transition revealed the presence of one individual (out of 13) carrying a D1j haplotype, thus suggesting a possible relation with CA (104). A recent study of Native American ancestry with a large set of autosomal markers assigned admixed individuals from both Central and Pampas regions to the same cluster, which also includes individuals belonging to Native American communities from the Gran Chaco, the Amazonia and Northern Andes (105). Thus, more studies of ancient samples from the Pampas will be needed to elucidate the nature and extent of the genetic links and to better understand the biological and cultural relatedness between this region and CA, suggested by craniometric analysis (106) and shared burial traditions (107).

Relationships with Patagonia

In this research, we have found no evidence linking CA to Patagonia. While a single CA sample was grouped with others from Northern Patagonia, defining the sub-haplogroup C1d1f, this sample represents an example of the migratory process to CA during the last decades according to the genealogical information. Another clade that deserves discussion is D1g5, a lineage found in two ancient and three modern samples of CA. Despite D1g being a clearly Patagonian clade (8), D1g5 is the most widespread subclade within D1g, ranging from Tierra del Fuego to the extreme north of Argentina (108). Moreover, only D1g5a, which accounts for four out of the five D1g5 CA samples, is present at lower latitudes, while the rest of the subclades within D1g and D1g5 are found almost exclusively in the central–southern part of Chile and Argentina (8,108). Thus, the presence of D1g5a in CA is unlikely to reflect a link with Patagonia. In addition, clades C1b13, B2i2 and D4h3a5, which are highly represented in Patagonia (7,8), have not been found among the CA samples.

Overall, our results do not support previous claims pointing to a possible link between the CA and the Patagonian region based on archeological (38) and morphological evidence (21,24,109) as well as mtDNA haplogroup distribution (26,44). The apparent genetic similarity is owing to the high frequency of C and D haplogroups in both regions; however, after increasing the resolution with complete mitochondrial sequences, we show that these regions do not share the same sub-haplogroups.

Demographic model for CA: a genetic variation mosaic

The results obtained in the present work reveal a mosaic of genetic diversity in CA with several local clades, most of them having great temporal depth, and demonstrate that a dense sampling scheme of a particular region allows to reach sufficient resolution of the genetic diversity to better understand the population dynamics. In this sense, it is possible to identify the genetic signals of population movements, either leading to gene flow among populations that occurred in the last centuries at different intensities and that have changed the mitochondrial variability found in modern CA populations, or to partial population replacements that occurred in more remote times.

Although it is not possible to know exactly the processes that led to the presence of foreign clades in CA, the archeological and historical evidence allows us to interpret the presence of Andean lineages as a result of major forced population movements into CA, which were caused by Inca and/or European conquest. In the case of clades shared with the Gran Chaco and NEA regions, they are probably the result of individual or kin-structured migration with these regions. Also, some isolated haplotypes could have arrived in the region as a consequence of individual or family migrations in the historic or modern times. Further investigation would be needed to evaluate the alternative hypotheses.

Based on our findings, we propose a demographic model for the population of CA within the general process of peopling of the Southern Cone of South America (Fig. 5). According to this demographic model, CA is considered as the center of origin and dispersal of many lineages, some of which were probably lost through genetic drift or extinction, but most having survived, revealing the maternal genetic continuity between ancient and modern populations.

In this work, we found that local clades from CA are not nested inside any other South American clade with a wide spatial distribution. As a consequence, it is not possible to discern whether the ancient populations from CA arrived from the Pacific or Atlantic routes, or both. Whatever their entry route was, the first inhabitants reached CA during the Pleistocene–Early Holocene transition as evidenced by the archeological record and our TMRCA estimates for local clades. In addition, the presence of some lineages with discontinuous distributions (C1b3 and C1c9) suggests the genetic persistence of relic groups that arrived in remote times, probably with the first founders. Following early settlements, significant and distinctive mitochondrial diversity emerged over millennia. By 8000 ybp, the archeological record reports a greater occurrence of residential sites in several regions of CA (38,39,110–112), which is consistent with the population growth during that time. Furthermore, lithic typology changed around 7500 years ago, which was associated with the appearance of triangular points. This new technology was established around 5000 ybp and reached a wider territorial expansion that continued until around 2000 ybp (37). The genetic data we present here provide further support to the archeological evidence of population size increase. Neutrality tests [Tajima’s D and Fu’s F(s)] hold significant negative values for modern mitogenomes, indicating population expansion, while the Bayesian skyline plot (BSP; Fig. 6) shows a rather sharp increase in the population size between 7500 and 5000 ybp. In agreement, the TMRCA for the outermost nodes of the local clades (A2bb, A2bc, B2ak1 and B2ak1a) also suggests a population growth between 7200 and 4200 ybp (Fig. 3 and Supplementary Material, Table S4). Finally, BSP analysis shows that a smooth demographic contraction occurred immediately after the population growth and lasted up to ~2500 ybp. This was followed by a long period of demographic stability until the present time (Fig. 6). Finally, for recent centuries, our results point to gene flow from the surrounding regions, mainly from Central Andes and, to a lesser extent, from Gran Chaco, NEA and the subtropical forests of South America.

Figure 4

Phylogenetic relationships between modern and ancient mitogenomes. Dendrogram that includes all the mitogenomes from the present study, with ancient samples shown as bold lines. The stars highlight the most recent coalescence of an ancient sample with a modern one.

In conclusion, given the deep temporality of the local CA clades and the high levels of molecular diversity, which do not differ between the modern and ancient mitogenomes, the central region of Argentina should no longer be seen as a mere ‘crossroads’ where the Andean, Patagonian, Gran Chaco and NEA components met. On the contrary, the populations of CA reveal themselves as having their own identity in genetic terms.

Figure 5

Demographic scenario for the peopling and settlement of CA. The arrows represent the relationships and migration routes of the clades discussed in the article.

Materials and Methods

Modern samples

Complete mtDNA sequences were obtained from 92 unrelated individuals randomly chosen among those carrying Native American lineages from a set of 812 biological samples collected in CA (27). Specifically, this modern dataset includes individuals from the provinces of Córdoba (N = 52), Santiago del Estero (N = 16) and San Luis (N = 24) (Fig. 1). All participants gave written consent to be part of this non-profit scientific investigation. Information gathered in the field regarding birthplaces of the subjects and their parents indicates that the inhabitants of the studied populations present low-to-moderate mobility and are restricted to their places of origin by at least two or three generations in most of the cases. Only 7 of the 92 individuals declared to have a maternal ancestor who came from a different region than CA (see Supplementary Material, Table S1a and Supplementary section S5). All associated information concerning modern samples is provided in Supplementary Material, Table S1a. Laboratory methods for analyses of modern samples (DNA extraction, preparation of Illumina sequencing libraries, massively parallel sequencing and raw data analysis), are detailed in Supplementary section S6 and Supplementary Material, Table S7.

Figure 6

Bayesian skyline plot of female effective population size based on a generation time of 25 years. We obtained a population size estimate over time for each sampled iteration from the Bayesian skyline analysis using the 113 CA mitogenomes. Black line and shaded area represent the average and 95% credibility interval of these estimates over time, respectively. Vertical dashed lines represent the time range for changes reported in the archeological record, which are likely linked to the major pulse of population growth (see main text).

Archeological samples

DNA extraction was performed from the bone or dental remains of 22 individuals of archeological origin, which have been preserved in the Museo de Antropología (Facultad de Filosofía y Humanidades, Universidad Nacional de Córdoba), and in the different public and private museums in the Córdoba province. The samples were originally obtained from different archeological sites located in the Córdoba province (Fig. 1). Macroscopic analysis of the archeological remains showed a general appearance of good preservation, with very little bone loss and postmortem fragmentation (25). Radiocarbon dating of 20 of these individuals span ~4000 years of local history, from 459 ± 40 to 4450 ± 80 uncalibrated ¹⁴C ybp. All associated information concerning archeological samples (including site location, current collection location and specimen description) is provided in Supplementary Material, Table S1b. Laboratory methods of the ancient samples (DNA extraction, preparation of Illumina sequencing libraries, mtDNA hybridization capture enrichment, massively parallel sequencing and raw data analysis) are detailed in Supplementary section S6.

Ethical statement

This study was performed following ethical guidelines for working with human remains, treating these deceased people with due respect. Our research program with modern populations and ancient human remains was approved by the Ethics Committee of the CEMIC (Comité de Ética en Investigación, Centro de Educación Médica e Investigaciones Clínicas ‘Norberto Quirno’). Skeletal samples were exported with proper permits from the Córdoba province and the Argentina government.

Data edition and Haplogroup assignment

Polymorphic sites were defined by differences from the rCRS (113) using mitoSAVE. This software converts data within the VCF into mtDNA haplotypes using phylogenetically established nomenclature for the human mitochondrial genome. Each position was only considered if there was a 10× minimum coverage. The data were verified manually by checking BAM files in IGV (Integrative Genomics Viewer) 2.3.34 (114). Once the debug and edited files were generated, .hsd files were uploaded to HaploGrep 2 (115) to assign mitochondrial haplogroups. Finally, all haplogroup assignments were manually confirmed by checking diagnostic positions as described in PhyloTree Build 17 (65).

Modern sample haplogroups obtained through NGS were verified by comparing with the previously Sanger-generated CR sequences. The data edition also implicated the exclusion for all calculations of the transitions at the hypervariable position 16 519 and any length variations in the C-stretch between nucleotides 16 181–16 193, 303–315 and 513–524. The mitogenomes in fasta format were aligned using the R package DECIPHER (116,117), sequentially calling functions AlignSeqs and AdjustAlignment. The alignments were checked manually in HaploGrep 2 (115).

Genetic diversity

Different measures of genetic diversity (sequence number, polymorphic sites, sequence diversity and nucleotide diversity), AMOVA, pairwise genetic distances (Φst and pairwise difference) and exact test of population differentiation were performed using Arlequin 3.5.1.3 (118). We also computed Weir and Cockerham’s F_ST (119) with the pp.fst function from the R package hierfstat, and its P-value was estimated randomly by permuting sample groups for 10 000 times.

Tests of neutral selection such as Tajima’s D (64) and Fu’s F(s) (120) were run using a bioinformatic pipeline (121). Statistical significance was estimated either from the expected distribution (for Tajima’s D) or from 10 000 simulations [for Tajima’s D and Fu’s F(s)]. Simulations were run using the three following steps: (i) generate 10 000 gene tree replicates with ms (122), containing the same number of sequences as the tested sample and with mutation parameter (4N₀μ) equal to its Watterson estimator (θ_W) (123); (ii) for each gene tree replicate, generate the sequences with Seq-Gen (124), following the HKY model with the number of nucleotides for the genomic region analyzed as the total sequence length l and the value θ_W/l as the branch length scale parameter and (iii) for each sequence sample thus obtained, calculate Tajima’s D and Fu’s F(s) scores as for the real data. P-values were then calculated as the number of simulations exhibiting a score lower than that observed for real data.

Phylogenetic analysis

We built maximum parsimony trees, one for each of the haplogroups A2, B2, C1b, C1c, C1d and D1 using the mtPhyl software (https://sites.google.com/site/mtphyl/home), which were hand-modified with reference to PhyloTree Build 17 (65). In total, 278 mitogenomes were included in the analysis (164 published and 114 from this study). The database is detailed on Supplementary Material, Table S5. Novel clades and subclades were defined when they included a minimum of two different sequences and shared at least one stable mutation. Mitogenomes used to build the phylogenetic trees were selected from an in-house database containing ~3700 sequences, which includes all the publicly available mitogenomes of Native American origin, both modern and ancient. The use of a large dataset allowed exploring the phylogenetic affinities between our samples and ancient and modern mitogenomes of the Americas. Because the capital cities of South American countries, such as Buenos Aires and Lima, disproportionately concentrate a large part of the nation’s populations owing to migrations of the last century (125), influencing inevitably the current lineages’ distribution (126), they were not considered when discussing the distribution/origin of lineages.

Demographic analysis

We performed a Bayesian skyline analysis using BEAST v1.8.3 to reconstruct past population dynamics. BEAST runs were performed only with the 113 full mtDNA genomes from CA (we excluded an undated ancient sample). For these analyses, we used GTR nucleotide substitution model, the best substitution model estimated with MrMTgui v1.0, and applied the corresponding substitution rates reported previously in Soares et al. (127). We assumed that every branch in a phylogenetic tree evolves according to the same evolutionary rate and therefore we used a strict clock. BEAST was then run with 100 000 000 iterations, sampling each 10 000^th iteration. We thus obtained a population size estimate over time for each sampled iteration. Finally, we reconstructed the consensus phylogenetic tree with TreeAnnotator, using default parameters, and plotted it using TreeStat v.18.3.

Haplogroup and sub-haplogroup coalescence age estimation

Coalescence times were estimated using a Bayesian approach. In addition to the 113 novel modern and ancient Native American mitogenomes from CA populations, we included a representative sample (N = 824) of the modern and ancient Native American mitogenomes from South America from the following publications: Fagundes et al. (85) (N = 47), Bodner et al. (5) (N = 45), de Saint Pierre et al. (7) (N = 46), Llamas et al. (66) (N = 92), Barbieri et al. (9) (N = 116), Brandini et al. (67) (N = 223), Sans et al. (128) (N = 5), Sevini et al. (78) (N = 124) and Simão et al. (88) (N = 90). Some of these samples (N = 8) were excluded from the analysis owing to an uncertain assignment in a phylogenetic tree which was performed to quickly test the usefulness of each downloaded fasta file. We additionally included all the published mitogenomes contained in the branches of interest for which the coalescence age was estimated, and these also were used for phylogenetic analysis (N = 44).

This set of mitogenomes was chosen so that the Bayesian analyses were computationally feasible while the dataset remained representative for both genetic diversity and geographical range. BEAST was then run on this 937 sequences alignment with 300 000 000 iterations, sampling each 10 000^th iteration, HKY evolution model, and by applying the corresponding substitution rates reported previously in Soares et al. (127) and by using a strict clock. Using a burn-in of 30 000 000 iterations, we then obtained the TMRCA with Tracer and reconstructed the consensus phylogenetic tree with TreeAnnotator with default parameters. The tree was finally plotted with iTOL online tool (129).

Phylogeographic analysis

For comparative purposes, we included in the analysis 7042 mtDNA CR published sequences from other populations of South America (Supplementary Material, Table S6a). Narrowing the mtDNA range to the CR allows us to build a more extensive database and to analyze the phylogeographic patterns with accuracy. CR sequences between rCRS np 16 032–16 544 and 51–555 were used. Polymorphic sites 16 182C, 16 183C, and 16 519 as well as indels’ variations at 16 181–16 193, 303–315 and 513–524 were not considered owing to their high mutation rates (127).

Haplotype networks

For selected sub-haplogroups (A2₊₁₅₀, B2_+16 142, C1b₊₁₄₆, C1b_[phone], C1c_+16 169, C1d1₊₁₉₅ and D1_[phone]), median joining networks were estimated with Network 4.6.1.3 (Fluxus Technology). Maximum parsimony was used as the post-processing option, using the three criteria proposed by Crandall and Templeton (130) based on predictions of the coalescence theory: frequency, topology and geography. The database used for median joining networks is detailed in Supplementary Material, Table S6b.

Frequency maps

Based on the frequency of these selected sub-haplogroups, we built isoline maps, following the Kriging procedure using Surfer 7 (Golden Software Inc., Golden, CO). Isolines start at a frequency of 1%, and the lower frequencies are not drawn on the map. The populations considered are detailed in Supplementary Material, Table S6c, and the frequencies of each sub-haplogroup for each population are given in Supplementary Material, Table S6d.

Data availability statement

The mitogenomes generated during this study are available in GenBank database (accession numbers MW886102 - MW886215, Supplementary Material, Table S1).

Acknowledgements

We especially thank the participants for their cooperation with our work during the collection of biological samples. We are grateful to the Museo de Antropología de Córdoba (FFyH, UNC), Museo Histórico Municipal (La Para, Córdoba), Museo de Ciencias Naturales de la Región de Ansenuza ‘Aníbal Montes’ (Miramar, Córdoba) and Museo Municipal ‘Capitán Juan de Zevallos’ (Valle Hermoso, Córdoba) for allowing us to access their collections. We thank Dr Minoru Yoneda and Dr Mai Takigami for radiocarbon dating at the Graduate School of Frontier Sciences, University of Tokyo, Japan. We thank Juan B. Berros and Hernán Dopazo for giving us unlimited access to the Biocodices S.L. computational cluster to perform the BEAST analyses. We also thank Ulises D’Andrea and Beatriz Nores who generously provided the ancient sample #199 from Alpa Corral, Córdoba.

Conflict of Interest statement. We declare no conflict of interest.

Funding

This work was supported by Agencia Nacional de Promoción de la Investigación, el Desarrollo Tecnológico y la Innovación (PICT 2007-1549, PICT 2012-711 and PICT 2015-3155), Secretaría de Ciencia y Tecnología (Universidad Nacional de Córdoba), Ministerio de Ciencia y Tecnología de la Provincia de Córdoba (PID 2018-79) and Consejo Nacional de Investigaciones Científicas y Técnicas (2015-11220150100953CO). M.P. is a postdoctoral fellow and A.G., R.N., J.M.B.M, C.M.B., M.F. and D.A.D. are research career members of CONICET, Argentina.

Author notes

Angelina García, Rodrigo Nores, Elizabeth A. Matisoo-Smith, David Comas and Darío A. Demarchi authors contributed equally to this work.