Genetic Diversity and Demographic History of the Shaggy Soft-Haired Mouse Abrothrix hirta (Cricetidae; Abrotrichini).

Genetic information on species can inform decision making regarding conservation of biodiversity since the response of organisms to changing environments depend, in part, on their genetic makeup. Territories of central-southern Chile and Argentina have undergone a varying degree of impact during the Quaternary, where the response of local fauna and flora was rather species-specific. Here, we focus on the sigmodontine rodent Abrothrix hirta, distributed from 35° S in Chile and Argentina to northern Tierra del Fuego. Based on 119,226 transcriptome-derived SNP loci from 46 individuals of A. hirta, we described the geographic distribution of the genetic diversity of this species using a maximum likelihood tree, principal component and admixture analyses. We also addressed the demographic history of the main intraspecific lineages of A. hirta using GADMA. We found that A. hirta exhibited four allopatric intraspecific lineages. Three main genetic groups were identified by a Principal Component Analysis and by Ancestry analysis. The demographic history of A. hirta was characterized by recent population stability for populations at the northernmost part of the range, while southern populations experienced a recent population expansion.

Introduction

Decision making aimed to conserve biodiversity requires information on several aspects of the biological system to be protected. In particular, genetic information can be key when conservation efforts are made to ensure the survival through time of natural populations and the processes in which they are involved (Supple and Shapiro, 2018). The response of individual organisms to changing environmental conditions can be challenging, and their acclimation and adaptive potential are conditioned by their genetic makeup (Fitzpatrick and Edelsparre, 2018). Given the uncertainty of future environmental conditions, as well as the lack of clear knowledge regarding the link between genes and phenotypes, it might be a safe bet to direct conservation efforts to ensure the preservation, when detected, of distinct intraspecific lineages (Moritz, 2002), particularly in a geographic area where climatic fluctuations modeled current genetic features of the local fauna. The concept of Evolutionary Significant Units (ESU) and Management Units (MU) were early advanced with the aim of conserving the evolutionary heritage and current genetic structure of populations, respectively (Moritz, 1994) (see also Ryder, 1986; Casacci et al., 2014). In fact, one of the main levels of Biodiversity is the genetic level, which emphasizes the need to preserve intraspecific lineages (Gaston and Spicer, 2004).

The biota of southern South America (south of 38°S), roughly corresponding to Patagonia, has been subjected to marked climate changes, glaciations, volcanism, tectonics, and seashore shifts over the past thousands of years (Heusser, 1984; Rabassa et al., 1989; McCulloch et al., 2000). Those historical events have driven diverse processes in the local biota, including population fragmentation and genetic divergence, secondary contact, and range shift. These processes left recognizable genetic footprints, whose characterization is the basis to infer the processes. Since the footprints get erased with time and with the occurrence of subsequent events, such inferences have limitations (Avise, 2009). The body of literature accumulated on the phylogeography of the biota from southern South America, including plants (e.g., Premoli et al., 2010; Achimón et al., 2018; Biersma et al., 2018; Nicola et al., 2019), reptiles (e.g., Victoriano et al., 2008; Femenias et al., 2020), amphibians (e.g., Nuñez et al., 2011), fishes (e.g., Ruzzante et al., 2006; Zemlak et al., 2008; Ruzzante et al., 2020), birds (Norambuena et al., 2018), and mammals (e.g., Lessa et al., 2010; Palma et al., 2012), reveals that the responses to Quaternary dynamics were varied and mostly species-specific (e.g., Lessa et al., 2010; Sersic et al., 2011). Two general patterns were described by Lessa et al. (2010) for Patagonian sigmodontine rodents. On the one hand, species would have persisted in Pleistocene refugia outside Patagonia and would have colonized southern areas from there when climatic conditions became favorable after the Last Glacial Maximum (LGM; ca. 20,000 to 10,000 ya). On the other hand, species would have persisted in one or more Patagonian refugia, whose location would vary according to the species, but may include the Valdivian Refugium at coastal and pre-Andean areas of Chile (from ca. 39°S to ca. 43°S), which has been indicated as an area where animal and plant species persisted during the last glaciation (e.g., see Valdez and D’Elía, 2018).

Abrothrix hirta is a sigmodontine rodent that was revalidated as a distinct species when it was removed from the synonymy of Abrothrix longipilis (Teta and Pardiñas, 2014; Teta et al., 2017). The geographic distribution of A. hirta is wide. Its northern limit in Chile is not clear, but it spans from ca 35° S at both sides of the Andes to northern Tierra del Fuego and from the Pacific coast to pre-Andean areas up to 2,000 masl, including some lowland areas of the Patagonia. Across this large distribution, the species inhabits a wide range of vegetal formations, including sclerophyllous forests in the north, pre-Andean shrublands, temperate rainy Valdivian and Magellan Forests, and the contrastingly arid Patagonian steppes east of the Andes. Most of the knowledge on the natural history and ecology of A. hirta was described referring to it as A. longipilis. However, based on the study areas it is possible, in most cases, to reassign the information to A. hirta (e.g., Patterson et al., 1990; Kelt et al., 1994). More recent studies provided information on A. hirta regarding cranial morphology (Canto et al., 2017), renal transcriptome (Valdez et al., 2015), and digestive morphology (Naya et al., 2014). Cranial variation is not geographically structured (Teta and Pardiñas, 2014), while variation on digestive morphology is associated with habitat type (i.e., forest vs. steppe; Naya et al., 2014). The phylogeographic pattern of A. hirta was addressed in two separate contributions. Palma et al. (2010) examined Chilean populations, where six intraspecific lineages were depicted. The same year, Lessa et al. (2010) studied Argentinean and Chilean specimens of A. hirta, identifying two allopatric clades distributed in northern and southern Patagonia, respectively. Considering both studies, there are no association of clades with the main habitat types, namely forest and steppe; lineages rather replace each other latidudinally. In addition, both Patagonian lineages exhibited signals of demographic expansion, which would have occurred before the LGM (Lessa et al., 2010). In any case, both studies, as well as most of those focused on other Patagonian rodents cited above, were based on a single mitochondrial gene. The single-gene approach is limited in the sense that they do not allow distinguishing between the effect of neutral demographic processes from that of natural selection (Avise, 2009). Considering that natural selection operates on individual or on a few genes at the time, while demographic processes affect the whole genome, widening the sampling to a genomic scale would inform on this issue. This is one of the reasons why it is necessary to transcend genetic studies based on one or a few loci to those at genomic scale (e.g., Lessa et al., 2014). In this line, here we use a wide panel of transcriptomic-derived single nucleotide polymorphisms (SNPs) to study the genetic diversity and demographic history of Abrothrix hirta. We also aimed to assess the intraspecific genetic structure of the species, and to test the hypotheses raised by Lessa et al. (2010), i.e., the existence of two Patagonian genetic lineages that underwent expansion prior to the LGM.

