Chromosomal diversity and molecular divergence among three undescribed species of Neacomys (Rodentia, Sigmodontinae) separated by Amazonian rivers

The Neacomys genus (Rodentia, Sigmodontinae) is distributed in the Amazon region, with some species limited to a single endemic area, while others may occur more widely. The number of species within the genus and their geographical boundaries are not known accurately, due to their high genetic diversity and difficulties in taxonomic identification. In this work we collected Neacomys specimens from both banks of the Tapajós River in eastern Amazon, and studied them using chromosome painting with whole chromosome probes of Hylaeamys megacephalus (HME; Rodentia, Sigmodontinae), and molecular analysis using haplotypes of mitochondrial genes COI and Cytb. Chromosome painting shows that Neacomys sp. A (NSP-A, 2n = 58/FN = 68) and Neacomys sp. B (NSP-B, 2n = 54/FN = 66) differ by 11 fusion/fission events, one translocation, four pericentric inversions and four heterochromatin amplification events. Using haplotypes of the concatenated mitochondrial genes COI and Cyt b, Neacomys sp. (2n = 58/FN = 64 and 70) shows a mean divergence of 6.2% for Neacomys sp. A and 9.1% for Neacomys sp. B, while Neacomys sp. A and Neacomys sp. B presents a medium nucleotide divergence of 7.4%. Comparisons were made with other published Neacomys data. The Tapajós and Xingu Rivers act as geographic barriers that define the distribution of these Neacomys species. Furthermore, our HME probes reveal four synapomorphies for the Neacomys genus (associations HME 20/[13,22]/4, 6a/21, [9,10]/7b/[9,10] and 12/[16,17]) and demonstrate ancestral traits of the Oryzomyini tribe (HME 8a and 8b, 18 and 25) and Sigmodontinae subfamily (HME 15 and 24), which can be used as taxonomic markers for these groups.


Introduction
The Amazon is one of the richest biomes in terms of Brazilian mammalian species [1]. However, this biodiversity is not homogeneously distributed. Theories that consider ecologic, morphologic, chromosomal and/or molecular analysis performed with terrestrial vertebrate taxonomic groups have shown the occurrence of distinct biogeographic regions in the Amazon [2][3][4].
The large Amazonian rivers are suggested as major geographic barriers to species distribution in the region [3,5,6]. The Riverine Barrier Hypothesis was first proposed by Wallace [7] and reviewed by many authors [8,9]. Silva et al. [5] recognize eight distinct endemic areas for Amazonian species limited by the major rivers: Belém, Guiana, Imeri, Inambari, Napo, Rondônia, Tapajós and Xingu. Each one of them has a distinct evolutionary history, with regard to species diversification.
Taxonomic studies in the Amazonian regions have been difficult because of overlapping characteristics between distinct species [10]. The taxonomy issue of rodents from the Sigmodontinae subfamily (Rodentia, Cricetidae) has been a problem, since this subfamily belongs to one of the most complex and diverse group of New World mammals [11,12].
Recently, genetic strategies have helped to solve problems related to evolution and taxonomy, such as the comparative analysis of mitochondrial gene sequences of Cytochrome C Oxidase-subunit I (COI) and Cytochrome b (Cytb), frequently used for the comparison of species in the same genus or the same family [13]. Also, chromosome painting has been useful for karyotypic evolution studies based on cross-species chromosome homology [14,15], but only 24 species from seven genera among the Sigmodontinae have been analyzed by this technique [16][17][18][19][20][21][22][23]. Neacomys has been one genus in which the understanding of karyotype evolution is complicated by the number of species, geographic boundaries and phylogenetic relationships.
In order to test the Riverine Barrier Hypothesis, the present study compared Neacomys from different banks of Amazon rivers. We defined the karyotypes of two undescribed species of Neacomys, Neacomys sp. A (NSP-A, 2n = 58/FN = 68) and Neacomys sp. B (NSP-B, 2n = 54/ FN = 66), collected from the right and left banks respectively, of the Tapajós River in Eastern Amazon. Whole chromosome probes of Hylaeamys megacephalus (HME) [20], were used to determine regions of chromosomal homology, and the mitochondrial genes COI and Cytb were used for molecular analysis. These results were compared with those from an undescribed species mentioned by da Silva et al. [28]. We present biogeographic inferences and discuss the chromosomal evolution of these taxa.

