Genome-wide data support recognition of an additional species of Neotropical river otter (Mammalia, Mustelidae, Lutrinae)

Abstract

Cryptic biodiversity continues to be revealed worldwide, even in apparently well-known groups such as carnivorans. The Neotropical Otter (Lontra longicaudis) presents shape variation in its nose pad, a character that has been used to differentiate species in this group. Based on this, 3 subspecies are recognized: L. l. annectens (Mexico, Central America, and South America west of the Andes), L. l. enudris (Amazon and Orinoco basins), and L. l. longicaudis (Paraná basin and remaining distribution). Previous studies partially supported their distinctness based on mitochondrial DNA markers, morphometrics, and ecological niche modeling. We analyzed genome-wide nuclear markers (ultraconserved elements) of 29 L. longicaudis individuals across the species’ range to assess its population structure. Phylogenomic analysis recovered L. longicaudis as paraphyletic with robust support, with 1 clade comprising samples from Mexico and Colombia (trans-Andean populations) and another encompassing the remaining samples (cis-Andean populations), which grouped with 2 other South American species, L. felina and L. provocax. Principal component and admixture analyses strongly differentiated the 2 main L. longicaudis groups, and distinguished the Amazonian individuals from the remaining cis-Andean samples. Our results support the recognition of trans-Andean populations of L. longicaudis as a distinct otter species, which should be recognized as Lontra annectens.

Resumo

Uma biodiversidade críptica continua a ser revelada em todo o mundo, mesmo em grupos aparentemente bem conhecidos, como os carnívoros. A lontra neotropical (Lontra longicaudis) apresenta variação de formato em sua almofada nasal, característica esta que tem sido utilizada para diferenciar as espécies de lontra. Com base nisso, três subespécies são reconhecidas: L. l. annectens (México, América Central e América do Sul, a oeste dos Andes), L. l. enudris (bacias dos rios Amazonas e Orinoco) e L. l. longicaudis (bacia do rio Paraná e restante de sua distribuição). Estudos anteriores apoiaram parcialmente a distinção das subespécies com base em marcadores de DNA mitocondrial, morfometria e modelagem de nicho ecológico. Neste estudo, nós analisamos marcadores nucleares ao longo de todo o genoma (elementos ultraconservados) de 29 indivíduos de L. longicaudis, distribuídos por toda a sua distribuição, para avaliar a estruturação populacional desta espécie. A análise filogenômica recuperou, com suporte robusto, L. longicaudis como parafilética, com um clado compreendendo exemplares do México e Colômbia (populações transandinas) e um segundo clado englobando os exemplares cisandinos que se agruparam com outras duas espécies sul-americanas, L. felina e L. provocax. As análises de componentes principais e de mistura genética diferenciaram fortemente os dois principais grupos de L. longicaudis e, além disso, distinguiram os indivíduos amazônicos das demais amostras cis-andinas. Nossos resultados apoiam o reconhecimento das populações transandinas de L. longicaudis como uma espécie distinta de lontra, que deve ser reconhecida como Lontra annectens.

Carnivora, genomics, Lontra, species delimitation, taxonomy, ultraconserved elements

Carnivora, delimitação de espécies, elementos ultraconservados, genômica, Lontra, taxonomia

Correct delimitation of species is a very relevant issue because it affects not only taxonomy, but also ecological analyses, conservation assessment, management programs (in situ and ex situ), and public policy (Zachos 2018). The complexity of this problem is exacerbated in the case of allopatry, cryptic species, and/or lineages that have diverged recently (Fujita et al. 2012), for which phenotypic differences may be subtle. In this context, genetic data provide robust approaches to delimit species that are difficult to recognize solely with morphological traits (Baker and Bradley 2006). In recent years, several additional species of mammals have been recognized in the Neotropics (see Burgin et al. 2018 for a review). Genetic data can also be applied at the intraspecific level, allowing the delimitation of evolutionarily significant units (ESUs; Ryder 1986; Moritz 1994), which are historically differentiated populations with a unique set of traits and thus warrant conservation and management actions as distinct entities (Frankham et al. 2002).