Materials and Methods

Specimen Sampling and Data Collection

The geographic sampling covers most of the distribution of the species, based on 46 individuals of Abrothrix hirta from 14 localities (Figure 1A and Table 1) present at sclerophyllous shrublands, Valdivian and Magellanic Forests, and Patagonian steppe. Data of 16 specimens from four localities were taken from Valdez et al. (2015), while data on the other 30 individuals from other 10 localities were newly gathered in this study (Table 1). Several of these newly sequenced specimens were already housed at scientific collections. Part of the Chilean specimens were newly collected in this study using Sherman-like traps baited with oatmeal. The taxonomic identity of captured individuals was first determined in the field based on the external morphology; this was later corroborated in the laboratory with sequences of the mitochondrial gene Cytochrome-b as in D’Elía et al. (2008). Immediately after euthanasia, the right kidney was removed and conserved in RNAlater stabilization solution (Invitrogen). Specimen vouchers and tissue samples were deposited at the Colección de Mamíferos, Universidad Austral de Chile, Valdivia, Chile. All procedures involving live animals were conducted following guidelines of the American Society of Mammalogists (Sikes and Animal Care and Use Committee of the American Society of Mammalogists, 2016) and were approved by the Uso de Animales en la Investigación committee of the Universidad Austral de Chile (permits 04/11 and 296/2017); the Servicio Agrícola y Ganadero authorized the collection of individuals (permits 1231/2017, 5611/2013, 2164 and 5165).

FIGURE 1

Phylogeographic pattern of Abrothrix hirta. (A) Map of collection localities of Abrothrix hirta included in this study. Collection localities are color-coded as in Figure 3. Circles and ellipses indicate the main lineages shown in panel (C). (B) Groups of localities and color-coded localities within each group: G1 = 2, G2 = 4; G3 = 5-6, G4 = 7, G5 = 8-10, G6 = 11-13, G7 = 14. Localities included in each group (G1–G7). Details on each locality are provided in Table 1. (C) Maximum likelihood tree for samples of A. hirta. Node labels correspond to Bootstrap support values. Terminals are color-coded as in panels (A) and (B). Lineages are labeled as described in the text.

TABLE 1

Specimens of Abrothrix hirta included in this study.

#	Locality	Specimen ID	Group for Fst-based analyses
1	Chile, Región del Maule, Constitución, Quivolgo; -35.315; -72.395.	UACH 8190	Not included
2*	Chile, Región del Maule, Cauquenes, Pelluhue; -35.80165; -72.5301.	UACH 7859 UACH 7863 UACH 7866	G1
3	Chile, Región de la Araucanía, Malleco, Cordillera Las Raíces, Malalcahuello; -38.419535; -71.505124.	UACH 8180	Not included
4*	Chile, Región de Los Ríos, Valdivia, San José de la Mariquina, Fundo San Martín; -39.651; -73.19.	UACH 7310 UACH 7311 UACH 7312 UACH 7313	G2
5*	Chile, Región de los Ríos, Valdivia, Valdivia, Chabelita; -39.776967; -73.304783.	UACH 7314 UACH 7315, UACH 8316, UACH 8317, UACH 8318	G3
6	Chile, Región de los Ríos, Valdivia, Isla Teja, Campus Universidad Austral de Chile; -39.795101; -73.260741.	UACH 8320	G3
7*	Argentina, Provincia de Chubut, Telsen, Gan Gan, Establecimiento La Maroma; -42.694353; -68.233911.	UACH 8329, UACH 8330, UACH 8331, UACH 8332, UACH 8333, UACH 8334	G4
8	Argentina, Santa Cruz, Río Chico, Cajón del Río Oro, 3 km al O por camino del puente Río Oro; -47.421050; -71.958233.	UACH 8335, UACH 8336, UACH 8338, UACH 8339, UACH 8340	G5
9*	Argentina, Santa Cruz, Río Chico, La península, 2 km al SSO por camino de Tío Camping (EA, La Península); -47.447417; -71.903483.	UACH 8337	G5
10	Argentina, Santa Cruz, Río Chico, Valle del Río Oro, 5 km por camino del puente Río Oro; -47.422150; -71.981567.	UACH 8341, UACH 8342, UACH 8343, UACH 8344	G5
11*	Chile, Región de Aysén, Sector Barrancoso 40 km al sur de Cochrane por ruta 7 (carretera austral); -47.496717; -72.808617.	UACH 8306, UACH 8307, UACH 8308, UACH 8309, UACH 8310, UACH 8311, UACH 8312	G6
12	Chile, Región de Aysén, Sector Barrancoso 41,5 km al sur de Cochrane por ruta 7 (carretera austral); -47.500267; -72.821050.	UACH 8313	G6
13	Chile, Región de Aysén, Sector Barrancoso, 35 km al sur de Cochrane por ruta 7 (carretera austral); -47.477717; -72.794483.	UACH 8314, UACH 8315	G6
14*	Chile, Región de Magallanes, Punta Arenas, Puerto del Hambre; -53.606; -70.937.	UACH 8323, UACH 8324, UACH 8325, UACH 8326, UACH 8346	G7

