Over the last few decades, bioanthropology and archaeological research has incorporated the use of genetic information recovered from skeletal remains to examine a variety of issues regarding the origin and worldwide dispersion of our species. DNA studies on past and current populations have been at the center of the latest models of occupation of the Americas and of the different microevolutionary processes shaping the genetic diversity of native populations (e.g., Brandini et al. Reference Brandini, Bergamaschi, Cerna, Gandini, Bastaroli, Bertolini and Cereda2017; Lindo et al. Reference Lindo, Haas, Hofman, Apata, Moraga, Verdugo and Watson2018; Perez et al. Reference Perez, Postillone, Rindel, Gobbo, Gonzalez and Bernal2016; Posth et al. Reference Posth, Nakatsuka, Lazaridis, Skoglund, Mallick, Lamnidis and Rohland2018). The use of massively parallel sequencing has facilitated the recovery of mitochondrial and nuclear genomes, enabling increasing precision in the definition of American founding mitochondrial lineages, the description of regional peopling processes and interpopulation dynamics, the identification of adaptations to local environments, and changes in demographic patterns (De la Fuente et al. Reference De la Fuente, Ávila-Arcos, Galimany, Carpenter, Homburger, Blanco and Contreras2018; Fehren-Schmitz and Georges Reference Fehren-Schmitz and Georges2016; Lindo et al. Reference Lindo, Haas, Hofman, Apata, Moraga, Verdugo and Watson2018; Posth et al. Reference Posth, Nakatsuka, Lazaridis, Skoglund, Mallick, Lamnidis and Rohland2018).
However, the application of such procedures to archaeological samples poses challenges and requires strategies that cannot be readily implemented by most researchers because of their high cost or the lack of adequate infrastructure, equipment, and specialized technical staff; one such difficult-to-implement strategy is the fabrication of special probes to capture the target DNA and the enrichment of sequencing libraries. Therefore, analyses on South American archaeological samples are usually conducted in the context of cooperation agreements between local researchers and nonlocal research centers mostly located in the United States and Europe (Wade Reference Wade2018). As a result, the involvement of South American researchers in ancient DNA (aDNA) studies is frequently subordinated to the research agenda of the nonlocal laboratories. This situation discourages the pursuit of locally relevant research objectives that often involve more restricted time or spatial scales (Arencibia et al. Reference Arencibia, Crespo, Guraieb, Russo, Dejean and Goñi2019; Cardozo et al. Reference Cardozo, Dejean, Russo, Mazza, Loponte, Acosta, Feuillet and Corne2021; Crespo et al. Reference Crespo, Cardozo, Tessone, Vázquez, Kisielinski, Arencibia, Tackney, Francisco Zangrando and Dejean2020; Motti et al. Reference Motti, Winingear, Valenzuela, Nieves-Colón, Harkins, García Laborde, Bravi, Guichón and Stone2020; Postillone et al. Reference Postillone, Martínez, Flensborg and Dejean2020; Russo et al. Reference Russo, Mendisco, Avena, Dejean and Seldes2016, Reference Russo, Dejean, Avena, Seldes and Ramundo2018; Tavella et al. Reference Tavella, Demarchi and Nores2020). The development of such asymmetrical relationships is common to many scientific fields that depend on costly technological development; this “subordinate integration” (Kreimer Reference Kreimer2006; Kreimer and Ugartemendía Reference Kreimer and Ugartemendía2007) is a recurring phenomenon in the global division of scientific labor today. In this context, the procurement of resources to support the pursuit of avenues of research that are consistent with local and regional goals is essential to the local advancement of archaeological science.
The Beagle Channel area holds a particular interest for genetic studies. The area was probably one of the last territories to be reached by humans, which is one reason why it has been intensively studied over many years by different research groups. These studies have yielded a large and diverse body of information that has provided valuable contextual information for framing genetic studies. The earliest archaeological signal in the region dates to about 8000 years BP (Zangrando et al. Reference Zangrando, Bjerck, Piana, Breivik, Tivoli and Negre2018, Reference Zangrando, Tivoli, Ponce, Alunni, Fernández Ropero and Martinoli2022) and the available zooarchaeological, technological, and isotopic evidence points to the existence of maritime adaptations in the area as early as 6400 years BP (Orquera and Piana Reference Orquera and Piana1999, Reference Orquera and Piana2009; Suby et al. Reference Suby, Zangrando and Piana2011; Zangrando Reference Zangrando2009). Additionally, Patagonia in general has exceptional preservation conditions that facilitate the recovery of genetic material, a fact made evident by the increasing number of Patagonian aDNA analyses performed in laboratories around the world (Arencibia et al. Reference Arencibia, Crespo, Guraieb, Russo, Dejean and Goñi2019; Crespo et al. Reference Crespo, Dubois, Gabriela Russo, Lanata and Dejean2017a, Reference Crespo, María Russo, Lanata and Dejean2017b, Reference Crespo, Cardozo, Tessone, Vázquez, Kisielinski, Arencibia, Tackney, Francisco Zangrando and Dejean2020; de la Fuente et al. Reference De la Fuente, Galimany, Kemp, Judd, Reyes and Moraga2015, Reference De la Fuente, Ávila-Arcos, Galimany, Carpenter, Homburger, Blanco and Contreras2018; Lalueza et al. Reference Lalueza, Pérez-Pérez, Prats, Cornudella and Turbón1997; Motti et al. Reference Motti, Hagelberg, Lindo, Malhi, Bravi and Guichón2015, Reference Motti, Sebastián Muñoz, Cruz, D´Angelo del Campo, Borrero, Bravi, Guichón, Otero, Svoboda and Banegas2019; Nakatsuka et al. Reference Nakatsuka, Luisi, Motti, Salemme, Santiago, D'Angelo del Campo and Vecchi2020; Parolin et al. Reference Parolin, Galimany, Otero, Dahinten, Gabriela Millán, Moraga, Otero, Svoboda and Banegas2019; Raghavan et al. Reference Raghavan, Steinrücken, Harris, Schiffels, Rasmussen, DeGiorgio and Albrechtsen2015).
From a genetic point of view, maternal lineage C1 has been found in high frequencies in the Beagle Channel area (Tierra del Fuego, Argentina; Figure 1) and D1g and D4h3a have also been identified, albeit in lower proportions (Crespo et al. Reference Crespo, María Russo, Lanata and Dejean2017b, Reference Crespo, Lanata, Cardozo, Avena and Dejean2018; de la Fuente et al. Reference De la Fuente, Galimany, Kemp, Judd, Reyes and Moraga2015, Reference De la Fuente, Ávila-Arcos, Galimany, Carpenter, Homburger, Blanco and Contreras2018; Lalueza Fox Reference Lalueza Fox1996). Similar frequencies of the C1 and D1g lineages have also been observed in ancient samples from nearby areas, such as Peninsula Mitre and Isla de los Estados; these findings could indicate different degrees of contact or shared origins among diverse native groups in Tierra del Fuego and the Beagle Channel (Crespo et al. Reference Crespo, María Russo, Lanata and Dejean2017b, Reference Crespo, Cardozo, Tessone, Vázquez, Kisielinski, Arencibia, Tackney, Francisco Zangrando and Dejean2020; Nakatsuka et al. Reference Nakatsuka, Luisi, Motti, Salemme, Santiago, D'Angelo del Campo and Vecchi2020).
Region of interest. The star indicates the location of the Paiashauaia 1 burial by the Beagle Channel (Tierra del Fuego, Argentina); the circle indicates the location of the IPY08 burial by Seno Año Nuevo (Chile; de la Fuente et al. Reference De la Fuente, Ávila-Arcos, Galimany, Carpenter, Homburger, Blanco and Contreras2018). Image generated with Google Earth (v7.3.3.7786).
Given the region's distinctive characteristics, we studied an individual recovered from an archaeological site located on the Beagle Channel coast by analyzing its lineage within the regional phylogenetic context. To this end, we implemented a process to build and enrich an aDNA library and to undertake the massively parallel sequencing methodologies and bioinformatic procedures that are necessary to obtain a mitogenome “made in Argentina” (i.e., only using locally available resources).
Materials and Methods
The Analyzed Individual
An upper left incisor was taken from human remains recovered from Paiashauaia I, a primary mortuary context located at a rocky outcrop next to the mouth of the Paraná River, a few meters from the Beagle Channel coast (Figure 1; laboratory code: PZ3). Bioarchaeological analysis determined that the individual was a female between 35 and 45 years old, dated at 1504 ± 46 years BP by radiocarbon AMS (Suby et al. Reference Suby, Zangrando and Piana2011; for more details on methodology, see Supplemental Text 1). These remains are currently preserved in the Museo del Fin del Mundo's repository in Tierra del Fuego; permits for sample transportation and analysis were obtained from the Secretaría de Cultura of Tierra del Fuego in 2016. In addition, the sample was treated with respect and scientific rigor as set forth by the code of ethics for the study of human remains of the Asociación de Antropología Biológica of Argentina (2007; Aranda et al. Reference Aranda, Barrientos and Del Papa2014).
The hypervariable regions 1 and 2 of PZ3's mitochondrial DNA were successfully recovered, amplified, and sequenced by Sanger (Crespo et al. Reference Crespo, Cardozo, Tessone, Vázquez, Kisielinski, Arencibia, Tackney, Francisco Zangrando and Dejean2020). According to this study, PZ3's mitochondrial lineage was a distinctive D1g exclusively found in the Beagle Channel area.
Laboratory Protocols
All precautions for preventing and detecting exogenous contamination in the ancient sample were followed (see the detailed protocol in Supplemental Text 1). DNA extraction and library preparation were carried out at the aDNA laboratory clean room facility at the Universidad Maimónides in Buenos Aires. Library amplification, bait preparation, and capture were done at the laboratories of the Instituto Nacional de Tecnología Agropecuaria (INTA). The sample was powdered for DNA extraction. Extraction blank controls were conducted and subjected to PCR for the hypervariable region, along with the archaeological samples, to monitor for possible contamination.
The aDNA library was constructed using the NEBNEX Ultra II DNA Library Prep Kit for Illumina following the manufacturer's instructions (Supplemental Text 1). DNA fragmentation was not performed before construction of the library. After USER enzyme (NEB) treatment, the adaptor-ligated DNA was purified without size selection using AMPureXP Beads (Beckman Coulter).
The indexed aDNA library was captured in solution using 3′ biotin-labeled probes with the mitochondrial genome of an H haplogroup blood donor, coupled to Dynabeads® M-280 Streptavidin Magnetic Beads (Invitrogen; Maricic et al. Reference Maricic, Whitten and Pääbo2010); the solution was incubated at 55°C for 48 hours. After washing, the captured library was eluted by heating, and the DNA was purified (Supplemental Text 1).
The captured aDNA library was sequenced by synthesis with MiSeq Reagent Kit v2 (2 × 150) at INTA's Consorcio Argentino de Tecnología Genómica (CATG). PCR and Sanger sequencing were used to confirm SNPs in two low-depth regions (primers specified in Supplemental Table 1).
Bioinformatic Analysis
Processing and Mapping of Sequenced Reads and Variant Calling. The quality of the Illumina reads was evaluated with FastQC (v0.10.1), and the reads were trimmed and end-clipped to a PHRED score of 30 using Cutadapt (v1.16; Martin Reference Martin2011; see the detailed protocol in the Supplemental Text 1). Paired-end reads were then merged using AdapterRemoval (v2.3.3; Schubert et al. Reference Schubert, Lindgreen and Orlando2016) and mapped to the revised Cambridge Reference Sequence (rCRS; Andrews et al. Reference Andrews, Kubacka, Chinnery, Lightowlers, Turnbull and Howell1999) using the Burrows Wheeler Aligner (BWA; v0.7.10; Li and Durbin Reference Li and Durbin2009). The output was converted to a BAM file using Samtools (v1.7; Li and Durbin Reference Li and Durbin2009), and one 3′ nucleotide was trimmed in the BAM file with bamUtil (v1.0.15; Jun et al. Reference Jun, Wing, Abecasis and Kang2015). Mapping statistics were calculated with Bamtools (v2.5.1; Barnett et al. Reference Barnett, Garrison, Quinlan, Str̈mberg and Marth2011) and Samtools, duplications were marked and discarded using Picard-tools (v2.18; http://broadinstitute.github.io/picard/), and aDNA damage patterns were assessed using MapDamage2 (v2.0; Jónsson et al. Reference Jónsson, Ginolhac, Schubert, Johnson and Orlando2013).
Variant calling was performed using BCFtools (v1.9; Danecek and McCarthy Reference Danecek and McCarthy2017; Narasimhan et al. Reference Narasimhan, Petr Danecek, Yali Xue and Durbin2016). The information yielded by the Sanger sequences and obtained for HVR1 in Crespo and colleagues (Reference Crespo, Cardozo, Tessone, Vázquez, Kisielinski, Arencibia, Tackney, Francisco Zangrando and Dejean2020) was incorporated into variant calling. Finally, variant filtering was also made using BCFtools.
Construction of the Consensus Mitogenome. A final VCF file was used to construct a consensus genome of the individual PZ3 using BCFtools (Supplemental Text 1). Positions with depth <2 were treated as missing data (N). The mitochondrial haplogroup was determined by using HaploGrep2 (v2.2; Weissensteiner et al. Reference Weissensteiner Hansi, Kloss-Brandstätter, Forer, Specht, Bandelt, Kronenberg, Salas and Schönherr2016).
Haplotype Polymorphisms and Phylogenetic Relationships
We compiled a database that included ancient and modern D1 sequences from Argentina, Brazil, and Chile to compare the PZ3 mitochondrial haplotype to other D1 mitogenomes. Sequences reported up to 2016 were obtained from Phylotree Build 17 (February 18, 2016; van Oven and Kayser Reference Van Oven and Kayser2009), and more recently reported sequences were retrieved from the literature. One representative of each of the D1 sublineages not reported in the three countries was also included for analysis, as well as one D4h3a and the RSRS (Behar et al. Reference Behar, Van Oven, Rosset, Metspalu, Loogväli, Silva, Kivisild, Torroni and Villems2012), yielding a final data matrix of 79 sequences (Supplemental Table 2).
Sequences were aligned using the FFT-NS-2 algorithm in MAFFT (v7.471; Katoh et al. Reference Katoh, Misawa, Kuma and Miyata2002). The absolute number of changes among haplotypes was computed in R (v3.4.4; R Development Core Team 2020) using the haplotypes package (v1.1.2; Aktas Reference Aktas2020) while considering gaps as a fifth character state. A haplotype median-joining network was constructed using Network (v10.2; http://www.fluxus-engineering.com). To facilitate visualization, only the D1 individuals from Argentina and Chile were included in the analysis. Each polymorphic site was weighted according to the mutational rates described by Soares and colleagues (Reference Soares, Ermini, Thomson, Mormina, Rito, Röhl, Salas, Oppenheimer, Macaulay and Richards2009). We did not consider mutations 309.1C, 315.1C, A16182c, A16183c, and 16193.1C, or positions with gaps and missing values.
The phylogenetic analysis of D1 mitogenomes was performed with BEAST (v2.6.2; Bouckaert et al. Reference Bouckaert, Vaughan, Barido-Sottani, Duchêne, Fourment, Gavryushkina and Heled2019). Two partitions were considered, one corresponding to the HVR (1–577; 16023–16569) and the other to the coding region (578–16022). Site models were averaged using bModelTest (Bouckaert and Drummond Reference Bouckaert and Drummond2017). As suggested for intraspecific data (Drummond and Bouckaert Reference Drummond and Bouckaert2015), analyses were conducted under a strict molecular clock model, with the prior distribution for the clock rate set as the default. Analyses were done under the coalescent constant population, coalescent exponential population, and coalescent extended Bayesian skyline models. The MCMC chain length was set to 15 million, 20 million, and 65 million iterations, respectively, and data were stored every 1,500th iteration. A clock rate of 0.024694 substitutions per site per million years (Llamas et al. Reference Llamas, Fehren-Schmitz, Valverde, Soubrier, Mallick, Rohland and Nordenfelt2016) and a ploidy of 0.5 were set for the coalescent extended Bayesian skyline. The subsequent model selection was conducted using nested sampling (Maturana Russel et al. Reference Maturana Russel, Brewer, Klaere and Bouckaert2019). For the Bayesian skyline plot, the y-axis transformation to effective population size was performed using a generational time of 25 years (as used in Brandini et al. Reference Brandini, Bergamaschi, Cerna, Gandini, Bastaroli, Bertolini and Cereda2017).
The run parameters were checked with Tracer (v1.7.1; Rambaut et al. Reference Rambaut, Drummond, Xie, Baele and Suchard2018), and trees maximizing the posterior probability were scored with TreeAnnotator (v2.6.2; Bouckaert et al. Reference Bouckaert, Vaughan, Barido-Sottani, Duchêne, Fourment, Gavryushkina and Heled2019). A qualitative analysis of the consensus trees obtained was performed with DensiTree (v2.2.7; Bouckaert et al. Reference Bouckaert, Vaughan, Barido-Sottani, Duchêne, Fourment, Gavryushkina and Heled2019). Graphical postprocessing was done with FigTree (v1.4.4; Rambaut Reference Rambaut2018), and Inkscape (v1.0; www.inkscape.org) was used for image design.
Results
No sign of exogenous DNA contamination was found in the extraction blank, and we ruled out any contamination of the sample with DNA from the researchers who had contact with it. The results of concentration, quantification, and fragment analyses indicate that the library was successfully generated. There are typical postmortem DNA damage patterns that help distinguish ancient sequences from modern contaminants, such as short sequences and cytosine deamination, observed as transitions from C to T at the 5′ end and from G to A at the 3′ end of the readings (Briggs et al. Reference Briggs, Stenzel, Johnson, Green, Kelso, Prüfer and Meyer2007; Jónsson et al. Reference Jónsson, Ginolhac, Schubert, Johnson and Orlando2013; Sawyer et al. Reference Sawyer, Krause, Guschanski, Savolainen and Pääbo2012). Supplemental Figure 1 shows the results obtained with MapDamage2. The excess of purines at the position immediately before the fragment start is a characteristic pattern in ancient DNA (Supplemental Figure 1A); Supplemental Figure 1B illustrates what is expected for a library treated with the USER enzyme, which removes most of these signs of damage.
About 43% of the reads from the enriched sequence PZ3 library passed bioinformatic filters and were mapped against the rCRS (Supplemental Table 3). The average depth achieved was 8X. The incorporation of Sanger sequences allowed us to increase the read depth for two variant positions. The generated consensus mitogenome had 99.7% coverage. A total of 37 SNPs were identified in the haplotype of individual PZ3 with respect to the rCRS (Table 1). Both this number of mutations and their distribution pattern (coding region: 27 mutations / 15,449 total bases; HVR: 10 mutations / 1,120 total bases) were consistent with expected values (Soares et al. Reference Soares, Ermini, Thomson, Mormina, Rito, Röhl, Salas, Oppenheimer, Macaulay and Richards2009), giving robustness to the results.
Variants Found in the PZ3 Consensus Mitogenome.
This haplotype was assigned to the haplogroup D1g and is very similar to two previously reported individuals from the Chilean Antarctic Province: YA2D (de Saint Pierre et al. Reference De Saint Pierre, Bravi, Motti, Fuku, Tanaka, Llop, Bonatto and Moraga2012) and IPY08 (De la Fuente et al. Reference De la Fuente, Ávila-Arcos, Galimany, Carpenter, Homburger, Blanco and Contreras2018). These three haplotypes are characterized by the lack of the A8116 G mutation, originally proposed as a defining characteristic of the D1g haplogroup along with the C16187T reversion (Bodner et al. Reference Bodner, Perego, Huber, Fendt, Röck, Zimmermann and Olivieri2012) but that was later removed in the clade redefinition of de Saint Pierre and colleagues (Reference De Saint Pierre, Bravi, Motti, Fuku, Tanaka, Llop, Bonatto and Moraga2012). The only differences between these and PZ3 are located at positions 309–315 and 3107 (Figure 2 and Supplemental Table 4), which are either hotspot regions or challenging positions for aDNA reconstruction and thus have low or no evolutionary weight.
Bayesian phylogeny depicting the relationships between individuals of the human mitochondrial haplogroup D1 and external groups. The x-axis denotes time into the past in years. The tree is based on an extended Bayesian skyline analysis of the mitogenomes considering two partitions (coding and hypervariable region) performed in BEAST (v.2.6.2; Bouckaert et al. Reference Bouckaert, Vaughan, Barido-Sottani, Duchêne, Fourment, Gavryushkina and Heled2019). Sequence sources and sample numbers are available in Supplemental Table 2. Nodes with posterior probability >0.7 have bars indicating age ranges accumulating 95% posterior probability. The table shows the age ranges and posterior probability for each node. Names of major well-supported sub-haplogroups are indicated. Circles: modern individuals; stars: ancient individuals.
The geographic information of the individuals included in the database (Supplemental Table 2) made possible the Bayesian phylogenetic reconstruction of the mitochondrial D1 lineage for Argentina, Brazil, and Chile. Nested sampling results pointed to the extended Bayesian skyline as the best of the three models considered (ML: −25,739.289; SD: 16.571; Figure 2). The topology obtained strongly suggests that PZ3 is closely related to YA2D and IPY08 (sample numbers 76 and 18, respectively; Supplemental Table 2; posterior probability of 1.00). The estimated age of their hypothetical ancestor is between 203 and 4,439 years. These three individuals were grouped within the D1g haplogroup, whose node age was estimated to be between 12,555 and 17,713 years but with low posterior probability (0.653), a characteristic that was generally observed among sublineage relationships (Figure 2; Supplemental Figure 2).
The haplotype network also supports the inclusion of individual PZ3 within haplogroup D1g, with a derived haplotype from that of YA2D and IPY08; the three were distanced by at least six mutations from the haplogroup's common ancestor (Figure 3). Moreover, this analysis reconstructs D1g sub-haplogroups with good support. The Bayesian skyline plot (BSP; Figure 4) shows an increase in population size between 10,000 and 15,000 years ago, overlapping with the age range obtained by Brandini and others (Reference Brandini, Bergamaschi, Cerna, Gandini, Bastaroli, Bertolini and Cereda2017) for South American mitogenomes.
Haplotype median-joining network of the D1 lineage including individuals from Argentina and Chile constructed using Network (v10.2). Positions 309.1, 315.1, 16182, 16183, and 16193.1, as well as positions with gaps and missing values, were not considered.
Bayesian skyline plot representing historical demographic trends in sampled D1 individuals. Time is reported on the x-axis as years BP, and the effective population size change is reported on the y-axis. Estimations were based on a mutation rate of 0.024694 substitutions per site per million years (Llamas et al. Reference Llamas, Fehren-Schmitz, Valverde, Soubrier, Mallick, Rohland and Nordenfelt2016), and y-axis transformation was performed considering a generational time of 25 years (as used in Brandini et al. Reference Brandini, Bergamaschi, Cerna, Gandini, Bastaroli, Bertolini and Cereda2017).
Discussion
This article presents the first mitogenome entirely generated in Argentina, made possible by the pooled resources and joint efforts of researchers from several local institutions. The analysis performed did not find signs of contamination, and the mitogenome obtained presented the damage patterns expected for ancient DNA samples, although a complete profile could not be achieved due to the use of the USER enzyme during the library assembly. In addition, the pattern of mutation distribution (Table 1) was consistent with expectations for the mitogenome (Soares et al. Reference Soares, Ermini, Thomson, Mormina, Rito, Röhl, Salas, Oppenheimer, Macaulay and Richards2009), attesting to the reliability of our results. The mitochondrial haplotype of individual PZ3 was assigned to the D1g haplogroup, one of the most frequently found in Patagonia (Bodner et al. Reference Bodner, Perego, Huber, Fendt, Röck, Zimmermann and Olivieri2012). In addition, the topology obtained by Bayesian phylogenetic reconstruction (Figure 2) supports this conclusion (posterior probability of 0.653). It also suggests that the PZ3 haplotype was closely related to a Yamana one—YA2D—from Navarino Island, as well as to an individual of unknown age, IPY08, from Hoste Island, located near the Seno Año Nuevo; both islands are in the Chilean Antarctic province (Figure 1).
Although with low posterior probability, the age of D1g node was estimated at between 12,555 and 17,713 years (Figure 2), concordant with the 20.9–11.7 thousands of years range estimated by Roca-Rada and colleagues (Reference Roca-Rada, Politis, Messineo, Scheifler, Scabuzzo, González and Harkins2021). It also partially overlaps with the 18,300 ± 2400 age range estimated for this clade by Bodner and colleagues (Reference Bodner, Perego, Huber, Fendt, Röck, Zimmermann and Olivieri2012), who suggested that the D1g and D1j lineages could have been part of the genetic diversity of the first humans to settle in South America. Our estimate for the age of D1g is also consistent with the genetic and archaeological evidence of the earliest human occupation of this subcontinent (Brandini et al. Reference Brandini, Bergamaschi, Cerna, Gandini, Bastaroli, Bertolini and Cereda2017; Fagundes et al. Reference Fagundes, Kanitz, Eckert, Valls, Bogo, Salzano and Smith2008; Tamm et al. Reference Tamm, Kivisild, Reidla, Metspalu, Smith, Mulligan and Bravi2007). Moreover, the population increase observed between 10,000 and 15,000 years ago in the Bayesian skyline plot (Figure 4) is consistent with the age range for population growth estimated by Brandini and colleagues (Reference Brandini, Bergamaschi, Cerna, Gandini, Bastaroli, Bertolini and Cereda2017) based on multiple South American mitogenomes. According to the authors, this expansion corresponds to the entry of human populations into South America.
From a more general perspective, it is worth noting that most of the D1g sub-haplogroups were reconstructed with posterior probability greater than 0.85 (Figure 2). However, relationships among these sub-haplogroups and between them and the PZ3-YA2D-IPY08 clade cannot be determined with certainty, given the low support of higher-level nodes (the deepest in time). This pattern can be observed throughout the entire phylogeny and can also be appreciated by qualitative analysis of the consensus trees obtained in this study (Supplemental Figure 2), where uncertainties regarding tree topology increase with the node level. This reveals a weak phylogenetic signal at the level of sub-haplogroup relationships, possibly due to the group's youth in evolutionary terms. Alternatively, the weak signal could be an artifact caused by poor sampling of the branches (only a few sequences are available for each sub-haplogroup). To our knowledge, there are no other published data regarding node support for this haplogroup. To provide more accurate answers regarding the population events that could be responsible for the observed pattern, it is necessary to pursue this line of inquiry and expand the database.
The age of the most recent common ancestor shared by PZ3, YA2D, and IPY08 was estimated at between 203 and 4,439 years (Figure 2). This wide age range is probably due to the paucity of information provided by the few and evolutionarily scarcely relevant differences between these three sequences. At that time, populations with maritime economies had already settled and stabilized in the southern Tierra del Fuego Archipelago, possibly as early as 6400 years BP (Orquera et al. Reference Orquera, Legoupil and Piana2011; Orquera and Piana Reference Orquera and Piana2009). The median haplotype network (Figure 3) showed that YA2D, IPY08, and PZ3 mitogenomes, the only insular individuals in our database, are separated from the hypothetical common ancestor shared with the rest of the D1g clade by at least six mutations: C10202T, T10724C, T13020C, T16086C, C16286T, and T16189C (the last one is considered to be a hotspot). The most noteworthy of the mutations are those in positions 16086 and 16286: they have also been reported in studies of the hypervariable region in both ancient and current individuals from the Beagle Channel area (Crespo et al. Reference Crespo, Cardozo, Tessone, Vázquez, Kisielinski, Arencibia, Tackney, Francisco Zangrando and Dejean2020; de la Fuente et al. Reference De la Fuente, Galimany, Kemp, Judd, Reyes and Moraga2015, Reference De la Fuente, Ávila-Arcos, Galimany, Carpenter, Homburger, Blanco and Contreras2018; de Saint Pierre et al. Reference De Saint Pierre, Bravi, Motti, Fuku, Tanaka, Llop, Bonatto and Moraga2012; Motti et al. Reference Motti, Winingear, Valenzuela, Nieves-Colón, Harkins, García Laborde, Bravi, Guichón and Stone2020). In addition, the absence of the mutation in position 8116 is characteristic of Yaghan populations (de Saint Pierre et al. Reference De Saint Pierre, Bravi, Motti, Fuku, Tanaka, Llop, Bonatto and Moraga2012).
These results may contribute to the discussion of at least two issues regarding population dynamics in the Tierra del Fuego archipelago. First, there are three types of archaeological evidence suggesting the existence of interactions between populations in the Beagle Channel and the Brunswick Peninsula-Seno Otway during the mid-Holocene (Figure 1): (1) the presence of green obsidian artifacts in Beagle Channel sites, a material that was frequently used in the Otway area archaeological sites (Álvarez Reference Álvarez2004; Manzi Reference Manzi, Civalero, Fernández and Guráieb2004; Morello et al. Reference Morello, Stern and Román2015; Pallo Reference Pallo2016) and that would probably be obtained in raw form from somewhere in the Brunswick Peninsula (Emperaire and Laming-Emperaire Reference Emperaire and Laming-Emperaire1961; Stern and Prieto Reference Stern and Prieto1991); (2) decoration designs and techniques found on similar artifact types in both areas (Fiore Reference Fiore2006); and (3) the distribution of subfoliar and lanceolate projectile points (Piana and Orquera Reference Piana, Orquera, Morello, Martinic, Prieto and Bahamonde2007). However, the material evidence of this superregional interaction is no longer found in the archaeological record after about 5000–4000 years BP, suggesting spatial circumscription and a reduction in mobility due to several factors (Zangrando Reference Zangrando2009). This change in the archaeological record coincides with our estimate for the maximum age of the most recent common ancestor of PZ3, YA2D, and IPY08. From this moment on, genetic drift may have fixated different lineages in the Beagle Channel and the islands of the western archipelago (de la Fuente et al. Reference De la Fuente, Ávila-Arcos, Galimany, Carpenter, Homburger, Blanco and Contreras2018). This process could also explain the preponderance of the D4h3a lineage over D1g and C1 lineages during the mid-Holocene in the western archipelago, an area associated with Alacalufe groups, and the dominance of C1 in the Beagle Channel population (de la Fuente et al. Reference De la Fuente, Galimany, Kemp, Judd, Reyes and Moraga2015; Moraga et al. Reference Moraga, de Saint Pierre, Torres and Ríos2010).
Second, the genetic proximity between PZ3 and IPY08 suggests the existence of gene flow or a recent common origin among the groups that occupied both sectors of the archipelago during the mid-Holocene. Although the sites where these individuals were recovered are only 90 km apart, there is very little information about the level of interaction between the Beagle Channel and southwestern areas of the archipelago during that period, mainly because there have been few archaeological explorations of Hoste Island. During historic times, Wedell (Reference Wedell2006 [1825]:159–162) indicated the presence of different groups in the Seno Año Nuevo area, yet Bridges (Reference Bridges1886) excluded the southern sector of Hoste Island from the Yaghan territory. Several authors have reported that the geographical distribution of Yaghan populations toward the western end of the archipelago was unclear during the nineteenth century: that area was also occupied by groups that could be ethnically and linguistically linked to Alacalufe populations (Orquera and Piana Reference Orquera and Piana1999). However, as has been suggested for the easternmost sector of the Fuegian archipelago (Crespo et al. Reference Crespo, Cardozo, Tessone, Vázquez, Kisielinski, Arencibia, Tackney, Francisco Zangrando and Dejean2020), the genetic diversity of populations before European contact indicates greater interaction between groups than what is suggested in historical documents. Even if we accept the existence of a genetic continuity between local populations pre- and post-European contact, the cultural geography of the channels and island south of Tierra del Fuego may have been highly dynamic.
Final Remarks
The mitogenome for individual PZ3 was obtained with an average depth of 8X and 99.7% coverage. The quality of the readings was high but there were few of them, which resulted in a relatively low average depth. Adjustments are being made to deal with these shortcomings. The haplotype was assigned to haplogroup D1g and was found to be closely related to YA2D and IPY08, a Yamana from Navarino Island and an ancient individual from Hoste Island. The estimated age for these mitogenomes’ most recent common ancestor was 203–4,439 years (Figure 2), a time in which available archaeological evidence suggests that human groups with maritime economies were already settled in the southern Tierra del Fuego archipelago. The mitogenomes of PZ3, YA2D, and IPY08 present mutations that have only been described in past and extant Yaghan populations from the Beagle Channel area (Figure 3). The maximum age for the common ancestor coincides with a reduction in mobility in Beagle and the western archipelago islands at approximately 4500 years BP.
The construction of the molecular phylogeny made possible the reconstruction of most of D1g sub-haplogroups (Figure 2) and the estimation of age intervals for their nodes that were consistent with results from previous studies. However, the relationships between these sub-haplogroups could not be determined with high certainty due to low support of the high-level nodes.
Finally, although massively parallel sequencing methodologies have been widely applied to aDNA studies, their implementation at the local level still poses many challenges. This first mitogenome entirely generated in Argentina faced some of these difficulties, but thanks to the collaboration of researchers from different local institutions it sets a solid precedent for future studies. This local workflow enabled local researchers to have greater control over sample processing and, in the future, will allow them to define and pursue locally relevant research goals and problems that are not necessarily aligned with research agendas in countries that traditionally led these kinds of studies.
Acknowledgments
We are grateful for the support of the Sección de Antropología Biológica (Instituto de Ciencias Antropológicas) and the Cátedra de Antropología Biológica y Paleoantropología (Departamento de Ciencias Antropológicas) of the Facultad de Filosofía y Letras (Universidad de Buenos Aires, Argentina), as well as for the support of CEBBAD and CCNAA (Universidad Maimónides), INTA, and IMPaM. The Secretaría de Cultura de la Provincia de Tierra del Fuego, Antártida e Islas del Atlántico Sur, granted permits for sample transportation and analysis in 2016. Finally, we thank the journal editor and three anonymous reviewers whose suggestions helped significantly improve this article.
Funding Statement
We appreciate the financial support provided by Fundación Científica Felipe Fiorellino, Fundación de Historia Natural Félix Azara, CONICET, UBA (UBACyT 20020150200233BA and 20020170200363BA), and ANPCyT (PICT 2014-3012).
Data Availability Statement
The consensus mitogenome of individual PZ3 was deposited in GenBank under accession number ON422324. The raw data, in fastq format, are available at the European Nucleotide Archive under accession number PRJEB58193.
Competing Interests
The authors declare none.
Supplemental Material
To view supplemental material for this article, please visit https://doi.org/10.1017/laq.2023.13.
Supplemental Text 1. Materials and Methodology.
Supplemental Figure 1. Plots generated with Mapdamage2 (v2.0; Jónsson et al. Reference Jónsson, Ginolhac, Schubert, Johnson and Orlando2013) to visualize damage patterns. (A) Plots showing the base frequency outside and in the read (the open gray box corresponds to the read); (B) plots with the positions’ specific substitutions from the 5′ (left) and the 3′ end (right). Color codes: red, C to T substitutions; blue, G to A substitutions: gray, all other substitutions; orange, soft-clipped bases; green, deletions relative to the reference; purple, insertions relative to the reference.
Supplemental Figure 2. DensiTree analyses showing the posterior distribution of the trees retrieved by BEAST. Dark regions indicate branch topology agreement among trees.
Supplemental Table 1. Illumina Adapters and Primers Used in this Article.
Supplemental Table 2. Database Compiled and Used in this Article. The information included regarding the haplogroup and the location or population of origin of the sample is the one provided by the paper that published the sequence (or its update in Phylotree17).
Supplemental Table 3. Sequence Data Statistics for the PZ3 Mitogenome.
Supplemental Table 4. Absolute Number of Changes between Haplotypes Performed in R (V3.4.4; R Core Team, 2020) with the Haplotypes Package (v1.1.2; Aktas Reference Aktas2020). Gaps were considered as a fifth character state. Accession numbers (or Sample ID when not available) and haplogroup are listed.