For otters (Mammalia, Mustelidae, Lutrinae), the shape of the rhinarium (nose pad around the nostrils) has been suggested to be an effective morphological character to differentiate species (Foster-Turley et al. 1990). This criterion applies well to all 13 extant otter species, except for the Neotropical Otter, Lontra longicaudis, in which this trait is variable throughout the species distribution (van Zyll de Jong 1972; Foster-Turley et al. 1990; Larivière 1999a), raising the question about its value as an informative character for South American otters (van Zyll de Jong 1972). Previous studies placed several otters from the Americas in the L. platensis group, comprising 6 species: L. annectens, L. enudris, L. platensis, L. provocax, L. felina, and L. canadensis (see van Zyll de Jong 1972 for revision). A later revision considered 4 species as valid: L. canadensis, L. felina, L. provocax, and L. longicaudis—the latter including L. annectens, L. enudris, and L. platensis (van Zyll de Jong 1972). These 3 units were then treated as subspecies (van Zyll de Jong 1972): L. l. annectens occurring in Mexico, Central America, and South America west of the Andes; L. l. enudris in the Amazon and Orinoco rivers and adjacent areas east of the Andes; and L. l. longicaudis (including L. platensis) in the Paraná river basin (and implicitly referring to the rest of the L. longicaudis range). This scheme has been subsequently followed by the main mammalian taxonomic sources up to the present (e.g. Larivière 1999a; Wozencraft 2005; Larivière and Jennings 2009; IUCN 2021; Fig. 1A).

Fig. 1.

(A) Lontra longicaudis geographic range, indicating the distribution of the proposed subspecies (adapted from Hernández‑Romero et al. 2018), and the locations of the samples used in this study; regional populations are identified by different dot colors and symbols. (B) Admixture plots based on 30,205 variable sites across the sampled UCE regions. Each individual is represented by a vertical bar (numbered as in panel A and with regions color/symbol-coded by circles on top). Each bar is colored to reflect its population assignment proportions assuming k = 2 (top) and k = 3 (bottom). (C) Maximum likelihood tree inferred from concatenated genome-wide data considering only sites with <40% missing data (40N dataset). Bootstrap proportions >95% are shown at the nodes of the tree. The asterisk indicates the individual represented by a complete genome sequence (see Methods). Numbers indicate the collection site of each sample as depicted in panel A. (D) Principal Component Analysis (PCA) plot of PC1 vs. PC2 for 29 L. longicaudis individuals (numbered as in panel A) using 30,205 variable sites; regional populations are identified by dot colors and symbols, and subspecies by shaded ellipses. The percentage of the variance explained by each PC is indicated on its respective axis.

Previous studies that investigated patterns of population differentiation throughout the range of L. longicaudis were based on analyses of mitochondrial DNA (mtDNA) markers and/or morphological data (Trinca et al. 2007, 2012; Guerrero et al. 2015; Hernández-Romero et al. 2015, 2018; Ruiz-García et al. 2018). Overall, these studies found patterns that partially agree with the proposed subspecies, with 4 geographically structured groups: (i) a trans-Andean group, including individuals from Mexico, Costa Rica, and Colombia, corresponding to L. l. annectens; and 3 cis-Andean groups, formed by samples from (ii) Bolivia, without a formal nomenclatural designation; (iii) Amazonia, corresponding to L. l. enudris; and (iv) eastern South America, corresponding to L. l. longicaudis (Trinca et al. 2012; Hernández-Romero et al. 2018). A considerably deep divergence was observed between L. l. annectens and the other groups, with the number of mutational steps and estimated divergence time being similar to, or even higher than, that between 2 other congeneric species, L. felina and L. provocax (Trinca et al. 2012; Hernández-Romero et al. 2018).

Divergence between the cis- and trans-Andean groups was also identified with cranial morphological characters (Hernández-Romero et al. 2015). A less striking phylogeographic break was also found between L. l. longicaudis and L. l. enudris (Trinca et al. 2007; Hernández-Romero et al. 2018). An analysis that assessed potential geographical barriers among these groups identified 2 barriers with strong support, one in the Andes and another between the Amazon and Paraná river basins (Hernández-Romero et al. 2018). Additionally, the subspecies presented differences in their ecological niches (Hernández-Romero et al. 2018), further supporting the inference that they are evolutionarily distinct.