FIGURE 3

Principal component analysis based on 119,226 SNP loci of 46 specimens of Abrothrix hirta. Each specimen was colored according to the locality at which it was collected as indicated in Figure 1 (details in Table 1). Lineages found in the phylogenetic analysis (Figure 1) are indicated with the ellipses. (A) PC1 vs PC2; (B) PC2 vs PC3; (C) PC1 vs PC3.

Transcriptome Sequencing and SNP Calling

Total RNA extraction, messenger-RNA purification, and complementary-DNA construction were conducted using RNeasy mini kit (Quiagen), PolyATract mRNA Isolation System II, and TruSeq RNA Sample Preparation Kit, respectively, following Valdez et al. (2015). Libraries constructed for each specimen were subjected to high-throughput sequencing following Illumina HiSeq2000 protocol (Illumina Inc.) to obtain paired-end reads with a length of 101 bp.

Transcriptome data from each of the 46 samples were analyzed individually. Raw reads were trimmed to remove Illumina adapters and low quality reads using Trimgalore tool[1]. Thresholds for minimum quality and read length were set to 20 phred and 25 bp, respectively. Reads derived from ribosomal RNA (rRNA) were excluded from the dataset; these reads were identified mapping trimmed reads using Bowtie version 2.2.6 (Langmead and Salzberg, 2012) against a reference containing rRNA sequences from Rodentia species available at GenBank[2]. Filtered read datasets were assembled de novo using Trinity version 2.8.6 (Grabherr et al., 2011). Assembled contigs were annotated using the Blastx alignment tool (Camacho et al., 2009); a dataset containing amino acid sequences of Mus musculus from OMA[3] was used as a reference for gene annotation. After annotating individual assemblies, the best contig was selected for each gene; i.e., the longest contig that spans at least 80% of the amino acid sequence of the corresponding reference protein. As in Giorello et al. (2014), custom scripts were used in this step. The set of such best contigs were used as a reference transcriptome for Single Nucleotide Polymorphism (SNP) calling.

For each sample, SNPs were identified using GATK v.4.1.4.0 (McKenna et al., 2010), following the analytical workflow recommended for RNA-seq data. First, the STAR’s two-pass mode (Dobin et al., 2013) was applied to individual read files against the reference transcriptome. This accounts for possible splice junctions in the dataset. Following this step, GATK’s SplitNCigarReads was used to split reads into exon segments and hard-clip any sequences overhanging into the intronic regions. Duplicated reads that are likely to have been generated due to artifactual processes (such as PCR duplicates) were removed using MarkDuplicates tool from Picard (Wouters et al., 2002). GATK’s HaplotypeCaller was used to generate individual vcf files for each sample. All 46 vcf files were then merged into one file using BCFtools (Danecek and McCarthy, 2017). Only SNP loci present in all 46 individuals were included in downstream analyses. Loci departing from Hardy-Weinberg equilibrium, multiallelic loci, and having a minor allele lower than 0.01 were also filtered out. A total of 119,226 SNP loci were kept after filtering.

Datasets Used in the Analyses

Three datasets were constructed; all of them containing biallelic SNP loci that were in Hardy-Weinberg equilibrium, present in frequencies higher than 0.01. Filtering of the data was conducted using bcftools[4]. The first one (hereafter dataset 1) contains 119,226 SNP loci from all 46 sampled individuals of Abrothrix hirta (Figure 1A and Table 1). The second dataset (hereafter dataset 2) includes 119,226 SNPs loci from 44 individuals from 12 localities; it excludes samples from localities 1 and 3 (Figure 1A and Table 1) because only a single specimen was included from these localities. Dataset 2 was used for those analyses that require multiple sampling per locality (see below). To improve statistical power, geographically close localities were merged into 7 groups in the dataset 2 (G1–G7; see Figure 1B and Table 1). G1 includes samples from the Chilean region of Maule; G2 and G3 contain samples from Los Rios Region. Groups G4 and G5 are distributed in the Argentinean provinces of Chubut and Santa Cruz, respectively (these are all the samples from open arid areas east of the Andes); G6 includes samples from the Chilean Region of Aysén. Finally, G7 is form by samples from the Magallanes Region of Chile (Figure 1C and Table 1). The third dataset (hereafter dataset 3) contains 220,080 SNPs loci from the 46 individuals of A. hirta (Table 1), and a one individual of the closely related species Abrothrix longipilis (UACH 8097), A. manni (UACH 7876), and A. sanborni (UACH 8322); the last three were used to form the outgroup, since they are the most closely related species to A. hirta (Cañón et al., 2014; Teta et al., 2017).

Genetic Diversity

Nucleotide diversity (Pi) was calculated using dataset 2, where a single value was obtained for each locality group (G1–G7). To calculate Pi, the average number of nucleotide differences between two randomly chosen samples from within groups of populations was estimated using PopGenome (Pfeifer et al., 2014); this number was then divided by the total number of SNP loci in the dataset. The observed heterozygosity (oHe) for each site and the mean oHe for each population group (G1–G7; Table 1) were calculated using the R package pegas (Paradis, 2010). The fixation index Fst was calculated among pairs of groups using PopGenome.