Material and methods
Animals collected during this study were handled following procedures recommended by the American Society of Mammalogists. JCP has a permanent field permit, number 13248 from "Instituto Chico Mendes de Conservação da Biodiversidade". The Cytogenetics Laboratory from UFPa has permit number 19/2003 from the Ministry of Environment for sample transport and permit 52/2003 for using the samples for research. The Ethics Committee (Comitê de Ética Animal da Universidade Federal do Pará) approved this research (Permit 68/2015). The rodents were maintained in the lab with food and water, free from stress, until their euthanasia using intraperitoneal injection of barbiturate (Pentobarbital, 120 mg/kg) after local anesthetic (lidocaine used topically).
We studied the karyotypes of nine specimens (five males and four females including one fetus) of Neacomys sp. A, collected from the right bank of the Tapajós River, in the Itaituba (Fig 1, localities 5, 6 and 8) and Jacareacanga municipalities (Fig 1, localities 7 and 9), Pará state, Brazil; seven specimens of Neacomys sp. B, three (one female and two males-one fetus) from the Itaituba municipality (Fig 1, locality 4) and four (three males and one female) from the Juruti municipality (Fig 1, locality 3), both from the left bank of the Tapajós River, Pará state, Brazil (S1 Table).
Samples were collected using Pitfall traps [33] and deposited at the mammal collection of Museu de Zoologia da Universidade Federal do Pará (UFPA), in Belém, Pará, Brazil. Chromosomal preparations were obtained from bone marrow [34] and by fibroblast cell culture made from two fetuses, established in the Centro de Estudos Avançados da Biodiversidade, Laboratório de Citogenética, Instituto de Ciências Biológicas, Universidade Federal do Pará, Belém, Pará, Brazil. G-banding and C-banding were performed following Sumner et al. [35] and Sumner [36], respectively. Fluorescent in situ Hybridization (FISH) studies were made using telomeric probes (All Telomere, ONCOR) and chromosome painting with whole chromosome probes of HME [20]. Twenty-one of the 24 probes of HME correspond to one HME chromosome pair, while three probes correspond to two pairs ( [9,10]; [13,22]; [16,17]). Digital images were obtained by Nis-Elements software and Nikon H550S microscopy. Chromosome classification followed Levan et al. [37].
The map was made using QUANTUM-GIS (QGIS) program version 2.10.1. Database were obtained from DIVA and IBGE.
We used sequences of 625 base pairs from 26 samples for Cytochrome C Oxidase-subunit I (COI), with 20 new sequences (Neacomys sp. A and Neacomys sp. B) and six sequences obtained from GenBank (S1 Table). For Cytochrome b (Cytb) we used sequences of 801 base pairs from 32 samples, with 18 new (Neacomys sp. A and Neacomys sp. B), and the others were kindly supplied by J.L. Patton (Museum of Vertebrate Zoology, Berkeley) or retrieved from GenBank (S1 Table). We included data from da Silva et al. [28] on an undescribed species (here mentioned as "Neacomys sp.") found in Marabá and Marajó Island (Fig 1, places 1 and 2).
The genetic distance between taxa was estimated with Molecular Evolutionary Genetics Analysis-MEGA v. 6.0 software [41], recovering K2P model.
The phylogenetic reconstructions were made using both the maximum likelihood (ML) method, run in RaxML v. 8 [42] with 1000 bootstrap replicates and Bayesian inference (BI) as implemented in MrBayes v. 3.2.1 [43]. In MrBayes, the analysis of substitution model parameters was unlinked across partitions. Two independent runs were initiated simultaneously with four independent Markov-Chain Monte Carlo (MCMC) chains (one cold and three heated). The MCMC algorithm was based on 700,000 cycles (generations), sampled every 5,000 cycles, with 20% of the samples being discarded as burn-in. Convergence was assessed by comparing the two runs. The MCMC output was visualized and diagnosed in Tracer v. 1.6 [44]. The run was considered satisfactory when, for all traces, the Effective Sample Size (ESS) values were over 200. Hylaeamys megacephalus, Oecomys rutilus, O. concolor, Deltamys and Thalpomys were used as outgroup. All these species belong to the Sigmodontinae subfamily and are phylogenetically close to Neacomys, according to Weksler [30].
Divergence time estimates were performed using BEAST 1.8.3 [45]. For calibration, we use three calibration points (4.4 Ma corresponding to separation time estimate between Oecomys and Hylaeamys [46]; 4.5 Ma corresponding to separation between Neacomys and Thalpomys; and 5 Ma corresponding to separation between Deltamys and the Neacomys/Thalpomys clade). Uncorrelated relaxed clock was assigned to the length rates among branches and Yule prior was used for the tree. Four independent runs were made of 20 5 generations, showing parameters and trees every 2,500 generations. The convergence of races was evaluated in Tracer v. 1.6 [44], assuming ESS values above 200 as satisfactory. Tree's and log file's results were summarized in TreeAnnotator v. 1.8.3 and LogCombiner v. 1.8.3 [47], respectively; we discard 20% as burn-in. The tree was displayed and edited in Figtree v. 1.4.2 (http://tree.bio.ed.ac.uk/ software/figtree/).
The other twelve autosomal probes show multiple signals in NSP-A, with ten (HME 1, 6, 7, 8, [9,10], 11, 14, [16,17], 19 and 23) hybridizing to two chromosomes each, where HME 1 and 8 hybridize to two whole distinct chromosomes each while the others hybridize to a chromosome and a portion of another chromosome. HME [13,22] show signals in three chromosomes and HME 5 in four different chromosomes.

Molecular phylogeny
The genus Neacomys was shown to be monophyletic in both analysis of maximum likelihood and Bayesian inference (Figs 5 and 6, respectively) for the concatenated mitochondrial genes (COI and Cytb; Table 3), supported by maximum values of bootstrap posterior probability. Six valid species were recovered for the clades : N. dubosti, N. guianae, N. minutus, N. musseri, N. paracou and N. spinosus. Besides, our data recovered two new clades (Neacomys sp. A and Neacomys sp. B), being monophyletic, with a high degree of divergence between them (Table 3), and other species in the Neacomys genus. NSP-A and NSP-B clades are described for the first time in this study.

Divergence time estimates and ancestral areas
Our divergence time estimates suggest that the diversification of species currently recognized for Neacomys genus occurred in the last 1.88 Ma. The last split was between Neacomys sp. and Neacomys sp. A about 0.45 Ma (Fig 7). Our data were unable to recover some Neacomys ancestor nodes with high statistical support values. However, they recovered with maximum support the ancestor of N. guianae, N. musseri and N. minutus that occurred in current areas of Inambari and Guiana. Thus, the ancestor of Neacomys sp. was endemic in the Marajó Island and Xingu area.

Karyotypic and phylogenetic analyses of Neacomys
Neacomys sp. A (2n = 58/FN = 68; Fig 1, localities 5-9) presents a similar karyotype to those described by Da Silva et al. [28] for an undescribed species, identified as Neacomys sp., for which two karyotypes are described: 2n = 58/FN = 64 for Marabá, in the southeastern portion of the state of Pará (Fig 1, locality 1) and 2n = 58/FN = 70 for specimens from Marajó Island (Fig 1, locality 2). Comparative analysis by G-and C-banding demonstrate that the differences  in FN among the three karyotypes are due to differences in heterochromatic blocks, in which CH forms the short arms of some bi-armed chromosomes. Neacomys sp. B (2n = 54/FN = 66) presents a new karyotype for the genus when compared with species with 2n between 56 and 64 (N. dubosti, N. guianae, N. paracou, Neacomys sp. and N. spinosus, Table 1). In NSP-B the bigger autosomes pairs are metacentric and submetacentric, while in other species they are medium-size acrocentrics, indicating multiple fusion/fission and/or translocation events.
According to Bradley and Baker [49] and Baker and Bradley [50], who made a meta-analysis of the genetic divergence in the Cytb gene for many groups of rodents, values of genetic divergence below 2% were present in different populations of the same species; over 5% are associated with potentially unrecognized species, and over 10% belongs to different species.
However, Neacomys sp. shows a mean divergence of 6.2% for Neacomys sp. A and 9.1% for Neacomys sp. B, while Neacomys sp. A and Neacomys sp. B present a medium nucleotide divergence of 7.4% from each other in concatenated mitochondrial genes (COI and Cyt b; Table 3). These three taxa present >10% of divergence from other Neacomys species in both analyses.
Thus, based on the genetic species concept [49,50] and the karyotypic and molecular data of this study, we conclude that Neacomys sp. A and Neacomys sp. B are two undescribed species within the genus and distinct from the undescribed species (Neacomys sp.) proposed by Da Silva et al. [28]. Moreover, these three undescribed species may represent cryptic species, which reinforces a taxonomic analysis to define their taxonomic status.

Chromosomal rearrangements and signatures
The comparison of Neacomys sp. A and Neacomys sp. B karyotypes show that these species had a karyotypic evolutionary history that involved complex rearrangements with some chromosomal signatures that differ them from other Sigmodontinae (see below), as also many autapomorphic characteristics for each species which confirm that this genus is very diverse even in karyotypes with not very distant 2n, and they differ from one another by 11 fusion/fission events and one translocation in 16 pairs of NSP-A and 14 pairs of NSP-B, plus four pericentric inversions in four autosomal pairs, and four CH amplification events in three autosomal pairs and the X chromosome. Only eight chromosomal pairs show conserved synteny with no detectable change ( Table 4, Fig 8). Fusion/Fission 3 (*HME 3) 6 (*HME 1) 2p (HME 3*1) 2q (HME 3*1) Fusion/Fission 22 (*HME 5) 9 (*HME [9,10]/7/[9,10]) 2 (*HME 2)

Fig 8. Comparative analysis by G-banding and ZOO-FISH with HME whole chromosome probes [20], between Neacomys sp. A and Neacomys
The absence of interstitial telomeric sequences (ITS; Fig 4B and 4D) may indicate that the rearrangements are old and that such sequences may have degenerated to the point of being undetectable by FISH [51,52]. Alternatively, the rearrangements may have occurred without involving telomeric sequences. Similar results are described for five karyotypes of three Neacomys species [28].
When we compare those authors' results with the NSP-A and NSP-B (Oryzomyini) karyotypes (Table 2), and extrapolating them using G-banding for the five karyotypes of three Neacomys species (Neacomys sp., N. dubosti and N. paracou) [28], we observe that the associations HME 20/ [13,22], 6/21 and 7/ [9,10], which are ancestral traits for Sigmodontinae according to Pereira et al. [23], are present also in Neacomys. However, these segments are rearranged in the genus and so they are synapomorphies: the first was due to a fusion with HME 4, generating HME 20/ [13,22]/4; the second was due to a fission in the segment that corresponds to HME 6, generating HME 6a/21 and 6b; the third was due to a fission in the HME 7 segment, followed by a paracentric inversion, generating HME [9,10]/7b/ [9,10] (Fig 9).
In our comparative analysis (Table 2), we observed another trait that could belong to the hypothetical ancestral karyotype of the Sigmodontinae subfamily: HME 15 non-associated, being a symplesiomorphy in TNI 19, CLA 12, NSP-A 15 and NSP-B 7. We consider that the acrocentric HME 24 is the ancestral form, being a symplesiomorphy in TNI 14, CLA 14, NSP-A 11 and NSP-B 11, and that the metacentric form (HME) and associated (AMO 6, ASP 3, NLA 9) are derivative. We also propose other ancestral traits for the Oryzomyini tribe: HME 8 disassociated in two fragments, HME 18 non-associated and HME 25 nonassociated.

Biogeography in Neacomys
The geographical barrier of the Amazonian rivers [6] explains the lack of gene flow between interpluvial regions in Amazon, and confines some species to a single endemic area [5], as described for several groups of terrestrial vertebrates, including primates [53], birds [54] and rodents [55]. In the Neacomys genus, some species occur in more than one endemic area [2,[24][25][26][27][28], which is in disagreement with the pattern observed for the three undescribed species of Neacomys, who have isolated distributions: Neacomys sp. within the Marajó island and Xingu endemic area [28], Neacomys sp. A within the Tapajós endemic area, and Neacomys sp. B within the Rondônia endemic area.
In contrast to the low node support observed for Neacomys sp. B, our divergence time estimates (Fig 7) recovered with high statistical support indicates that the geographical area of the ancestor of Neacomys sp. + Neacomys sp. A was the current Xingu and Tapajós areas of endemism and Marajó Island, while the Neacomys sp. ancestor area was the current endemic area of Xingu and Marajó Island.
Our divergence time estimates (Fig 7) suggest that the diversification of Neacomys sp. B and Neacomys sp. A + Neacomys sp. occurred about 0.74 Ma, and the last split was between Neacomys sp. A and Neacomys sp. about 0. 45 Ma. Based on speciation events in genus Psophia (Aves) and not on geological analyses, Ribas et al. [56] proposed that the Tapajós river drainage system was developed approximately 1.3-0.8 Ma, whereas Tocantins and Xingu rivers drainage systems were established about 0.8-0.3 Ma, acting as isolating barriers and creating the Tapajós, Xingu and Belém endemic areas.
Those divergence time estimates and diversification are within the range and in agreement with the gradient of chromosomal and molecular differentiation (see Karyotypic and phylogenetic analyses of Neacomys and Chromosomal Rearrangements and Signatures), which shows that Neacomys sp. [28], Neacomys sp. A and Neacomys sp. B form a monophyletic group, while the first two are sister species (Figs 5 and 6) and share more chromosomal similarities with each other than with Neacomys sp. B, that presents derivative chromosomal forms. Therefore, our data supports the hypothesis that the common ancestor from these taxa was distributed through the eastern Amazon and the Tapajós and Xingu rivers formation and also Marajó Island separation of the continent act as isolating barriers to gene flow and determine the pattern of diversification of these three undescribed species. Thus, our data provide strong support for the Riverine Barrier Hypothesis [7][8][9].
We emphasize that NSP-A and NSP-B were collected in areas not yet related to any other previously described species or distribution areas corresponding to them, thus enlarging the geographic distribution of the Neacomys genus [2,[24][25][26][27][28], for the southwestern region of the Pará state (Brazil). The number of species within the genus and their geographical boundaries remain uncertain [2].

Conclusions
The comparative chromosomal and molecular analyses in this study demonstrate that the Xingu and Tapajós Rivers act as geographic barriers for these three undescribed Neacomys species, delimiting Neacomys sp. distribution within the Marajó Island and Xingu endemic areas, NSP-A within the Tapajós endemic area and NSP-B within the Rondônia endemic area. In addition, we establish four synapomorphies for Neacomys (associations HME 20/ [13,22]/4, 6a/21, [9,10]/7b/[9,10] and 12/ [16,17]) and ancestral traits for the Oryzomyini tribe (HME 8a and 8b, 18 and 25) and Sigmodontinae subfamily (HME 15 and 24). It is important to continue using HME probes as taxonomic markers in other Sigmodontinae, for the definition of each tribe's chromosomal signatures and for the elucidation of taxonomic and phylogenetic relationships.