However, another mtDNA study including a larger number of samples for the northwestern portion of the distribution—encompassing Colombia, Peru, and Bolivia—did not detect genetic differentiation among populations of the Orinoco, Magdalena, and western Amazonia in the Colombian–Brazil–Peruvian border (Ruiz-García et al. 2018). That study only included a single sample from the Paraná River and no samples from the eastern and southern portions of the distribution, limiting its comparison with the other studies targeting this species.

Given the unique features of mtDNA (including lack of recombination, matrilineal inheritance, and reduced effective population size), its use in phylogeographic inference must be complemented with markers drawn from the nuclear genome to provide a comprehensive understanding of organismal patterns of population structure and history (Avise 1994; Hare 2001). Furthermore, L. longicaudis appears to have male-biased dispersal, with females tending to be philopatric (Trinca et al. 2013). Such a dispersal pattern may exacerbate the difference in phylogeographic inferences drawn from mtDNA versus nuclear markers (Toews and Brelsford 2012). For example, in other animal species with male-biased gene flow for which results from nuclear and mitochondrial markers could be compared, the latter yielded more striking subdivisions, which were often not detected when analyzing nuclear data (e.g. Eizirik et al. 2001; Rüppell et al. 2003; Tchaicka et al. 2007; Caparroz et al. 2009; Lopes et al. 2015). Such observations highlight the need to employ nuclear markers as a complement to mtDNA when assessing historical population subdivisions. Furthermore, a study using microsatellite markers revealed genetic structure of L. longicaudis among adjacent Mexican basins (Latorre-Cárdenas et al. 2020) that corresponded with topographical characteristics of the basins (Latorre-Cárdenas et al. 2021).

Thus, considering the results of previous studies, known biogeographic patterns identified for other mammalian species, and limitation of the genetic markers used so far, our goal was to assess L. longicaudis population structure throughout its distribution by applying genome-wide nuclear data. We compared our results with the previously reported subspecies, seeking to clarify their taxonomic status and to delimit ESUs or management units within this species.

Materials and methods

Samples

We selected 29 blood and tissue samples from the 3 known L. longicaudis subspecies, aiming to encompass the full geographic distribution of this species (Fig. 1A; Appendix  I). All samples were collected from animals that were found dead or kept in ex situ institutions, precluding the need for collection permits according to local laws and regulations. We extracted genomic DNA using standard approaches (e.g. phenol/chloroform; Sambrook et al. 1989, or the QIAamp DNA Kit [QIAGEN]), assessed its quality on 1% agarose gels stained with GelRed (Biotium), and roughly quantified the DNA content of each sample by comparing it to a Low DNA Mass Ladder (Invitrogen). The DNA was then used for library construction, target capture of ultraconserved elements (UCEs), and sequencing on a HiSeqX Illumina platform with 2 × 150 bp paired-end reads. The capture experiment employed 5,472 baits targeting 5,060 tetrapod UCEs (Faircloth et al. 2012), including their flanking regions.

Nuclear DNA data