Intraspecific Lineages

In order to detect and describe the main intraspecific lineages of Abrothrix hirta, a maximum likelihood tree was inferred using dataset 3 and IQ-TREE (Nguyen et al., 2015). The best-fit substitution model for dataset 3 was TM + F + I + G4 + ASC, which was selected based on the Bayesian Information Criterion and implemented in the ML tree reconstruction. Branch support was estimated using 1000 replicates and ultrafast bootstrap (Hoang et al., 2018). An additional step was run to further optimize ultrafast bootstrap trees by nearest-neighbor interchange to account for model violations.

Genetic Structure

Two approaches were implemented to explore the genetic structure of A. hirta. The first one, using dataset 1, consisted in a principal component analysis (PCA) conducted with Plink 1.9 (Purcell et al., 2007). The top 20 principal components of the variance-standardized relationship matrix were extracted. The eigenvectors with the first three principal components were plotted. The second approach consisted in a maximum likelihood estimation of individual ancestry, which was conducted using the R package LEA (Frichot and François, 2015). Individual admixture coefficients were estimated from the genotypic matrix assuming K ancestral populations. Ten independent runs were performed with K = 1 to K = 10 and calculating individual ancestry proportion and a cross-entropy criterion. The run yielding the value of cross-entropy criterion has the number of ancestral population (K) that best explains the genetic dataset (Alexander and Lange, 2011).

In order to test for a pattern of isolation by distance, Mantel statistics were calculated based on Pearson’s product-moment correlation using Fst estimates between pairs of groups (G1–G7), calculated using dataset 2, and geographic distances among one representative of each group of localities and estimated in decimal degrees using the package sp (Pebesma and Bivand, 2005).

Demographic Modeling

Here the units of analysis were the four uncovered main lineages (see below). One of these (the so-called Araucarian lineage; see results, Figure 1C) was excluded because it is composed by a single sample (UACH 8180) precluding its consideration in this analysis. The site frequency spectrum (sfs) was estimated using easySFS[5]. A projection of 8, 8, 8 was used to generate a joint sfs file, which was the input file for the demographic analysis using GADMA (Noskova et al., 2020). The method of ordinary differential equations was applied using moments software (Jouganous et al., 2017) implemented in GADMA. Running parameters where Theta0: 0.37976, which was calculated as 4∗mu∗L, as in Gutenkunst et al. (2009); time for generation: 1.0, number of repeats: 30; the remaining parameters were set as default. Several combinations of initial and final structure were previously run and the combination of Initial structure: 1, 1, 1, Final structure: 2, 2, 2; showed the best likelihood value. Demographic models with this structure were evaluated calculating parameters of (a) time of population divergence, (b) effective population size, and (c) migration, for three a priori defined populations, namely, the Maule lineage, the Valdivian Lineage, and the Patagonian lineage. Heuristic searches by means of the genetic algorithm (GA) implemented by GADMA was used to find the best-fit parameter values in order to automatically infer the best demographic model from the joint sfs data. Composite Likelihood Akaike Information Criterion (CLAIC) was used to choose the best model considering the three populations in one.

Results

Genetic Diversity of Abrothrix hirta

The genetic diversity within locality groups (G1–G7; Figure 1B and Table 1) was higher at northern localities and gradually decreased toward the south (Figure 2); this pattern was detected with both Pi and oHet estimates. Pi values ranged between 0.01 (for G7) and 0.12 (for G1; Figure 2A and Supplementary Table 1); oHet for individual loci ranged from 0 to 1, but the frequency of each value varied (Figure 2B) and the mean oHet decreased at higher latitudes (Supplementary Table 2).

FIGURE 2

Molecular diversity for seven groups of populations of Abrothrix hirta. (A) Molecular diversity index (Pi) calculated for seven groups of localities (G1- G7; Figure 1; Table 1). The values for this graph are given in Supplementary Table 1. (B) Observed heterozygosity across individuals within groups. The data for this graph is given in Supplementary Table 2.

Pairwise estimates of Fst between group pairs ranged from 0.04 for the comparison of G2 and G3, which are relatively close to each other in the Valdivian forests of the Los Ríos Region, to 0.28 for the comparison between the distantly located groups G4 in the Patagonian steppes of Gan Gan, Chubut, Argentina and G7 in the Magellanic forests of Puerto del Hambre, Magallanes, Chile (Supplementary Table 3).

Intraspecific Lineages in Abrothrix hirta

The ML tree of 46 individuals of Abrothrix hirta (Figure 1C) strongly supports the monophyly of the species (BS = 100). Four main intraspecific lineages were recovered, which are geographically segregated and latitudinally ordered. Specimens from the northernmost localities in the Region del Maule (localities 1 and 2: G1; Figure 1 and Table 1) form a strongly supported clade (BS = 100), referred hereafter as the Maule lineage. Specimens from Región de Los Ríos (localities 4, 5, and 6: G2 and G3; Figure 1A and Table 1) form a strongly supported clade (BS = 100), referred here as the Valdivian lineage. A single sample from Malalcahuello (locality 3, Figure 1 and Table 1) comprises the Araucarian lineage. Finally, unlike the previously mentioned lineages, the last main lineage, referred here as the Patagonian lineage, has a broad geographic distribution; it includes individuals from locality 7 (G4) in Gan Gan, Argentina, southwards to localities 8-10 (G5) in Santa Cruz, Argentina, 11-13 (G6), Aysen, Chile and 14 (G7), Puerto del Hambre, Magallanes Region, Chile. Interestingly, specimens of each of these four locality groups form monophyletic groups. Relationships among main lineages are as follows (Maule (Valdivian (Araucarian, Patagonian))); all relationships are highly supported (BS = 100). Branch lengths in the ML tree showed that the Maule and Valdivian lineages contain samples that are more divergent among them in comparison to those in the Patagonian lineage (Figure 1A).

