YMPEV 4915
No. of Pages 7, Model 5G
7 June 2014 Molecular Phylogenetics and Evolution xxx (2014) xxx–xxx 1
Contents lists available at ScienceDirect
Molecular Phylogenetics and Evolution journal homepage: www.elsevier.com/locate/ympev 5 6
A primer on the phylogeography of Lagothrix lagotricha (sensu Fooden) in northern South America
3 4 7 8 9 10 12 11 13 1 2 5 9 16 17 18 19 20 21 22 23 24 25 26 27 28
Q1
Sergio Botero a,⇑, Pablo R. Stevenson b, Anthony Di Fiore c a
Laboratory of Cellular Biophysics, The Rockefeller University, New York, NY 10065, USA Departamento de Ciencias Biológicas, Universidad de los Andes, CO-4976 Bogotá, Colombia c Department of Anthropology, University of Texas at Austin, Austin, TX 78712, USA b
a r t i c l e
i n f o
Article history: Received 30 September 2013 Revised 29 April 2014 Accepted 15 May 2014 Available online xxxx Keywords: Woolly monkey Lagothricha Lugens Poeppiggi Haplotype network D-loop
a b s t r a c t The taxonomic history of the genus Lagothrix is complex, with molecular and morphological assessments giving conflicting results for the separation between its taxa. Phylogeographic studies of the most widely distributed species, Lagothrix lagotricha, have only been attempted recently and are limited to few individuals per collection site, many of which were captive making their geographical origin dubious. There is debate regarding the possibility of raising subspecies of Lagothrix lagotricha to the species level, therefore the geographical origin of samples is particularly relevant. In the present work we revisit the intraspecific phylogeography of L. lagotricha from northwestern South America, including the subspecies L. l. poeppiggi, L. l. lagotricha and L. l. lugens (sensu Fooden, 1963), using DNA sequence data from hypervariable region I of the mitochondrial control region (D-loop HVI). Our results suggest a complex picture in which there are well delimited evolutionary units that, nonetheless, do not correlate well with the morphological variation used to support the current delimitation of taxa. Additionally, we corroborate previous results showing a lack of reciprocal monophyly between the putative subspecies of Lagothrix lagotricha, and we propose that this may be due to ancestral polymorphism that has been maintained following the recent spread of woolly monkeys throughout the western Amazonian lowlands and into the inter-Andean region of Colombia. Ó 2014 Elsevier Inc. All rights reserved.
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
48 49
1. Introduction
50
Lagothrix species, commonly known as woolly monkeys, are large Neotropical frugivores from the subfamily Atelinae distributed across the central and western Amazon basins in Brazil, Bolivia, Peru, Ecuador, western Venezuela, and Colombia, extending north along eastern cordillera of the Andes (Fooden, 1963). Wooly monkeys are important as seed dispersers (Stevenson, 2000, 2007; Stevenson and Guzmán-Caro, 2010), but their large size and specialization as frugivores makes them particularly vulnerable to anthropogenic effects, especially hunting (Peres and Palacios, 2007; Stevenson and Aldana, 2008). The taxonomic history of the genus is complex. In the first major taxonomic revision of the woolly monkeys, Fooden (1963) recognized two species within the genus, Lagothrix flavicauda, the yellow tailed woolly monkey and Lagothrix lagotricha, the lowland
51 52 53 54 55 56 57 58 59 60 61 62 63
⇑ Corresponding author. Address: 1230 York Avenue Box 18, New York, NY 10065, USA. E-mail addresses: [email protected] (S. Botero), [email protected]. co (P.R. Stevenson), anthony.difi[email protected] (A. Di Fiore).
or common woolly monkey, with four subspecies. Four decades later, Groves (2001) undertook a comparative cladistic analysis of all of the atelids using craniodental morphology and concluded that the yellow-tailed woolly monkey was most closely related to spider monkeys, not to other woolly monkeys. He thus resurrected a former subgenus name, Oreonax, for O. flavicauda, and raised all of the subspecies of Lagothrix lagotricha to species status: Lagothrix lagotricha, L. lugens, L. cana, and L. poeppigii. More recently, the resurrection of Oreonax has been shown to be an artifact of sampling choice when replicating the methodology that Groves (2001) employed (Matthews and Rosenberger, 2008). Moreover, molecular analyses of both a large mitogenomic dataset (Di Fiore et al., this issue) and a large suite of nuclear markers (Chaves et al., 2012) have shown that Lagothrix flavicauda is indeed the sister group to other woolly monkeys, with a divergence time estimated at roughly 2 Ma, but have not resolved the relationships among the putative Lagothrix lagotricha subspecies. One study using the mitochondrial cytochrome oxidase II gene (COII) found a lack of reciprocal monophyly between the currently-recognized subspecies (Ruiz-Garcia and Pinedo-Castro, 2010), and this was corroborated for Colombian L. l. lugens and
http://dx.doi.org/10.1016/j.ympev.2014.05.019 1055-7903/Ó 2014 Elsevier Inc. All rights reserved.
Please cite this article in press as: Botero, S., et al. A primer on the phylogeography of Lagothrix lagotricha (sensu Fooden) in northern South America. Mol. Phylogenet. Evol. (2014), http://dx.doi.org/10.1016/j.ympev.2014.05.019
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
YMPEV 4915
No. of Pages 7, Model 5G
7 June 2014 2 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
S. Botero et al. / Molecular Phylogenetics and Evolution xxx (2014) xxx–xxx
L. l. lagotricha using the hypervariable region I of the mitochondrial control region (D-loop HVI), COII, and karyotyping (Botero et al., 2010). Thus, there seems to be little rationale to consider the different forms of lowland woolly monkeys as different species, and we will use Fooden’s (1963) subspecific nomenclature for the rest of the paper. A major limitation of the aforementioned studies is that they were based on either a single or few individuals for each collection locale. In many cases, too, the samples used for analysis came from individuals of unconfirmed origin such as pets or animals hunted for subsistence (Ruiz-Garcia and Pinedo-Castro, 2010), or from captive individuals whose geographical provenience was unknown (Botero et al., 2010). The use of samples with unconfirmed origin complicates the interpretation of results since there is significant phenotypical variability among and within each of the putative subspecies and intermediate phenotypes are common (Defler, 2004). The only phylogenetic study to include wild individuals directly sampled by the researchers included six Colombian populations and confirmed the lack of monophyly between L. l. lagotricha and L. l. lugens using the D-loop HVI (Botero and Stevenson, 2014). This study also found a population of woolly monkeys that shows the brown phenotype characteristic of L. l. lagotricha and that is located within that subspecies’ geographical range, but is nonetheless indistinguishable at the molecular level for the D-loop HVI marker from the closest population of L. l. lugens included in the study. When this population is considered to be L. l. lugens, as its mitochondrial D-loop haplotype implies, Colombian woolly monkeys indeed segregate into two well separated evolutionary units. Botero and Stevenson (2014) thus suggest that these molecularly defined demes might be sufficiently different to be considered subspecies, but that the current taxonomic scheme based on pelage coloration is not appropriate to distinguish them. While these results remain to be corroborated with other molecular markers and are limited to L. l. lagotricha and L. l. lugens, they suggest a complex picture for the phylogeography of Lagothrix that warrants further investigation. Our goal in this paper, then, is to expand previous work by revisiting the intraspecific phylogeography of Lagothrix lagotricha from northwestern South America, incorporating samples from L. l. poeppigii, as well as L. l. lagotricha and L. l. lugens into the analysis. Samples of the fourth subspecies, L. l. cana, which has the largest and southernmost distribution were not available for our analysis.
128
2. Materials and methods
129
2.1. Populations and samples included
130
We included previously published genetic samples from 8 wild populations of woolly monkeys, including two populations of L. l. lagotricha, four of L. l. lugens (Botero and Stevenson, 2014) and two of L. l. poeppiggi (Di Fiore and Fleischer, 2005; Di Fiore et al., 2009) (Fig. 1). We only included populations represented by samples from multiple wild individuals where we were certain of the geographic origin of the samples. Each population location fell within the geographical limits of an ascribed subspecies and showed the characteristic phenotype of that subspecies (Fooden, 1963). However, we assigned samples from the Guaviare population in Colombia as L. l. lugens for our analysis, since in a previous study (Botero and Stevenson, 2014) these samples were indistinguishable in terms of their mtDNA haplotypes at the HVI locus from the Macarena, Colombia population of L. l. lugens. In the previous studies sampling was performed non-invasively by collecting fecal matter during extended follows of social groups (Botero and Stevenson, 2014; Di Fiore et al., 2009), except for some
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
samples of L. l. poeppigii, which were collected using tissue biopsy darting (Di Fiore and Fleischer, 2005). All available samples for each population were included, thus the number of samples reflects the success in DNA extraction and sample collection. No new samples were collected for the present study. The coordinates and sampling sizes for each population are indicated in Table 1. Fig. 1 shows their location.
147
2.2. Molecular marker
154
We used a 431 bp fragment of the D-loop HVI as the molecular marker since it was available from the previous study of the L. l. lagotricha, L. l. lugens (Botero and Stevenson, 2014), and L. l. poeppiggi (Di Fiore, 2009) populations. Accession numbers for the sequences in GenBank are given in Table 1. While the L. l. poeppiggi sequences have been used in previous studies (Di Fiore, 2009), they had not been submitted to GenBank previously and conditions for their amplification were not provided in that publication. These sequences were amplified using primers H16340 (50 -CCTGARGTAGGAACCARATG-30 ) and L15926 (50 -SAATTACCCCGGYCTTGTAAACC-30 ) which amplify a fragment of approximately 610 bp (Collins and Dubach, 2000; Kocher et al., 1989). PCR amplification was performed in 25 ll reactions with final concentrations of 0.2 mM each dNTP, 1.5 mM MgCL2, 0.48 lM of each primer, 1 ll of a 1:10 dilution of DNA template, 1 of PCR reaction buffer, and 0.125 ll of Promega GoTaqÒ polymerase at 5 U/ll. Cycling conditions consisted of an initial denaturation step of 5 min at 94 °C followed by 40 cycles of 1 min denaturation at 94 °C, 1 min annealing at 55 °C and 1.5 min extension at 72 °C, after which a final extension step of 5 min at 72 °C and a 64 °C incubation step were used to stop the reaction. In cases where the selected primers did not provide a clear band or failed to amplify, 3 alternate PCR reactions were attempted with primer pairs that amplify shorter fragments within the HVI region: H16340 – R264 (50 -ACAACAAGCTTACAAGCAAGTAC-30 ), F222 (50 -ATGGATTTACGTGTTAGATGGC-30 ) – R423 (50 -AGGCTTTGCCACAAAGTACC-30 ), and F373 (50 -AATGCACTAATTACATAGGG-30 ) – L15926. For the last fragment, the annealing temperature was reduced to 53 °C, and in all cases the reaction was supplemented with BSA to a final concentration of 1 lg/ll. Note that for the present study we only included sequences for which the complete 431 bp region corresponding to fragment sequenced by Botero and Stevenson (2014) was available. Other details about the L. l. poeppiggi sequences have been published elsewhere (Di Fiore, 2009). Traditionally, woolly monkeys have been thought to be characterized by having philopatric males and dispersing females. However, evidence of male dispersal (Di Fiore et al., 2009; Maldonado and Botero, 2009) and great flexibility in the grouping patterns in some populations has been observed (Stevenson et al. unpublished). If only females disperse, mitochondrial estimates would be an accurate estimate of the total gene flow, whereas if males disperse too, the total gene flow is expected to be greater. As a result, estimates based on mitochondrial markers do not fully represent the patterns of gene flow occurring, but a lower bound to it. Our conclusions should be treated with caution until replicated with nuclear markers.
155
2.3. Analyses
201
We used the software DNAsp v5.00.07 (Rozas et al., 2003) to calculate basic descriptive statistics: nucleotide diversity p (Nei and Li, 1979), number of haplotypes and segregating sites per population and per taxa. We then calculated Tajima’s D (Tajima, 1989), Fu’s F (Fu, 1997), and Ramos-Onsins and Rosas’ R2 (Ramos-Onsins and Rozas, 2002) statistics to examine whether, for any of the taxa, the patterns of molecular diversity were suggestive of deviation
202
Please cite this article in press as: Botero, S., et al. A primer on the phylogeography of Lagothrix lagotricha (sensu Fooden) in northern South America. Mol. Phylogenet. Evol. (2014), http://dx.doi.org/10.1016/j.ympev.2014.05.019
148 149 150 151 152 153
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
203 204 205 206 207 208
YMPEV 4915
No. of Pages 7, Model 5G
7 June 2014 3
S. Botero et al. / Molecular Phylogenetics and Evolution xxx (2014) xxx–xxx
Fig. 1. Map indicating the areas where the populations included in the study were sampled. Countries and states limits are visible in grey and light grey respectively. Map modified from the http://geocommons.com/ map creator tool.
Table 1 Populations’ coordinates, sample sizes, descriptive statistics and GenBank accession numbers.
209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228
Population
Taxon
N
Nucleotide diversity
Haplotypes present
Segregating sites
Coordinates
GenBank accession numbers
Yasuní Tiputini Amazonas Vaupés Guaviare Meta Caquetá Huila
L. L. L. L. L. L. L. L.
37 17 16 10 9 17 9 11
0.034 0.033 0.029 0.029 0.007 0.005 0.015 0.025
18 11 11 6 7 4 3 6
61 43 38 27 9 7 13 26
0°420 S, 6°280 W 0°380 S, 76°90 W 3°230 S, 70°90 W 1°40 S, 69°300 W 2°220 N, 72°410 W 2°370 N, 74°40 W 1°200 N, 74°530 W 1°360 N, 76°60 W
KF704183–KF704219 KF704220–KF704236 GU212712–GU212727 GU212774–GU212783 GU212737–GU212745 GU212757–GU212773 GU212728–GU212736 GU212746–GU212756
l. l. l. l. l. l. l. l.
poeppigii poeppigii lagothricha lagothricha lugens lugens lugens lugens
from a neutral evolution model under constant population size. We calculated coalescent significance statistics with 10,000 replicates for these analyses. We performed an analysis of molecular variance (AMOVA) to characterize the population structure (Excoffier et al., 1992). This analysis is analogous to a conventional statistical analysis of variance, but in which the variance is divided among relevant evolutionary units. In this case, it would correspond to variance within populations (Fst), variance between populations within subspecies (Fsc), and variance between subspecies (Fct). We used 10,000 permutations to estimate confidence values. We also estimated pairwise Fst (Holsinger and Weir, 2009) values between all populations with 1000 permutations to calculate their significance. We used the Holm–Bonferroni correction for multiple testing (Holm, 1979) to determine the significance of both the AMOVA and pairwise Fst calculations for an a level of 0.05. We ran a Mantel test comparing the pairwise Fst matrix among all populations to a matrix of linear distances between the populations. For this analysis we used the distances in hundreds of kilometers (100 km) and 10,000 permutations to estimate confidence values. These analyses
were performed in the software Arlequin v3.1 (Excoffier et al., 2005). As a final analysis, we estimated a network for the entire set of D-loop HVI haplotypes recovered in our study using a statistical parsimony approach (Templeton et al., 1992). For this we used the software TCS (Clement et al., 2000) with a fix connection limit of 100 steps, in large excess of the required steps, to guarantee that the network was fully connected.
229
3. Results and discussion
237
Descriptive statistics are presented in Table 1 for all populations included in the study and in Table 2 for populations assigned to the same putative taxa. Overall, L. l. poeppigii and L. l. lagotricha had similar nucleotide diversities, while that of L. l. lugens was about half that of the other putative subspecies, even with the larger number of populations and individuals included. No HVI haplotypes were shared between taxa. Table 2 also shows the results of the neutrality tests for all putative subspecies, which were
238
Please cite this article in press as: Botero, S., et al. A primer on the phylogeography of Lagothrix lagotricha (sensu Fooden) in northern South America. Mol. Phylogenet. Evol. (2014), http://dx.doi.org/10.1016/j.ympev.2014.05.019
230 231 232 233 234 235 236
239 240 241 242 243 244 245
YMPEV 4915
No. of Pages 7, Model 5G
7 June 2014 4
S. Botero et al. / Molecular Phylogenetics and Evolution xxx (2014) xxx–xxx
Table 2 Descriptive statistics and neutrality tests for all putative subspecies.
246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293
Taxon
N
Nucleotide diversity
Haplotypes present
Segregating sites
Fu’s Fs
Ramos-Onsins and Rozas’ R2
Tajima’s D
L. l. poeppigii L. l. lagothricha L. l. lugens
54 26 46
0.034 0.031 0.016
27 17 17
65 43 33
1.76 p = 0.32 1.15 p = 0.33 1.32 p = 0.37
0.11 p = 0.63 0.15 p = 0.81 0.1 p = 0.44
0.04 p = 0.59 0.47 p = 0.74 0.37 p = 0.42
non-significant in all cases, indicating that the null hypothesis of neutral evolution and lack of recent demographic contraction or expansion cannot be rejected. We caution the reader that these tests are interrogating deeper population history and would thus not detect recent anthropogenic effects. The low nucleotide diversity observed in L. l. lugens, even with the higher sample size, contrasts with the results of Ruiz-Garcia and Pinedo-Castro (2010) where L. l. lugens showed the highest nucleotide diversity of all four putative subspecies. Botero et al. (2010), in their study of captive individuals, included both COII and D-loop HVI and found much higher diversity in the later marker, which was again higher for L. l. lagotricha in both cases but with a much smaller difference for COII. Table 3 compares the nucleotide diversities found in these three studies; while there is not enough information to describe a pattern, it is clear that sampling effects can alter the results significantly, thus any results should be interpreted with caution until both the number of samples and the molecular markers evaluated increase. It is also possible that the COII marker used by Ruiz-Garcia and Pinedo-Castro (2010) is not experiencing neutral evolution, given its relevance for aerobic metabolism (Tsukihara et al., 1996), but this was not detected in previous studies. Results for the AMOVA analysis (Excoffier et al., 1992) are shown in Table 4. All indexes calculated were significant, reflecting significant structuring of Lagothrix mtDNA diversity at the within population, within subspecies, and between subspecies levels. The estimated value for Fsc was small, 0.09, but the fact that it is significant suggests that there is some level of structure present for populations within taxa. Overall, the AMOVA results indicate that most of the differentiation observed is due to differentiation among taxa. Population pairwise Fst values are shown in Table 5. Estimated Fst values were significantly different from 0 in all comparison except among 5 population pairs: Guaviare-Meta, Huila-Caquetá, Huila-Guaviare, Vaupés-Amazonas and Tiputini-Yasuní. This result implies that, on a genetic basis, those pairs of populations are indistinguishable from each other. For each of these cases, both members of the population pair belong to the same putative subspecies, suggesting a high level of mitochondrial gene flow among populations that are considered to belong to the same taxon, which is in agreement with the results of the AMOVA analysis. Our Mantel test estimated a correlation of 0.52 between the matrix of pairwise Fst values and the corresponding matrix of geographical distance (p = 0.007). Geographic distance explained 27% of the variance in genetic differentiation among the populations sampled. Still, more than two thirds of the observed differentiation among populations is due to processes other than isolation by distance, and, at present, there is no clear evidence suggesting
what these other processes might be. If geographical barriers were present, the relevant distances to include in the distance matrix would be different and by including this correction potentially more variance would be explained. However, the forests woolly monkeys inhabit are presumed to have been continuous until recent anthropogenic intervention (Etter et al., 2006). The pairwise Fst calculations suggest a possible mechanism, since the Caquetá population is differentiated from other populations of its subspecies except the nearby Huila population. The pattern observed would be better explained with a stepping stone model, with an as-yet-to-be identified barrier to gene flow existing between the Amazonian (Caquetá) and Orinoquian (Meta and Guaviare) populations, but all of which are still connected through the Andean populations (Huila). Major rivers could represent one of these barriers to gene flow (Ayres and Clutton-Brock, 1992), but a better sampling scheme is required to adequately address this possibility. The haplotype network among the full set of mtDNA haplotypes sampled across putative subspecies is shown in Fig. 2. The different taxa seem to be well delimited into clusters and intermixed within, as the AMOVA analysis suggests. However, as expected from the lack of reciprocal monophyly among putative subspecies observed in previous studies (Botero et al., 2010; Botero and Stevenson, 2014; Ruiz-Garcia and Pinedo-Castro, 2010), there are a few haplotypes that break this pattern. Interestingly, the haplotypes that break the monophyly tend to be several mutations away. This result supports the idea that paraphyly is due to maintenance of ancestral polymorphisms; if a dispersion event was responsible for this pattern one would expect the paraphyletic haplotypes to be more closely related to those of their population of origin. L. l. lugens showed the lowest level of haplotype diversity even when it was represented in our study by the most populations sampled, while L. l. lagotricha and L. l. poeppigii showed overall higher diversity and greater distance among haplotypes. This could suggest a more recent origin for L. l. lugens, but the data does not allow for more precise analysis. It is worth noting that we attempted to evaluate the levels of gene flow using an isolation with migration model, which is ideal for the problem at hand (Hey, 2010), but the data does not contain enough information and we were unable to converge on any estimator. This analysis worked previously in the two taxon case (Botero and Stevenson, 2014), but we would have had to make an assumption about the topology of the relationships among the 3 taxa to constrain the model and be able to converge. Given that such topology is part of the problem at hand we consider such an approach inadequate. It is important to note that the Guaviare and Meta populations cluster together and share haplotypes, as expected from the fact that the Fst between these populations was low and non-significant.
Table 3 Nucleotide diversity (p) comparison between present and previous studies. Study
Ruiz-Garcia and Pinedo-Castro (2010)
Botero et al. (2010)
Taxon n Marker
COII
COII
D-loop HVI
This study D-loop HVI
L. l. lagothricha L. l. lugens L. l. poeppigii
0.005 0.041 0.026
0.01 0.008
0.042 0.018
0.031 0.016 0.034
Please cite this article in press as: Botero, S., et al. A primer on the phylogeography of Lagothrix lagotricha (sensu Fooden) in northern South America. Mol. Phylogenet. Evol. (2014), http://dx.doi.org/10.1016/j.ympev.2014.05.019
294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341
YMPEV 4915
No. of Pages 7, Model 5G
7 June 2014 5
S. Botero et al. / Molecular Phylogenetics and Evolution xxx (2014) xxx–xxx Table 4 Fixation indices obtained in the AMOVA analysis. Hierarchical level
Index
Structure within populations Structure among taxa Structure among populations within taxa
Fst Fct Fsc
0.65 0.61 0.09
Genetic variation explained (%)
p Value
34.9 61.5 3.6
No. of Pages 7, Model 5G
7 June 2014 Molecular Phylogenetics and Evolution xxx (2014) xxx–xxx 1
Contents lists available at ScienceDirect
Molecular Phylogenetics and Evolution journal homepage: www.elsevier.com/locate/ympev 5 6
A primer on the phylogeography of Lagothrix lagotricha (sensu Fooden) in northern South America
3 4 7 8 9 10 12 11 13 1 2 5 9 16 17 18 19 20 21 22 23 24 25 26 27 28
Q1
Sergio Botero a,⇑, Pablo R. Stevenson b, Anthony Di Fiore c a
Laboratory of Cellular Biophysics, The Rockefeller University, New York, NY 10065, USA Departamento de Ciencias Biológicas, Universidad de los Andes, CO-4976 Bogotá, Colombia c Department of Anthropology, University of Texas at Austin, Austin, TX 78712, USA b
a r t i c l e
i n f o
Article history: Received 30 September 2013 Revised 29 April 2014 Accepted 15 May 2014 Available online xxxx Keywords: Woolly monkey Lagothricha Lugens Poeppiggi Haplotype network D-loop
a b s t r a c t The taxonomic history of the genus Lagothrix is complex, with molecular and morphological assessments giving conflicting results for the separation between its taxa. Phylogeographic studies of the most widely distributed species, Lagothrix lagotricha, have only been attempted recently and are limited to few individuals per collection site, many of which were captive making their geographical origin dubious. There is debate regarding the possibility of raising subspecies of Lagothrix lagotricha to the species level, therefore the geographical origin of samples is particularly relevant. In the present work we revisit the intraspecific phylogeography of L. lagotricha from northwestern South America, including the subspecies L. l. poeppiggi, L. l. lagotricha and L. l. lugens (sensu Fooden, 1963), using DNA sequence data from hypervariable region I of the mitochondrial control region (D-loop HVI). Our results suggest a complex picture in which there are well delimited evolutionary units that, nonetheless, do not correlate well with the morphological variation used to support the current delimitation of taxa. Additionally, we corroborate previous results showing a lack of reciprocal monophyly between the putative subspecies of Lagothrix lagotricha, and we propose that this may be due to ancestral polymorphism that has been maintained following the recent spread of woolly monkeys throughout the western Amazonian lowlands and into the inter-Andean region of Colombia. Ó 2014 Elsevier Inc. All rights reserved.
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
48 49
1. Introduction
50
Lagothrix species, commonly known as woolly monkeys, are large Neotropical frugivores from the subfamily Atelinae distributed across the central and western Amazon basins in Brazil, Bolivia, Peru, Ecuador, western Venezuela, and Colombia, extending north along eastern cordillera of the Andes (Fooden, 1963). Wooly monkeys are important as seed dispersers (Stevenson, 2000, 2007; Stevenson and Guzmán-Caro, 2010), but their large size and specialization as frugivores makes them particularly vulnerable to anthropogenic effects, especially hunting (Peres and Palacios, 2007; Stevenson and Aldana, 2008). The taxonomic history of the genus is complex. In the first major taxonomic revision of the woolly monkeys, Fooden (1963) recognized two species within the genus, Lagothrix flavicauda, the yellow tailed woolly monkey and Lagothrix lagotricha, the lowland
51 52 53 54 55 56 57 58 59 60 61 62 63
⇑ Corresponding author. Address: 1230 York Avenue Box 18, New York, NY 10065, USA. E-mail addresses: [email protected] (S. Botero), [email protected]. co (P.R. Stevenson), anthony.difi[email protected] (A. Di Fiore).
or common woolly monkey, with four subspecies. Four decades later, Groves (2001) undertook a comparative cladistic analysis of all of the atelids using craniodental morphology and concluded that the yellow-tailed woolly monkey was most closely related to spider monkeys, not to other woolly monkeys. He thus resurrected a former subgenus name, Oreonax, for O. flavicauda, and raised all of the subspecies of Lagothrix lagotricha to species status: Lagothrix lagotricha, L. lugens, L. cana, and L. poeppigii. More recently, the resurrection of Oreonax has been shown to be an artifact of sampling choice when replicating the methodology that Groves (2001) employed (Matthews and Rosenberger, 2008). Moreover, molecular analyses of both a large mitogenomic dataset (Di Fiore et al., this issue) and a large suite of nuclear markers (Chaves et al., 2012) have shown that Lagothrix flavicauda is indeed the sister group to other woolly monkeys, with a divergence time estimated at roughly 2 Ma, but have not resolved the relationships among the putative Lagothrix lagotricha subspecies. One study using the mitochondrial cytochrome oxidase II gene (COII) found a lack of reciprocal monophyly between the currently-recognized subspecies (Ruiz-Garcia and Pinedo-Castro, 2010), and this was corroborated for Colombian L. l. lugens and
http://dx.doi.org/10.1016/j.ympev.2014.05.019 1055-7903/Ó 2014 Elsevier Inc. All rights reserved.
Please cite this article in press as: Botero, S., et al. A primer on the phylogeography of Lagothrix lagotricha (sensu Fooden) in northern South America. Mol. Phylogenet. Evol. (2014), http://dx.doi.org/10.1016/j.ympev.2014.05.019
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
YMPEV 4915
No. of Pages 7, Model 5G
7 June 2014 2 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
S. Botero et al. / Molecular Phylogenetics and Evolution xxx (2014) xxx–xxx
L. l. lagotricha using the hypervariable region I of the mitochondrial control region (D-loop HVI), COII, and karyotyping (Botero et al., 2010). Thus, there seems to be little rationale to consider the different forms of lowland woolly monkeys as different species, and we will use Fooden’s (1963) subspecific nomenclature for the rest of the paper. A major limitation of the aforementioned studies is that they were based on either a single or few individuals for each collection locale. In many cases, too, the samples used for analysis came from individuals of unconfirmed origin such as pets or animals hunted for subsistence (Ruiz-Garcia and Pinedo-Castro, 2010), or from captive individuals whose geographical provenience was unknown (Botero et al., 2010). The use of samples with unconfirmed origin complicates the interpretation of results since there is significant phenotypical variability among and within each of the putative subspecies and intermediate phenotypes are common (Defler, 2004). The only phylogenetic study to include wild individuals directly sampled by the researchers included six Colombian populations and confirmed the lack of monophyly between L. l. lagotricha and L. l. lugens using the D-loop HVI (Botero and Stevenson, 2014). This study also found a population of woolly monkeys that shows the brown phenotype characteristic of L. l. lagotricha and that is located within that subspecies’ geographical range, but is nonetheless indistinguishable at the molecular level for the D-loop HVI marker from the closest population of L. l. lugens included in the study. When this population is considered to be L. l. lugens, as its mitochondrial D-loop haplotype implies, Colombian woolly monkeys indeed segregate into two well separated evolutionary units. Botero and Stevenson (2014) thus suggest that these molecularly defined demes might be sufficiently different to be considered subspecies, but that the current taxonomic scheme based on pelage coloration is not appropriate to distinguish them. While these results remain to be corroborated with other molecular markers and are limited to L. l. lagotricha and L. l. lugens, they suggest a complex picture for the phylogeography of Lagothrix that warrants further investigation. Our goal in this paper, then, is to expand previous work by revisiting the intraspecific phylogeography of Lagothrix lagotricha from northwestern South America, incorporating samples from L. l. poeppigii, as well as L. l. lagotricha and L. l. lugens into the analysis. Samples of the fourth subspecies, L. l. cana, which has the largest and southernmost distribution were not available for our analysis.
128
2. Materials and methods
129
2.1. Populations and samples included
130
We included previously published genetic samples from 8 wild populations of woolly monkeys, including two populations of L. l. lagotricha, four of L. l. lugens (Botero and Stevenson, 2014) and two of L. l. poeppiggi (Di Fiore and Fleischer, 2005; Di Fiore et al., 2009) (Fig. 1). We only included populations represented by samples from multiple wild individuals where we were certain of the geographic origin of the samples. Each population location fell within the geographical limits of an ascribed subspecies and showed the characteristic phenotype of that subspecies (Fooden, 1963). However, we assigned samples from the Guaviare population in Colombia as L. l. lugens for our analysis, since in a previous study (Botero and Stevenson, 2014) these samples were indistinguishable in terms of their mtDNA haplotypes at the HVI locus from the Macarena, Colombia population of L. l. lugens. In the previous studies sampling was performed non-invasively by collecting fecal matter during extended follows of social groups (Botero and Stevenson, 2014; Di Fiore et al., 2009), except for some
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
samples of L. l. poeppigii, which were collected using tissue biopsy darting (Di Fiore and Fleischer, 2005). All available samples for each population were included, thus the number of samples reflects the success in DNA extraction and sample collection. No new samples were collected for the present study. The coordinates and sampling sizes for each population are indicated in Table 1. Fig. 1 shows their location.
147
2.2. Molecular marker
154
We used a 431 bp fragment of the D-loop HVI as the molecular marker since it was available from the previous study of the L. l. lagotricha, L. l. lugens (Botero and Stevenson, 2014), and L. l. poeppiggi (Di Fiore, 2009) populations. Accession numbers for the sequences in GenBank are given in Table 1. While the L. l. poeppiggi sequences have been used in previous studies (Di Fiore, 2009), they had not been submitted to GenBank previously and conditions for their amplification were not provided in that publication. These sequences were amplified using primers H16340 (50 -CCTGARGTAGGAACCARATG-30 ) and L15926 (50 -SAATTACCCCGGYCTTGTAAACC-30 ) which amplify a fragment of approximately 610 bp (Collins and Dubach, 2000; Kocher et al., 1989). PCR amplification was performed in 25 ll reactions with final concentrations of 0.2 mM each dNTP, 1.5 mM MgCL2, 0.48 lM of each primer, 1 ll of a 1:10 dilution of DNA template, 1 of PCR reaction buffer, and 0.125 ll of Promega GoTaqÒ polymerase at 5 U/ll. Cycling conditions consisted of an initial denaturation step of 5 min at 94 °C followed by 40 cycles of 1 min denaturation at 94 °C, 1 min annealing at 55 °C and 1.5 min extension at 72 °C, after which a final extension step of 5 min at 72 °C and a 64 °C incubation step were used to stop the reaction. In cases where the selected primers did not provide a clear band or failed to amplify, 3 alternate PCR reactions were attempted with primer pairs that amplify shorter fragments within the HVI region: H16340 – R264 (50 -ACAACAAGCTTACAAGCAAGTAC-30 ), F222 (50 -ATGGATTTACGTGTTAGATGGC-30 ) – R423 (50 -AGGCTTTGCCACAAAGTACC-30 ), and F373 (50 -AATGCACTAATTACATAGGG-30 ) – L15926. For the last fragment, the annealing temperature was reduced to 53 °C, and in all cases the reaction was supplemented with BSA to a final concentration of 1 lg/ll. Note that for the present study we only included sequences for which the complete 431 bp region corresponding to fragment sequenced by Botero and Stevenson (2014) was available. Other details about the L. l. poeppiggi sequences have been published elsewhere (Di Fiore, 2009). Traditionally, woolly monkeys have been thought to be characterized by having philopatric males and dispersing females. However, evidence of male dispersal (Di Fiore et al., 2009; Maldonado and Botero, 2009) and great flexibility in the grouping patterns in some populations has been observed (Stevenson et al. unpublished). If only females disperse, mitochondrial estimates would be an accurate estimate of the total gene flow, whereas if males disperse too, the total gene flow is expected to be greater. As a result, estimates based on mitochondrial markers do not fully represent the patterns of gene flow occurring, but a lower bound to it. Our conclusions should be treated with caution until replicated with nuclear markers.
155
2.3. Analyses
201
We used the software DNAsp v5.00.07 (Rozas et al., 2003) to calculate basic descriptive statistics: nucleotide diversity p (Nei and Li, 1979), number of haplotypes and segregating sites per population and per taxa. We then calculated Tajima’s D (Tajima, 1989), Fu’s F (Fu, 1997), and Ramos-Onsins and Rosas’ R2 (Ramos-Onsins and Rozas, 2002) statistics to examine whether, for any of the taxa, the patterns of molecular diversity were suggestive of deviation
202
Please cite this article in press as: Botero, S., et al. A primer on the phylogeography of Lagothrix lagotricha (sensu Fooden) in northern South America. Mol. Phylogenet. Evol. (2014), http://dx.doi.org/10.1016/j.ympev.2014.05.019
148 149 150 151 152 153
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
203 204 205 206 207 208
YMPEV 4915
No. of Pages 7, Model 5G
7 June 2014 3
S. Botero et al. / Molecular Phylogenetics and Evolution xxx (2014) xxx–xxx
Fig. 1. Map indicating the areas where the populations included in the study were sampled. Countries and states limits are visible in grey and light grey respectively. Map modified from the http://geocommons.com/ map creator tool.
Table 1 Populations’ coordinates, sample sizes, descriptive statistics and GenBank accession numbers.
209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228
Population
Taxon
N
Nucleotide diversity
Haplotypes present
Segregating sites
Coordinates
GenBank accession numbers
Yasuní Tiputini Amazonas Vaupés Guaviare Meta Caquetá Huila
L. L. L. L. L. L. L. L.
37 17 16 10 9 17 9 11
0.034 0.033 0.029 0.029 0.007 0.005 0.015 0.025
18 11 11 6 7 4 3 6
61 43 38 27 9 7 13 26
0°420 S, 6°280 W 0°380 S, 76°90 W 3°230 S, 70°90 W 1°40 S, 69°300 W 2°220 N, 72°410 W 2°370 N, 74°40 W 1°200 N, 74°530 W 1°360 N, 76°60 W
KF704183–KF704219 KF704220–KF704236 GU212712–GU212727 GU212774–GU212783 GU212737–GU212745 GU212757–GU212773 GU212728–GU212736 GU212746–GU212756
l. l. l. l. l. l. l. l.
poeppigii poeppigii lagothricha lagothricha lugens lugens lugens lugens
from a neutral evolution model under constant population size. We calculated coalescent significance statistics with 10,000 replicates for these analyses. We performed an analysis of molecular variance (AMOVA) to characterize the population structure (Excoffier et al., 1992). This analysis is analogous to a conventional statistical analysis of variance, but in which the variance is divided among relevant evolutionary units. In this case, it would correspond to variance within populations (Fst), variance between populations within subspecies (Fsc), and variance between subspecies (Fct). We used 10,000 permutations to estimate confidence values. We also estimated pairwise Fst (Holsinger and Weir, 2009) values between all populations with 1000 permutations to calculate their significance. We used the Holm–Bonferroni correction for multiple testing (Holm, 1979) to determine the significance of both the AMOVA and pairwise Fst calculations for an a level of 0.05. We ran a Mantel test comparing the pairwise Fst matrix among all populations to a matrix of linear distances between the populations. For this analysis we used the distances in hundreds of kilometers (100 km) and 10,000 permutations to estimate confidence values. These analyses
were performed in the software Arlequin v3.1 (Excoffier et al., 2005). As a final analysis, we estimated a network for the entire set of D-loop HVI haplotypes recovered in our study using a statistical parsimony approach (Templeton et al., 1992). For this we used the software TCS (Clement et al., 2000) with a fix connection limit of 100 steps, in large excess of the required steps, to guarantee that the network was fully connected.
229
3. Results and discussion
237
Descriptive statistics are presented in Table 1 for all populations included in the study and in Table 2 for populations assigned to the same putative taxa. Overall, L. l. poeppigii and L. l. lagotricha had similar nucleotide diversities, while that of L. l. lugens was about half that of the other putative subspecies, even with the larger number of populations and individuals included. No HVI haplotypes were shared between taxa. Table 2 also shows the results of the neutrality tests for all putative subspecies, which were
238
Please cite this article in press as: Botero, S., et al. A primer on the phylogeography of Lagothrix lagotricha (sensu Fooden) in northern South America. Mol. Phylogenet. Evol. (2014), http://dx.doi.org/10.1016/j.ympev.2014.05.019
230 231 232 233 234 235 236
239 240 241 242 243 244 245
YMPEV 4915
No. of Pages 7, Model 5G
7 June 2014 4
S. Botero et al. / Molecular Phylogenetics and Evolution xxx (2014) xxx–xxx
Table 2 Descriptive statistics and neutrality tests for all putative subspecies.
246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293
Taxon
N
Nucleotide diversity
Haplotypes present
Segregating sites
Fu’s Fs
Ramos-Onsins and Rozas’ R2
Tajima’s D
L. l. poeppigii L. l. lagothricha L. l. lugens
54 26 46
0.034 0.031 0.016
27 17 17
65 43 33
1.76 p = 0.32 1.15 p = 0.33 1.32 p = 0.37
0.11 p = 0.63 0.15 p = 0.81 0.1 p = 0.44
0.04 p = 0.59 0.47 p = 0.74 0.37 p = 0.42
non-significant in all cases, indicating that the null hypothesis of neutral evolution and lack of recent demographic contraction or expansion cannot be rejected. We caution the reader that these tests are interrogating deeper population history and would thus not detect recent anthropogenic effects. The low nucleotide diversity observed in L. l. lugens, even with the higher sample size, contrasts with the results of Ruiz-Garcia and Pinedo-Castro (2010) where L. l. lugens showed the highest nucleotide diversity of all four putative subspecies. Botero et al. (2010), in their study of captive individuals, included both COII and D-loop HVI and found much higher diversity in the later marker, which was again higher for L. l. lagotricha in both cases but with a much smaller difference for COII. Table 3 compares the nucleotide diversities found in these three studies; while there is not enough information to describe a pattern, it is clear that sampling effects can alter the results significantly, thus any results should be interpreted with caution until both the number of samples and the molecular markers evaluated increase. It is also possible that the COII marker used by Ruiz-Garcia and Pinedo-Castro (2010) is not experiencing neutral evolution, given its relevance for aerobic metabolism (Tsukihara et al., 1996), but this was not detected in previous studies. Results for the AMOVA analysis (Excoffier et al., 1992) are shown in Table 4. All indexes calculated were significant, reflecting significant structuring of Lagothrix mtDNA diversity at the within population, within subspecies, and between subspecies levels. The estimated value for Fsc was small, 0.09, but the fact that it is significant suggests that there is some level of structure present for populations within taxa. Overall, the AMOVA results indicate that most of the differentiation observed is due to differentiation among taxa. Population pairwise Fst values are shown in Table 5. Estimated Fst values were significantly different from 0 in all comparison except among 5 population pairs: Guaviare-Meta, Huila-Caquetá, Huila-Guaviare, Vaupés-Amazonas and Tiputini-Yasuní. This result implies that, on a genetic basis, those pairs of populations are indistinguishable from each other. For each of these cases, both members of the population pair belong to the same putative subspecies, suggesting a high level of mitochondrial gene flow among populations that are considered to belong to the same taxon, which is in agreement with the results of the AMOVA analysis. Our Mantel test estimated a correlation of 0.52 between the matrix of pairwise Fst values and the corresponding matrix of geographical distance (p = 0.007). Geographic distance explained 27% of the variance in genetic differentiation among the populations sampled. Still, more than two thirds of the observed differentiation among populations is due to processes other than isolation by distance, and, at present, there is no clear evidence suggesting
what these other processes might be. If geographical barriers were present, the relevant distances to include in the distance matrix would be different and by including this correction potentially more variance would be explained. However, the forests woolly monkeys inhabit are presumed to have been continuous until recent anthropogenic intervention (Etter et al., 2006). The pairwise Fst calculations suggest a possible mechanism, since the Caquetá population is differentiated from other populations of its subspecies except the nearby Huila population. The pattern observed would be better explained with a stepping stone model, with an as-yet-to-be identified barrier to gene flow existing between the Amazonian (Caquetá) and Orinoquian (Meta and Guaviare) populations, but all of which are still connected through the Andean populations (Huila). Major rivers could represent one of these barriers to gene flow (Ayres and Clutton-Brock, 1992), but a better sampling scheme is required to adequately address this possibility. The haplotype network among the full set of mtDNA haplotypes sampled across putative subspecies is shown in Fig. 2. The different taxa seem to be well delimited into clusters and intermixed within, as the AMOVA analysis suggests. However, as expected from the lack of reciprocal monophyly among putative subspecies observed in previous studies (Botero et al., 2010; Botero and Stevenson, 2014; Ruiz-Garcia and Pinedo-Castro, 2010), there are a few haplotypes that break this pattern. Interestingly, the haplotypes that break the monophyly tend to be several mutations away. This result supports the idea that paraphyly is due to maintenance of ancestral polymorphisms; if a dispersion event was responsible for this pattern one would expect the paraphyletic haplotypes to be more closely related to those of their population of origin. L. l. lugens showed the lowest level of haplotype diversity even when it was represented in our study by the most populations sampled, while L. l. lagotricha and L. l. poeppigii showed overall higher diversity and greater distance among haplotypes. This could suggest a more recent origin for L. l. lugens, but the data does not allow for more precise analysis. It is worth noting that we attempted to evaluate the levels of gene flow using an isolation with migration model, which is ideal for the problem at hand (Hey, 2010), but the data does not contain enough information and we were unable to converge on any estimator. This analysis worked previously in the two taxon case (Botero and Stevenson, 2014), but we would have had to make an assumption about the topology of the relationships among the 3 taxa to constrain the model and be able to converge. Given that such topology is part of the problem at hand we consider such an approach inadequate. It is important to note that the Guaviare and Meta populations cluster together and share haplotypes, as expected from the fact that the Fst between these populations was low and non-significant.
Table 3 Nucleotide diversity (p) comparison between present and previous studies. Study
Ruiz-Garcia and Pinedo-Castro (2010)
Botero et al. (2010)
Taxon n Marker
COII
COII
D-loop HVI
This study D-loop HVI
L. l. lagothricha L. l. lugens L. l. poeppigii
0.005 0.041 0.026
0.01 0.008
0.042 0.018
0.031 0.016 0.034
Please cite this article in press as: Botero, S., et al. A primer on the phylogeography of Lagothrix lagotricha (sensu Fooden) in northern South America. Mol. Phylogenet. Evol. (2014), http://dx.doi.org/10.1016/j.ympev.2014.05.019
294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341
YMPEV 4915
No. of Pages 7, Model 5G
7 June 2014 5
S. Botero et al. / Molecular Phylogenetics and Evolution xxx (2014) xxx–xxx Table 4 Fixation indices obtained in the AMOVA analysis. Hierarchical level
Index
Structure within populations Structure among taxa Structure among populations within taxa
Fst Fct Fsc
0.65 0.61 0.09
Genetic variation explained (%)
p Value
34.9 61.5 3.6