We used the Phyluce package (Faircloth 2016) to remove adapters and the PALEOMIX 1.2.13.2 pipeline (Schubert et al. 2014) to: (i) filter reads, discarding reads shorter than 40 bp and with quality scores lower than 30; (ii) trim reads and filter out PCR duplicates with Picard MarkDuplicates (http://broadinstitute.github.io/picard/); and (iii) map remaining reads against the Lutra lutra reference genome (Mead et al. 2020) with BWA-backtrack, recommended for short sequences (Li and Durbin 2009). We then generated consensus sequences of all genomes using ANGSD 0.921 (Korneliussen et al. 2014) with parameters doFasta = 2, doCounts = 1, explode = 1, setMinDepth = 10, and minMapQ = 20.

For the final alignment, we included whole-genome sequences of L. canadensis, L. felina, L. provocax, and L. longicaudis that were previously reported by our group (de Ferran et al. 2022). We then filtered sites according to the amount of missing data, using trimAl v1.4 (Capella-Gutierrez et al. 2009) and enforcing 3 different levels of stringency: (i) keeping only sites with <40% missing data (referred to as the “40N” data set henceforth); (ii) keeping only sites with <20% missing data (“20N” henceforth); and (iii) keeping only sites with no missing data (“0N” henceforth). For each level of stringency, the remaining sites were concatenated into a single supermatrix. We constructed a maximum likelihood (ML) tree for each supermatrix with RaxML HPC-PTHREADS 8.2 (Stamatakis 2014), using a GTR+GAMMA model of nucleotide substitution and 100 bootstrap replicates to assess nodal support. Lutra lutra was designated as the outgroup to root phylogenetic trees.

mtDNA data

We downloaded sequences of 3 mtDNA gene fragments reported by Trinca et al. (2012) from GenBank (see Table S2 in Trinca et al. 2012 for accession numbers), which included sequences from L. longicaudis (n = 37), L. felina (n = 1), L. provocax (n = 1), L. canadensis (n = 1), and Amblonyx cinereus (n = 1). Fragments were originally amplified and sequenced from the ATPase8 + ATPase6 and ND5 genes, as well as the control region. We extracted orthologs of the same 3 fragments from a mitochondrial genome assembly of L. lutra (GenBank accession LR822067), which is derived from the mLutLut1 chromosome-length genome assembly (Mead et al. 2020). We also downloaded sequences of the NADH5 gene (645 bp) of 2 individuals of L. longicaudis sampled from Mexico (Hernández-Romero et al. 2018; GenBank accessions KY271064 and KY271066).

Sequences were imported into Geneious Prime 2022.1.1 (https://www.geneious.com) and concatenated, as the mitochondrial genome is inherited as a single locus. MAFFT version 7.490 (Katoh and Standley 2013) with default options (algorithm = AUTO, scoring matrix = 200 PAM/k = 2, gap open penalty = 1.53, offset value = 0.123) was used to produce a multiple sequence alignment of the 44 sequences. This alignment was 1,521 bp in length; however, because 25 of the 44 sequences were missing the last 21 bp from the control region fragment, we trimmed these nucleotides, resulting in a final alignment length of 1,500 bp. We then constructed an ML phylogeny using the RaxML version 8.2.11 plugin in Geneious Prime 2022.1.1 with the following settings: nucleotide substitution model = GTR+GAMMA; algorithm = Rapid Bootstrapping and search for best-scoring ML tree; number of starting trees or bootstrap replicates = 1,000; and parsimony random seed = 3. The final ML and bootstrap trees were visualized in Geneious Prime 2022.1.1 and rooted using A. cinereus and L. lutra as outgroups.

Population structure

To characterize patterns of population differentiation, we first estimated genotype likelihood (GL) at each site from mapped reads using the GATK model (McKenna et al. 2010) in ANGSD, keeping only sites that were present in at least 14 individuals using the following parameters: -GL 2, -doMajorMinor 1, -doMaf 1, -SNP_pval 2e-6, -minMapQ 20, -minQ 20, -minInd 14, -minMaf 0.05, -doGlf 2. We then performed a principal component analysis (PCA) using the pcangsd.py script. With the GLs, we also estimated the best-fit population structure and signatures of recent admixture proportions for each individual using NgsAdmix (Skotte et al. 2013), testing for 2 and 3 genetically distinct clusters (k). We chose NgsAdmix because it performs well with medium- and low-coverage data (Skotte et al. 2013).

Results

All analyses performed here led to consistent results regarding the genetic structure of the included samples (Fig. 1B–D). All phylogenomic trees recovered L. longicaudis as paraphyletic with robust support. For the 20N (Supplementary Data SD1) and 40N trees (Fig. 1C), we identified 3 main clades, with the most divergent comprising L. longicaudis samples from Mexico and Colombia (trans-Andean populations). Of the 2 internal sister groups, 1 comprised L. felina and L. provocax, while the other included all cis-Andean populations of L. longicaudis samples. Within the cis-Andean group, Amazonian samples occupied basal positions, suggesting a dispersal pattern starting from the Amazon southwards toward the rest of its range, which should be further investigated in future studies. The paraphyly of Amazonian samples was consistent among our different data sets; however, bootstrap support for the relevant nodes was limited, indicating that additional data will be required to dissect population substructure within the cis-Andean group.

The 0N tree presented similar results to the 20N and 40N trees (Supplementary Data SD2), but the relationship between the L. felina + L. provocax clade and the cis-Andean longicaudis group was unresolved, as the former clade was contained within the latter. This tree presented very low support for most nodes, which is explained by the very small number of variable and parsimony-informative sites present in this data set relative to the 40N and 20N data sets (Table 1).

Table 1.

Summary information of alignments used for the phylogenomic analysis, including size, number of variable sites, singletons, and parsimony-informative sites.

Filter	Alignment size	Variable sites	% variable sites	Singletons	Parsimony-informative sites
40N	2,054,111	35,910	1.75%	26,197	9,713
20N	1,641,359	26,228	1.60%	18,939	7,289
0N	127,478	646	0.51%	477	169

PCA using 30,205 variable sites identified the Mexican + Colombian samples as a distinct cluster along PC1, while the Northeast, South–Southeast, and Central cis-Andean populations were clustered together and distinct from the Amazon cluster along PC2 (Fig. 1D). PC1 (explaining 35.5% of the total variation) showed a striking segregation between the trans-Andean group (L. l. annectens) and the cis-Andean groups (L. l. enudris and L. l. longicaudis), while PC2 (explaining 3.95% of the variation) segregated the Amazon samples from the others. PC3 (explaining 3.79% of the variation) only segregated the 2 northernmost samples (16, 17) from the Northeast region from all remaining individuals (Supplementary Data SD3).

Admixture analyses also differentiated the trans-Andean population (L. l. annectens) from the others with both clustering scenarios (Fig. 1B). With k = 2, the trans-Andean samples formed a distinct group from all other samples. With k = 3, 3 clusters were identified, 1 formed by all the trans-Andean samples, another by the Amazonian plus 2 Northeastern samples (16 and 17), and the third by 1 Northeastern sample (18) plus the Central and South–Southeastern samples. Individuals 16 and 17 are from NE Brazil (see Fig. 1A), suggesting genetic connectivity between the Amazon and northernmost Atlantic coastal region. Two individuals (19 and 20) from the Central region showed admixture between the Amazonian and the Northeast/South–Southeast groups.

Discussion

Our analyses of genome-wide nuclear data consistently supported a strong genetic separation between trans-Andean and cis-Andean Neotropical river otters traditionally assigned to L. longicaudis. Evidence includes: (i) the PCA results, with PC1 strongly differentiating the trans-Andean L. l. annectens from the cis-Andean groups, while PC2 segregated the latter into 2 groups corresponding to samples from the Amazon and the Northeast, Central, and South–Southeast regions of South America (see Fig. 1D); (ii) admixture analysis, in which the trans-Andean population appears as a distinct unit; and (iii) phylogenetic trees, which recovered (with high support) L. longicaudis as a paraphyletic species, with its trans-Andean group as the most divergent clade and its cis-Andean population as the sister group to (L. felina + L. provocax).

Based on these results, we recommend that the trans-Andean group should be considered a valid species, for which the oldest available name is Lutra annectensMajor, 1897. Considering nomenclatural changes that support the use of Lontra for the American river otter species instead of Lutra (van Zyll de Jong 1972), we suggest use of Lontra annectens for the trans-Andean group. The type locality for this species is San Juan, approximately 24 km west of Huigra (River of Tepic, Nayarit, Mexico—literal translation from Larivière 1999a: “San Juan, a unos 24 km al oeste de Huigra,” [Rio de Tepic, Nayarit, Mexico]). Further studies are needed to define the geographic range of this species, but it is described to occur in Mexico, Central America, and South America west of the Andes (van Zyll de Jong 1972).

This conclusion is also supported by previous studies on L. longicaudis, which found substantial mtDNA divergence between these groups and which we have recapitulated in our results (Supplementary Data SD4), identified the Andes as an important barrier, and found relevant differences in their morphology and ecological niches (van Zyll de Jong 1972; Trinca et al. 2012; Hernández-Romero et al. 2015, 2018).

In addition to revealing this species-level partition, our analyses also indicated substructure within the cis-Andean L. longicaudis group. The admixture analysis and the observed phylogeographic patterns agree with previous morphological data (van Zyll de Jong 1972) that suggested that the cis-Andean clade of L. longicaudis presents clinal variation from north to south. Despite this potentially clinal pattern, our results partially support population structure between the Amazon and remaining populations, which must be tested with additional sampling in intermediate areas, and if confirmed, would agree with the previously reported subspecies. These results also agree with presence of a putative geographical barrier between the Amazon and Paraná river basins suggested in previous work (Hernández-Romero et al. 2018). This barrier could be the Brazilian Central Plateau, an important watershed, with rivers located to the north of it (e.g. Xingu, Tapajós, Araguaia, and Tocantins) flowing northward toward the Amazon region, and those in the south (e.g. Paraguay and Paraná) flowing southward (Buckup 2011). Interestingly, L. longicaudis seems to be rare or absent in most of the southern Amazon region and most of the Brazilian Cerrado (Rheingantz et al. 2017), occurring only in large rivers in those areas. However, it is still unclear if this pattern is real or may be due to limited or biased sampling in those regions. Although the Brazilian Central Plateau seems to be a barrier to dispersal, it appears to be permeable to some extent as some dispersal of individuals occurs between the Amazon and the Northeast Atlantic coast. Migration of L. longicaudis individuals among adjacent hydrographic basins in Mexico has been suggested to occur through their central portions (Latorre-Cárdenas et al. 2021). Should the genetic pattern of southern populations detected in this study be confirmed by additional sampling, presence of the Brazilian Central Plateau could provide an explanation for limited gene flow between the 2 cis-Andean subspecies. Thus, considering data presently available, we suggest that, for conservation and management purposes, 2 units should be considered within the cis-Andean population: Amazon (roughly equivalent to L. l. enudris) and Central–Eastern (roughly equivalent to L. l. longicaudis).

Furthermore, previous studies have indicated the presence of a potential additional mtDNA lineage from Bolivia, basal to the clade comprising most of the Amazonian samples (Trinca et al. 2012; Hernández-Romero et al. 2018). Therefore, considering the possibility that this is a distinct population, it is particularly important that future genetic studies target this region for sampling.

An additional aspect pertaining to mtDNA-based analyses is assessment of the extent to which the species-level nuclear partition detected here between L. longicaudis and L. annectens is also observed with mitochondrial markers. Because previous mtDNA-based phylogeographic studies used L. felina, L. provocax, and/or L. canadensis as outgroups, it was not possible to directly compare their results to ours regarding the relationship between these species and L. longicaudis. To determine if the difference between our results and those of previous studies was due to outgroup choice, we performed ML phylogenetic analyses using a concatenated data set comprised of 3 mitochondrial fragments (ATP6/8, ND5, control region) with L. lutra and A. cinereus as outgroups. Our analyses (Supplementary Data SD4) recovered the traditional L. longicaudis (trans- and cis-Andean) as monophyletic with 100% support, and as a sister species to (L. felina + L. provocax). Within L. longicaudis, the trans-Andean and cis-Andean groups formed distinct clades, as previously reported (Hernández-Romero et al. 2018). Different topologies found in the mitochondrial and nuclear data sets may be caused by incomplete lineage sorting, which is more likely to occur in groups with rapid divergence and large effective population sizes (Pamilo and Nei 1988; Maddison 1997). This is indeed the case for L. longicaudis, which was shown to have a large and increasing historical effective population size based on a coalescent rate analysis of the genome of this species (de Ferran et al. 2022). It may also be due to ancient hybridization between the cis- and trans-Andean groups, causing mtDNA introgression and thus mito-nuclear discordance, as reported previously for other mammals (e.g. Hailer et al. 2012; Toews and Brelsford 2012; Trigo et al. 2013; Seixas et al. 2018). Overall, although the mtDNA data do not retrieve paraphyly of L. longicaudis, they nevertheless support the deep division between its trans-Andean and cis-Andean groups, and are thus consistent with the recognition of L. annectens as a distinct species.

The fact that nuclear data recovered L. annectens as a sister species to the (L. longicaudis/L. felina/L. provocax) clade supports and extends the previous hypothesis that Neotropical Lontra species radiated from a single ancestor as part of the Great American Biotic Interchange (GABI; Koepfli and Wayne 1998; Koepfli et al. 2008; Eizirik 2012). The present results indicate a sequential, north-to-south diversification process, beginning with the split between Nearctic L. canadensis and the Neotropical clade. The latter was then split between trans-Andean L. annectens and the South American clade, which was in turn subdivided between tropical L. longicaudis and the southernmost pair L. provocax + L. felina. Such a sequence provides a remarkably clear biogeographic scenario of colonization and speciation events in this group as it occupied the Neotropics during and after the GABI. Interestingly, our results also suggest a potentially north-to-south pattern of diversification within the cis-Andean L. longicaudis clade. Such a pattern can be further investigated with expanded sampling, especially targeting contact areas between the Amazon and adjacent biomes.

The phylogenetic pattern and present-day biogeography of this group suggest a major role for allopatric speciation in this genus. The exception may be divergence between L. felina and L. provocax, which exhibit a partially sympatric distribution (on a broad scale), although they do not occur in the same habitats (Larivière 1999b), leading a previous study to suggest the possibility of parapatric speciation in this case (Vianna et al. 2010). Dissecting the speciation history of these remarkable organisms, including their adaptive divergence as they colonized the Neotropics, is an exciting prospect that should be further empowered in the future by the integration of genomic, morphological, and ecological studies.

Supplementary data

Supplementary data are available at Journal of Mammalogy online.

Supplementary Data SD1.—Maximum likelihood tree inferred from concatenated genome-wide data considering only sites with <20% missing data (20N data set). Bootstrap proportions >95% out of 100 replicates are shown at nodes of the tree. The asterisk indicates the individual represented by a complete genome sequence (see Materials and methods).

Supplementary Data SD2.—Maximum likelihood tree inferred from concatenated genome-wide data considering only sites with no missing data (0N data set). Bootstrap proportions >80% out of 100 replicates are shown at nodes of the tree. The asterisk indicates the individual represented by a complete genome sequence (see Materials and methods).

Supplementary Data SD3.—Principal component analysis (PCA) plot of PC3 versus PC4 for 29 Lontra longicaudis individuals using 30,205 variable sites with populations identified by dot colors and subspecies by ellipses.

Supplementary Data SD4.—Maximum likelihood tree based on 1,500 bp from concatenated sequences of 3 mtDNA markers (ATP8/6, ND5, and control region). Bootstrap proportions out of 1,000 replicates are shown at nodes of the tree. Labels for Lontra longicaudis (e.g. Ll-ANC1) correspond to haplotype designations used by Trinca et al. (2012). The 2 L. longicaudis samples from Mexico (GenBank accessions KY271064 and KY271066) are from the study by Hernández-Romero et al. (2018). The tree is rooted using Amblonyx cinereus and Lutra lutra. The scale bar at the bottom indicates the number of substitutions per site.

Acknowledgments

We thank Paulo H. Ott, Helen Waldemarin, Patrick Colombo, Caroline Zank, Gabriele Volkmer, Lídio França do Nascimento, Enrique Gonzalez, Raquel Velozo, Felipe G. Grazziotin, Carolina C. Cheida, Kátia Cassaro, Letícia Koprovski, Helio Rubens J. Pereira, Ana Lízia Brito, Fernando Rosas, Projeto Pequenos Cetáceos, Miriam Marmontel, Ana Cristina Mendes de Oliveira, Benoit de Thoisy, Acuario de Veracruz, A.C., and Zológico Benito Juárez in Morelia for access to samples; Giovanna M. T. Oliveira and Gabriele Z. Lazzari for laboratory assistance; the Smithsonian Institution High Performance Cluster (SI/HPC) for access to computational resources; and Drs. A. R. Percequillo, F. R. Santos, and S. L. Bonatto for constructive comments on a previous version of this manuscript (part of VdF’s doctoral dissertation at PPG-EEB/PUCRS).

Author contributions

VdF contributed in Conceptualization, Data Curation, Formal Analysis, Funding Acquisition, Investigation, Methodology, Resources, Validation, Visualization, Writing—original draft, Writing—review and editing; HVF contributed in Data Curation, Formal Analysis, Investigation, Methodology, Validation; CST contributed in Conceptualization, Data Curation, Funding Acquisition, Resources; PCH-R contributed in Data Curation, Resources, Writing—review and editing; GPL contributed in Resources, Writing—review and editing; CG-R contributed in Data Curation, Funding Acquisition, Project administration, Resources, Supervision, Writing—review and editing; K-PK contributed in Conceptualization, Formal analysis, Funding acquisition, Methodology, Project administration, Supervision, Validation, Writing—original draft, Writing—review and editing; EE contributed in Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing—original draft, Writing—review and editing.

Funding

This study was supported by CNPq/Brazil (grants 141172/2017-7 [awarded to VdF] and 309068/2019-3 [awarded to EE]), the PUCRS/CAPES-PrInt Program (fellowship 88887.370464/2019-00 awarded to VdF), and the Office of Naval Research Global/USA (award N62909-15-1-N107 to EE). This study is a contribution of the National Institutes for Science and Technology (INCT) in Ecology, Evolution and Biodiversity Conservation, supported by MCTIC/CNPq/Brazil (proc. 465610/2014-5) and FAPEG/Brazil (proc. 201810267000023).

Conflict of interest

None declared.

Data availability

Newly generated raw DNA sequences (in fastq format) have been deposited at NCBI and are publicly available as of the date of publication.

Appendix I

Location of Lontra longicaudis samples used in this study (see Fig. 1A). Region corresponds to the regional sampling localities indicated in the legend of Fig. 1A.

Sample	Latitude	Longitude	Country	Locality	Region
1	18.9874	−96.9850	Mexico	San José Neria, Chocamán, Veracruz	Mexico
2	19.0940	−96.1408	Mexico	Playa de Vacas, Medellín, Veracruz	Mexico
3	19.0692	−96.1390	Mexico	Paso Colorado, Boca del Río, Veracruz	Mexico
4	18.8417	−95.9494	Mexico	Rincón de la Palma, Alvarado, Veracruz	Mexico
5	18.6094	−95.6586	Mexico	Tlacotalpan, Tlacotalpan, Veracruz	Mexico
6	18.2371	−95.8234	Mexico	Paso del Cura, Chacaltianguis, Veracruz	Mexico
7	18.2371	−95.8234	Mexico	Paso del Cura, Chacaltianguis, Veracruz	Mexico
8	17.9569	−94.5511	Mexico	Capoacan, Minatitlan, Veracruz	Mexico
9	6.5506	−74.7873	Colombia	Maceo, Antioquia	Colombia
10	2.8162	−60.6660	Brazil	Boa Vista, Roraima	Amazon
11	5.6680	−53.7789	French Guiana	Maná, Ouest Guyanais	Amazon
12	−0.0407	−51.0789	Brazil	Amapá	Amazon
13	−3.3059	−64.7903	Brazil	Igarapé de Mineroa, Alvarães, Amazonas	Amazon
14	−2.4754	−60.8177	Brazil	Anavilhanas, Amazonas	Amazon
15	−3.8209	−60.3604	Brazil	Careiro, Amazonas	Amazon
16	−5.7098	−35.2792	Brazil	Extremoz Lagoon, Extremoz, Rio Grande do Norte	Northeast
17	−8.0301	−34.9103	Brazil	Recife metropolitan region, Pernambuco	Northeast
18	−14.6088	−39.0530	Brazil	Almada river, Ilhéus, Bahia	Northeast
19	−15.8386	−58.4507	Brazil	Porto Esperidião, Mato Grosso	Central
20	−18.9942	−57.6602	Brazil	Corumbá, Mato Grosso do Sul	Central
21	−20.3290	−46.3713	Brazil	Vargem Bonita, Minas Gerais	Central
22	−19.9245	−43.9352	Brazil	Belo Horizonte region, Minas Gerais	Central
23	−23.0086	−43.3626	Brazil	Rio de Janeiro, Rio de Janeiro	South–Southeast
24	−23.7651	−45.7604	Brazil	Barra do Uma river, Bertioga, São Paulo	South–Southeast
25	−23.6912	−53.9901	Brazil	Ilha Grande National Park, Paraná	South–Southeast
26	−29.9251	−50.2263	Brazil	Osorio, Rio Grande do Sul	South–Southeast
27	−30.1228	−50.6053	Brazil	Capão da Porteira, Rio Grande do Sul	South–Southeast
28	−31.7151	−55.9838	Uruguay	Tacuarembó, Tacuarembó	South–Southeast
29	−31.2114	−60.1624	Argentina	Cayastá, Santa Fé	South–Southeast
*	−2.0579	−60.0257	Brazil	Presidente Figueiredo, Amazonas	Amazon

*Individual represented by a whole-genome sequence (see text for details).