Genetic Structure of Abrothrix hirta

The PCA conducted with 119,226 SNP loci from all 46 individuals indicates that the species genetic variation is geographically structured (Figure 3 and Supplementary Tables 4, 5). The first three principal components explained 32.41% of the observed variance (Supplementary Table 4). In the two-dimensional space composed by PC1 and PC2 (Figure 3A), four groups of specimens are well segregated in the space; these groups are concordant with the main clades obtained in the ML tree (Figure 1C). The first PCA group, which is the one most segregated from the other groups along PC1, contains samples from localities 1 and 2, corresponding to Maule lineage. This group has a higher dispersion along PC3 than along PC1 and PC2 (Figure 3C). At low PC1 and high PC2 values (Figure 3A) is found a second group of samples from localities 3-6, corresponding to the Valdivian lineage. Samples from localities 7-14, corresponding to the Patagonian lineage, are closely grouped at low PC1 and PC2 values. The single sample from the Araucarian lineage (locality 3) approaches those of the Patagonian lineage (Figure 3A), but in PC2 vs. PC3, it appears closer to samples of the Maule lineage. The segregation of Valdivian and Patagonian samples is low along PC3 (Figure 3C). It is noteworthy that samples of the Patagonian lineage are closely grouped along the first three principal components.

A pattern of IBD was detected by the Mantel test; there is a positive and significant correlation (r = 0.6232; p-value = 0.007996) of Fst between pairs of groups (G1-G7) and their geographic distances (Supplementary Tables 3, 6, respectively).

The analysis of population structure with LEA found that the best number of ancestral population (i.e., genetic clusters) in the dataset is K = 3, which yielded the lowest value of cross-entropy (Supplementary Figure 1). A scheme of K = 4 was the second best. Note that the single individual representing the Araucarian lineage was not included in this frequency-based analysis. Individuals from G1 share high coefficients of ancestry corresponding to cluster 1 (Figure 4); moreover, cluster 1 is almost restricted to the Maule lineage (G1), except for the low proportions of this cluster present in individuals from G2, G3, G4, and G7. Individuals from G2 and G3 present high proportion of ancestry corresponding to cluster 2 and small proportions of clusters 1 and 3 (i.e., small sections of red and green in otherwise blue sections; Figure 4). Despite this, cluster 2 roughly approximates the Valdivian lineage (Figure 1C). Similarly, cluster 3, which is present in high proportions of ancestry in individuals from G2 to G7, roughly corresponds to the Patagonian lineage (Figure 1C). In the scheme of K = 4, the single difference consists in a split of cluster 1 into two. The index of fixation, Fst, among the three clusters were very similar; Fst of cluster 1 (red in Figure 4) vs. cluster 2 (blue in Figure 4) = 0.271; Fst of cluster 1 vs cluster 3 (green in Figure 4) = 0.277; and Fst of cluster 2 vs. Cluster 3 = 0.270.

FIGURE 4

Proportions of ancestry for specimens of Abrothrix hirta from seven groups of localities into three genetic clusters. The three genetic clusters (K = 3; see Supplementary Figure 1) are differently colored (red, blue, and green). The black bar below indicates the limits of each of the seven population groups (Figure 1). Samples and localities included in each group are specified in Table 1.

Historical Demography of Abrothrix hirta

The best demographic scenario estimated with GADMA for the three main lineages of A. hirta for which more than one specimen was included in the analysis, involves periods of population size stability and change, as well as events of gene flow between lineage pairs; it is described below in six-time intervals (separated by tick marks in the x-axis in Figure 5). The ancestral population of A. hirta was stable, with an initial Ne of 43,604 individuals until 207 Kya when it underwent a linear growth up in size until 198 Kya when this ancestral population split into two lineages. On the one hand arose the Maule lineage (Figure 5) and on the other hand the ancestral population of the current Valdivian and Patagonian lineages. The Maule lineage experienced a steady reduction of size until 42 Kya, to persist at minimal values until it underwent a sudden population growth 15 Kya ago, since when it maintained a stable size until the present. The ancestral population of Valdivian and Patagonian lineages also experienced a marked and continuous reduction in size until 42 Kya when it split into current Valdivian and Patagonian lineages. Since its origin, the Ne of the Valdivian lineage remained stable, while the Patagonian lineage experienced exponential population growth during the last ca. 7,500 years. Additionally, during the last ca. 40 Kya, events of bidirectional and asymmetric gene flow between Maule and Patagonian lineages, as well as between Valdivian and Patagonian lineages have taken place (black arrows in Figure 5).

FIGURE 5

Demographic history for three main lineages of Abrothrix hirta. (A) Demographic model of the three lineages of A. hirta: ML, Maule lineage, VL, Valdivian lineage, and PL, Patagonian lineage. The number at the top-right scales the y-axis that corresponds to effective population size. In the x-axis, time is indicated in thousands of years before present. Main glaciation events are indicated in the time axis; MQG, Middle Quaternary Glaciation, ELG, Early Late Quaternary Glaciation; MLG, Middle Late Quaternary Glaciation, LGM, Late Glacial Maximum. (B) Allele frequency spectrum distribution of sfs files (first line); theoretical model sfs (second line); residuals (third line); model adjustment for each lineage (fourth line).

Discussion

Abrothrix hirta is one of the most abundant rodents in southern South America (Patterson et al., 2015). However, the northern limit of its distribution remains unclear. As a result of our field work, here we extend the known distribution of the species toward the norther up Constitución, Maule Región include (locality 1; Table 1 and Figure 1) where it almost reaches the southern known locality of A. longipilis (see Valdez et al., 2020). A. hirta inhabits highly contrasting environmental conditions, including sclerophyllous forests, the humid temperate Valdivian and Magellan rainforest, and the arid Patagonian steppe. This wide geographic distribution includes areas that during the Quaternary were impacted, at varying degrees of intensity, by climatic fluctuations. Using a wide SNP panel derived from transcriptomes, we identified intraspecific lineages and described the genetic structure and demographic history of A. hirta. With this approach, we are contributing to undertake the transition promoted by Lessa et al. (2014) of studies on South American rodents toward the genomic era.

Genetic Variation in Abrothrix hirta Across Its Geographic Distribution

The fact that Abrothrix hirta presents its genetic variation geographically structured was consistently depicted by three analytical approaches, namely, ML reconstruction, PCA, and Ancestry analysis. These lineages are the Maule, Araucarian, Valdivian, and Patagonian lineages (Figures 1C, 3, 4). The Araucarian lineage is represented by A single specimen (UACH 8180) from Malalcahuello (locality 3), whose distinction was shown in both the ML tree and the PCA (Figures 1C, 3). This specimen was not included in the Ancestry analysis because it requires multiple samples per each locality. Further sampling in pre-Andean areas of the Araucanía Region will help test the genetic distinction of this geographic area which is shown by specimen UACH 8180.

Other small mammals that are co-distributed with Abrothrix hirta also exhibit geographic structure. Albeit in a much smaller distribution area, the congeneric Valdivian forest endemic A. manni presents two main clades (Valdez and D’Elía, 2018). Another congeneric species, the olive soft-haired mouse A. olivacea presents three main clades in the study area, one in central and northern Chile, and the others in Patagonia and in Tierra del Fuego (Lessa et al., 2010). The southern pericote Loxodontomys micropus also presents segregated lineages in Patagonia (Cañón et al., 2010). Unlike these, several other small mammals within the geographic area of Abrothrix hirta, such as Eligmodontia morgani, Reithrodon auritus, and Phyllotis xanthopygus (Lessa et al., 2010), do not exhibit phylogeographic breaks. However, all these studies were based on unilocus mitochondrial genetic variation, with the consequent limitation of uncertainty about the possible distortion of demographic inferences due to the effect of natural selection. Our SNP-based approach allows us to confidently state that A. hirta presents distinct genetic lineages across its distribution.

Several phylogeographical breaks were detected in this study; however, given the gaps in our geographic sampling, we focused attention on two breaks located approximately in the same areas as other species also show phylogeographic breaks. The northernmost break separates the Maule lineage from the Araucarian, Valdivian and Patagonian lineages on the other (Figure 1). A finer-grain sampling is needed to precisely locate this break, but it roughly coincides with the general area of the break described for the small cat Leopardus guigna (Napolitano et al., 2014), the frog Eupsophus roseus (Correa et al., 2017), and the lizard Liolaemus pictus (Vera-Escalona et al., 2012). Since divergence times are not provided for all of the above studies, it is not possible to establish whether the divergences occurred at the same time, and therefore whether they were produced by the same event. At this latitude lie the rivers Itata and Bio Bio. For a small mammal such as A. hirta, the riverine barrier might have acted as a primary barrier, bisecting a previously continuous population, or as a secondary barrier, preventing the gene flow between already differenced populations (Patton and da Silva, 1998). In the case of Abrothrix hirta, considering that the rivers are most likely older than the divergence in the basal dichotomy (ca. 198 Mya; Figure 5) it is possible that a founder effect, likely from Central Chile southwards across these rivers would have promoted the differentiation between the Maule lineage and the Araucarian, Valdivian, and Patagonian lineages.

Another phylogeographic break of Abrothrix hirta is located at ca. 43°S, at Los Ríos Region in Chile and Río Negro and northern Chubut in Argentina (Figure 1). A break at this latitude was also observed in the olive soft-haired mouse Abrothrix olivacea (Lessa et al., 2010) and in the huemul Hippocamelus bisulcus (Marín et al., 2013a). The latitude 43° S is also an area of contact of two widely distributed lineages of the guanaco Lama guanicoe (Marín et al., 2013b). The frogs Batrachyla leptopus (Nuñez et al., 2020), Pleurodema thaul (Barria et al., 2020), and the lizard Liolaemus petrophilus (Fontanella et al., 2012) also present phylogeographic breaks in this area. Several of the species mentioned above, diverged into intraspecific lineages and subsequently showed demographic expansion of at least one of them toward Patagonia after the LGM. It is possible that the split depicted at this phylogeographical break would have originated by genetic divergence of populations isolated in at least two refugial areas during the Pleistocene glaciations in Central-southern Chile and Argentina, from where they expanded southwards. This hypothesis can be tested using denser spatial sampling and demographic modeling approaches.

Within the study area, several processes have been proposed as possible causes of the generation and/or maintenance of intraspecific lineages in terrestrial vertebrates, such as a riverine effect in species of mammals, reptiles, frogs, and plants (e.g., Torres-Pérez et al., 2007; Sallaberry-Pincheira et al., 2011; Vásquez et al., 2013; Palma et al., 2014; Viruel et al., 2014; Muñoz-Mendoza et al., 2017), habitat fragmentation prompted by climatic changes (e.g., Lessa et al., 2010; Valdez and D’Elía, 2018) and isolation-by-distance followed by population differentiation (Sersic et al., 2011). Based on the phylogeographic pattern described here, a combination of the mentioned processes is plausible to explain the genetic footprints of Abrothrix hirta. In particular, the fact that the genetic variation of A. hirta is nested within that of the Maule populations suggests a history of colonization from this area toward the south.

The main dichotomy within the Patagonian lineage is consistent with the scheme of northern-southern Patagonian clades noted in Lessa et al. (2010). Samples from Gan Gan (locality 7) comprise a northern lineage that extends northwards in Argentina up to Lago Quillén, Chubut, an area that is not covered by our sampling. Lessa’s souththern clade coincided with our Patagonian clade (excluding samples from locality 7). This general pattern of latitudinally replacing clades contrasts with the variation in A. hirta in terms of morphological variation that are associated with water economy. Naya et al. (2014) showed that populations from mesic Valdivian and Magellanic forests at the west slope of the Andes have larger body size and longer intestines relative to those from the xeric Patagonian steppe, east of the Andes. The same contrast between genetic structure and eco-morphological traits related to water economy was documented in the olive soft-haired mouse Abrothrix olivacea (Lessa et al., 2010; Giorello et al., 2018).

The range of Abrothrix hirta faced marked climatic changes during the Quaternary (Heusser, 1984; Rabassa et al., 1989). Glaciation cycles strongly impacted the species range, particularly in Patagonian. For instance, during the LGM the ice cap covered part of the Chilean territory, from sea level in southern Chiloé, northwards through the regions of Los Ríos and Los Lagos up to ca. 33°S (McCulloch et al., 2000). Coastal areas of Los Lagos, Los Ríos, and lowlands up to Bio Bio were ice-free, although areas that are currently forested had undergone a transition to subantarctic parklands and woodlands (Heusser et al., 1996; Moreno et al., 2001). In our analyses, the genetic footprints of Abrothrix hirta can be traced back in time up to ca. 210 Kya (Figure 5A), i.e., before the Middle Quaternary Glaciation (190–130 Kya). At that time, the ancestral population experienced an exponential growth before the divergence that originated the Maule lineage and the ancestor of the other lineages. During the Middle Quaternary Glaciation (190–130 Kya), the Early Late Quaternary Glaciation (98–75Kya), and the Middle Late Quaternary Glaciation (43–33 Kya) and the interglacial intervals among them, a progressive demographic reduction marked the history of the Maule lineage as well as that of the ancestor of the other lineages. Before and during the LGM population sizes persisted in markedly small populations in the Maule, Valdivian, and Patagonian lineages. Noteworthy is the fact that there has been gene flow during the last ca. 40 Ky, with strong levels roughly before and after the LGM, although genetic interchange between the Valdivian and Patagonian lineages also occurred during the last glaciation. Taking into consideration the geographic distribution of the Maule and Valdivian lineages in areas that were not impacted by late Pleistocene climate as in southern areas, it is expectable to register gene flow among those lineages. It is noteworthy that higher intensity of gene flow was inferred from the Valdivian lineage toward the Maule (see arrows in Figure 5). Regarding the Patagonian lineage, it is likely that the gene flow detected during the LGM would have taken place at the northern edge of its range. Again, finer-scale sampling would clarify these events. However, although geographically limited, the dataset used in this study provided information that single-locus studies cannot; for this reason, the scenario presented here represents a new generation of hypotheses about the demographic history of Patagonian rodents.

Since the LGM, populations in the northern part of the range (ca. 34° S–42°S) remained stable, while in the south (ca. 43°S–54°S), populations have expanded. The Maule and Valdivian lineages maintained their effective population size over the last ca. 15 Ky (Figure 5). In a neutral framework, a population that has remained stable through time accumulates more genetic variants, which increases its divergence due to random mutations. Concordant with this expectation are the results showing that the genetic diversity was higher at northern localities (Figure 2), and that branch lengths in the ML tree are longer in Maule and Valdivian lineages than those of the Patagonian lineage (Figure 1C), despite the geographic distribution of the latter is markedly larger than the distributions of Maule and Valdivian lineages (Figure 1A). Additionally, the genetic dispersions of the specimens that compose the Maule and Valdivian lineages are larger than that those shown by the specimens of the Patagonian lineage (Figure 3). Conversely, our results suggest a process of post-glaciation colonization of open and Magellanic forest areas of Patagonia, as population in these areas experienced a recent population increase (Figure 5), showed a pattern of isolation-by-distance), and are genetically less diverse (Figure 2; see also branch lengths in Figure 1C).

Putative Pleistocene refugial areas, from which populations expanded toward their current distribution range, have been identified in central-southern Chile. This is the case for small mammals, including the mice Oligoryzomys longicaudatus (Palma et al., 2012) and Abrothrix manni (Valdez and D’Elía, 2018), the amphibians Batrachyla leptopus (Nuñez et al., 2020) and Rhinella arunco (Vásquez et al., 2013); the reptiles Liolaemys pictus (Victoriano et al., 2008; Vera-Escalona et al., 2012) and L. tenius, L. lemniscatus (Victoriano et al., 2008) and the plants Hypochaeris palustris (Muellner et al., 2005) and Nothofagus pumilio (Mattera et al., 2020). Most of the studies reported demographic expansions after the LGM, indicating progressive population growth in the area over the past ca. 15,000 years. Moreover, besides the population increase, the pattern of colonization of Patagonian areas from northern refugial areas was inferred in other sigmodontine rodents, including Abrothrix olivacea (sensu Lessa et al., 2010; Cañón et al., 2014), Eligmodontia morgani, Phyllotys xanthopygus, and (Smith et al., 2001; Lessa et al., 2010). Although the latter also presents signs of population persistence in another refugium at Tierra del Fuego. Also, plant species exhibited similar historical patterns (e.g., Hypochaeris incana Tremetsberger et al., 2009; Echinopsis chiloensis Ossa et al., 2019). Taken as whole, three phylogeographic studies focused Patagonia revealed that the response of the biota to climatic oscillations is mostly species-specific. This fact justifies the efforts to unravel species genetic patterns not only to increase knowledge on local dynamics but also to make informed conservation decisions.

Taxonomic Implications

The genetic patterns depicted in this study identified distinct lineages (Figure 1) among which there is gene flow (Figure 5); this is consistent with a taxonomical scheme where a single species is recognized. However, a scenario of multiples hybridizing species cannot be ruled out. Considering a single species framework, intraspecific lineages identified here might be recognized at the subspecies level. A study of the morphological and morphometric variation within and among the lineages described here will help determine the most appropriate taxonomical arrangement for these populations. Several taxa (e.g., apta, moerens, suffusus) are included under the synonymy of A. hirta. Our geographic coverage is limited and as such, it does not allow a proper association between available names and our identified lineages. Further studies on morphological and genetic variation are required to address this issue.

Implications for Conservation

Phylogeographic and population genomic studies provide the way to uncover the geographic distribution of genetic diversity and the demographic history of target species. As early stated by Moritz (1994), when advancing the notion of Evolutionarily Significant Units, understanding the evolutionary history of the species, the number, composition, and distribution of intraspecific lineages are relevant for conservation purposes (see also Ryder, 1986; Supple and Shapiro, 2018). Each lineage, which besides taxonomic status can be considered as distinct ESU, has its own combination of genetic variants (including those that discriminate the genetic groups in Figures 1, 3, 4), therefore a distinct genetic pool, which in turn would allow different responses to future environmental changes. Additionally, it has been shown that the distinctive properties of intraspecific lineages affect the projections of species distribution models (i.e., the sum of the areas resulting from models built for each lineage separately is distinct to the area resulting when modeling the species as a whole) and therefore the accuracy of predictions in the context of the environmental change due to global warming and anthropogenic impact (Barria et al., 2020).

While the conservation status of Abrothrix hirta in Argentina is considered as of low concern (Teta and D’Elía, 2019), it has not been categorized by neither the International Union for Conservation of Nature (IUCN) nor the Chilean authorities via the Ministerio del Medio Ambiente. Populations of this species are still included under Abrothrix longipilis, following an old taxonomic scheme abandoned seven years ago (Teta and Pardiñas, 2014). This is problematic not only because it greatly overestimates the distribution of A. longipilis, but also because the true genetic diversity and distribution of the intraspecific lineages of A. hirta remain hidden from conservation concerns. Our results indicate that two of the four main lineages of A. hirta occur Chilean in protected areas. This is the case of the Araucarian lineage, which is sheltered by the national protected area Reserva Nacional Malalcahuello (locality 6; Figure 1A). It is also possible that individuals from this lineage are present in other nearby Chilean and Argentinean protected areas. The Valdivian lineage is also being protected; we collected individuals at Fundo San Martín, a small private conservation property of the Universidad Austral de Chile (locality 4; Figure 1A). Again, it is also possible that other protected areas harbor populations of this lineage in Chile and/or Argentina. We did not collect the Patagonian lineage at any protected area but given its large distribution and the several large protected areas existing in both Argentinean and Chilean Patagonia, we suggest that the status of this lineage to be of least concern. However, it is of concern to determine whether the Maule lineage is safeguarded. This is the most divergent lineage of Abrothrix hirta and the one with the highest genetic diversity (Figures 1, 2). We have not collected individuals from this lineage within any protected areas.

Central-southern Chile is listed as one of the 25 biodiversity hotspots for conservation priorities; it harbors a high level of endemism of species and genera (Myers et al., 2000; Arroyo et al., 2006). It is also true that central-southern Chile suffers high anthropogenic impact since the Spanish colonization five centuries ago (Armesto et al., 2001; Smith-Ramírez, 2004; Silva and Saavedra, 2018). Currently, intensive deforestation coupled with grazing by cattle, human-caused fires, and expansion of invasive species are some of the strongest threats for the local biodiversity (Armesto et al., 2001; Smith-Ramírez, 2004; Escobar et al., 2015). This situation rises the necessity of categorizing not only at the species level, but also at the intraspecific lineages, either in a context of recognized subspecies or not. This might be the case of several species, particularly those with wide geographic distribution, that are categorized as of least concern in southern South America. In this genomic era, with the astonishing technological advances at hand, it is possible to individualize intraspecific lineages, exposing their unique genetic legacy to avoid losing the richness of this heritage.

Data Availability Statement

The data presented in the study are deposited in the EMBL-EBI European Variation Archive (EVA) repository, accession numbers PRJEB43642, ERZ1758066, ERZ1758067, and ERZ1758068; Bioproject PSUB014331; Biosamples SAMD00278635–SAMD00278678.

Ethics Statement

The collection of some Chilean specimens was authorized by the Servicio Agrícola y Ganadero (permits 1231/2017, 5611/2013, 2164 and 5165); collection procedures were approved by the Universidad Austral de Chile (permits 04/11 and 296/2017). Other studied specimens were already housed in museum collections.

Author Contributions

LV conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and tables, authored and reviewed drafts of the manscript, and approved the final draft. GD’E conceived and designed the experiments, authored and reviewed drafts of the manscript, and approved the final draft. Both authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer BP declared a past co-authorship with one of the authors, GD’E, to the handling editor.