Oecomys catherinae (Sigmodontinae, Cricetidae): Evidence for chromosomal speciation?

Among the Oryzomyini (Sigmodontinae), Oecomys is the most speciose, with 17 species. This genus presents high karyotypic diversity (2n = 54 to 2n = 86) and many taxonomic issues at the species level because of the presence of cryptic species and the overlap of morphological characters. For these reasons the real number of species of Oecomys may be underestimated. With the aim of verifying if the taxon Oecomys catherinae is composed of more than one species, we made comparative studies on two populations from two regions of Brazil, one from the Amazon and another from the Atlantic Forest using both classical cytogenetics (G- and C-banding) and comparative genomic mapping with whole chromosome probes of Hylaeamys megacephalus (HME), molecular data (cytochrome b mitochondrial DNA) and morphology. Our results confirm that Oecomys catherinae occurs in the southeast Amazon, and reveal a new karyotype for the species (2n = 62, FNa = 62). The comparative genomic analysis with HME probes identified chromosomal homeologies between both populations and rearrangements that are responsible for the different karyotypes. We compared our results in Sigmodontinae genera with other studies that also used HME probes. These chromosomal differences together with the absence of consistent differentiation between the two populations on morphological and molecular analyses suggest that these populations may represent cryptic species.

Introduction similarity allows precise hybridization of probes onto chromosomes, and the production of reliable homology maps of the species.
The present study aimed to investigate if O. catherinae represents more than one species. For this, both classical cytogenetics (G-and C-banding) and comparative genomic mapping using whole chromosome probes of H. megacephalus [25] were used in addition to molecular analysis (mitochondrial cytochrome b DNA) and morphological studies.
Ethics Committee (Comitê de Ética Animal da Universidade Federal do Pará) approved this research (Permit 68/2015). The rodents were maintained in the laboratory for 48 hours at most, free from stress in suitable rat cages, since OCA are a little smaller than rats. Oecomys is granivorous and frugivorous and they were feed with seeds and regional fruits bought in the city markets. Animals received food and water ad libitum. The room temperature was regulated to 25ºC and the light cycle was the usual day and night period of time since this is a region close to the Equator line and there is no variation on this cycle all around the year. The animals were euthanized by intraperitoneal injection of barbiturate (Pentobarbital, 120 mg/ kg) after local anesthetic (lidocaine used topically).

Sample
Twenty-eight Oecomys catherinae individuals were studied, from which eight were karyotyped, 23 were included in the molecular phylogenetic analysis and 16 were used in the morphological analysis. All showed the characteristics described originally for the species by Thomas [30] (S1 File). Fig 1 shows the locations where the specimens were collected. More details are provided in S1 Table and S1 File.

Karyotype analysis
Chromosomal preparations were obtained from bone marrow of five specimens (three males and two females) collected in the municipality of Parauapebas, Pará (05˚21 0 54@S; 49˚07 0 24@W) and three specimens collected in the state of Rio de Janeiro, one male from the Estação Ecológica do Paraíso, Cachoeiras de Macacu municipality (22˚32'35''S; 42˚48'19''W), and a male and a female from the Fazenda Samburá, Cambuci municipality (21 o 29'35"S; 41 o 52'20"W) (Fig 1, S1 Table and Table 1). Chromosomal preparations of specimens from Pará were made according to Ford & Hamerton [31] and those from Rio de Janeiro were obtained from bone marrow short time culture. Cell suspensions were cultured for two hours at 37ºC in culture medium with RPMI 1640 supplemented by 20% bovine fetal serum, ethidium bromide (0.01 ml/ml in the culture medium) and colchicine 10 -6 M. After receiving a hypotonic shock for 40 min in KCl, cells in suspension were fixed in Carnoy (three parts of methanol and one of acetic acid). Staining was carried out with Giemsa and G-and C-banding followed Verma and Babu [32] and Sumner [33], respectively. FISH with telomeric probes (All Telomere Probes, Oncor) followed the manufacturer's protocol.

Molecular analysis
We used 25 O. catherinae tissue samples from the Atlantic Forest and Amazonia, of which 23 were produced in the present study. DNA was extracted using phenol-chloroform and proteinase K-RNAse protocol [35]. Partial cytochrome b (Cytb) sequences were amplified with polymerase chain reaction (PCR) with primers MVZ05 5'-CGAAGCTTGATATGAAAAACC ATCGTTG-3' and MVZ16 5'AAATAGGAARTATCAYTCTGGTTTRAT-3' [36]. The amplification protocol consisted of initial denaturation at 94˚C for 3 minutes, followed by 35 cycles of 30 seconds of denaturation at 94˚C, 1 minute of annealing at 45˚C and 2 minutes of extension at 72˚C, with a final extension at 72˚C by 7 minutes.
Sequences were edited using BioEdit 7.0.5.2 [37] and aligned with ClustalX 2.0.9 [38], following the proposed parameters of Schneider [39], with posterior manual rectification with BioEdit. Using jModeltest version 2.1.4 [40] we found GTR+I+G with a substitution rate equal to 6, gamma distribution parameter equal to 0.6270, and invariable sites in proportion equal to 0.46 as the best model for our sequences. Bayesian analysis (BI) was conducted on MrBayes 3.2.0 [41] using the evolutionary model described above, two runs, four chains, 10 million generations and sample frequency equal to 1000. The Maximum Likelihood (ML) was estimated with Garli 2.0 [42], using the above referred evolutionary model.
We also used sequences from 11 other species of Oecomys, H. megacephalus and Thomasomys andersoni, all available in GenBank (S1 Table). The latter two were used as outgroups in the phylogenetic analyses (following Rocha et al. [43]). Nonparametric bootstrapping based on 1000 replicates was performed to calculate node support values for ML. The genetic distance was estimated by MEGA 5.2 using the Kimura 2-parameter substitution model (K2P).

Morphological analysis
We examined external and craniodental characters of 16 specimens of O. catherinae, nine of which were from the Atlantic Forest of southeastern Brazil and seven from southeastern Amazonia (Fig 1; S1 Table). Among these specimens, 15 were considered adults because they exhibited fully erupted dentition. We extracted 12 craniodental measurements with digital calipers from the adult specimens based on Voss [44]. We used the Student t-test to compare mean craniodental dimensions between the specimens from Atlantic Forest and Amazonia. Statistical analyses were performed with the software SPSS. 13.0. for Windows, at a 5% significance level. A detailed description of the morphological analysis is provided in S1 File.

Classical cytogenetics
The diploid number of O. catherinae from Pará (OCA-PA) was 2n = 62 and FNa = 62, with a chromosomal complement of 29 pairs of acrocentric and one small pair of metacentric chromosomes (pair 27); the X chromosome is a large submetacentric and the Y is a medium sized acrocentric (Fig 2A). The Constitutive Heterochromatin (CH) is present at the pericentromeric region in all autosomal pairs; the short arm of the X chromosome is heterochromatic as is almost all the Y chromosome ( Fig 2B). Oecomys catherinae from Rio de Janeiro (OCA-RJ) possesses 2n = 60 and FNa = 62, composed of two pairs of chromosomes with two arms (pairs 1 and 27) and 27 pairs of acrocentric chromosomes. The X chromosome is a large submetacentric and the Y is a large acrocentric (Fig 3A). The CH is present in the pericentromeric region of all autosomal pairs; the small arm of the X chromosome is all heterochromatic and the Y is almost all heterochromatic (Fig 3B). The results of chromosome specific probe hybridization of H. megacephalus (HME; 2n = 54 e NFa = 62, [25]) to the karyotypes of OCA-PA (2n = 62 and FNa = 62) and OCA-RJ (2n = 60 and FNa = 62) are shown in Table 2. CH regions did not show hybridization signals. In both karyotypes of OCA (PA and RJ) the 24 probes showed 38 segments of homology (Figs 2A and  3A).

Morphological analysis
We were not able to find any consistent differences in qualitative traits between the Atlantic Forest and the Amazonian populations of O. catherinae herein examined. In contrast, Amazonian specimens exhibited a longer incisive foramen and a wider braincase, whereas specimens from the Atlantic Forest exhibited a wider rostrum and a deeper upper incisive, with no overlapping values in this latter trait (S3 Table). A detailed description of the morphological analysis is given as supplementary material (S1 File).

Geographic distribution of O. catherinae
The populations of O. catherinae here studied are from the southeast Atlantic Forest, previously known as the area of occurrence of the species, and from the southeast Amazon (municipalities of Marabá and Parauapebas, state of Pará and Vila Rica, state of Mato Grosso), in a different area according to the previously known locality of the species (Fig 1). Our registers of O. catherinae in Parauapebas and Vila Rica confirm the occurrence of this species in the southeast Amazon region, as proposed by Lambert et al. [13,14], who reported the possible occurrence of this species in the Pinkaití Research Station, in the southeast portion of the state of Pará, located 300 km from Parauapebas. In turn, Asfora et al. [4] considered a sample from Rondônia referred as Oecomys cf. concolor by Andrades-Miranda et al. [19] as O. catherinae, because of the karyotype of 2n = 60 and FNa = 62. However, this statement needs to be confirmed by morphological and molecular analyses.

Chromosomal variation in O. catherinae
The karyotype of OCA-RJ (2n = 60, FNa = 62, Fig 3A) is similar to the one already described in the literature (for a revision see Asfora et al. [4]) and different from the karyotype with 2n = 60/FNa = 64 caused by a pericentric inversion of one small pair of acrocentric chromosomes in the karyotype with FNa = 62 to a metacentric in the karyotype with FNa = 64. OCA-PA presents a new karyotype for the species, with 2n = 62/FNa = 62 (Fig 2A). Our comparative genomic mapping with HME probes allows the identification of all chromosomal homeologies, and rearrangements that differentiate the karyotypes. The results show that the karyotypes of OCA-RJ and OCA-PA differ by two chromosomal rearrangements, which can explain the difference in the 2n of 60 to 62 chromosomes, where HME6, which corresponds to two pairs (9 and 16; Figs 3A and 4C) in the karyotype of OCA-RJ and to three pairs (9, 25 and 26) in the karyotype of OCA-PA (Figs 2A and 4C). Possibly this difference arose from a tandem fusion or fission, where chromosomes 25 and 26 of OCA-PA correspond to chromosome 16 of OCA-RJ (Fig 4). In addition a translocation explains the difference in the number of metacentric pairs, in which HME2, homologous to pair 4 in the karyotype of OCA-PA (Figs  2A and 4B), corresponds to pair 4 of OCA-RJ plus the short arm of pair 1 (Figs 3A and 4B), that forms the HME 2/3 association (Fig 4A). Also, both karyotypes differ from each other by the size of the X and Y chromosomes, because of the difference in the amount of CH (Figs 2B and 3B). Such rearrangements may lead to problems during gametogenesis, producing gametes with low viability, and thus reduced fertility [45]. Similar results with similar conclusions were found by other authors on rodent chromosomal studies [9,46]. The reduction in fertility between populations may result in reproductive isolation, and this may be reinforced by the probably allopatric distribution of these populations. Alternatively populations were first geographically separated and, following this separation, chromosomal divergence occurred/ emerged and got fixed.
In studies of other Sigmodontinae genera using HME probes, it is possible to verify the organization of the syntenic group that corresponds to chromosomes 2 and 6 of H. megacephalus [25]. In Cerradomys langguthi [25], Akodon montensis, Thaptomys nigrita [29], Necromys lasiurus and Akodon sp. [24], the HME2 chromosome is divided into two syntenic blocks of similar size that correspond to chromosomes 4 and 5 of Mus musculus [24], the probable ancestral form. In OCA-PA these two blocks are united as in H. megacephalus, but in OCA-RJ they are fragmented in two blocks of unequal size, the smaller fragment resulting from a translocation (described previously). Thus, our results indicate that in Oecomys the original chromosomal form might be two united blocks (OCA-PA) and a translocation in one portion of this block to the homologous chromosome of HME3 (OCA-RJ) would be the derived form, probably an autapomorphic trait. In turn, the HME6 chromosome in those species mentioned above is present as a unique block, associated with HME21 (with the exception of H. megacephalus), corresponding to an ancestral character, the association HME6/21 and chromosome 2 of M. musculus [24]. Therefore, the syntenic groups corresponding to HME6 in Oecomys, both in OCA-PA as in OCA-RJ, are derivative. A more robust analysis of the genus Oecomys will be necessary to define with certainty which karyotype in OCA is the derived form.
branch of Oecomys catherinae the alphanumeric identity and the geographic locality is provided. Bold branches refer to the samples for which karyotypic analysis was performed. https://doi.org/10.1371/journal.pone.0181434.g005 In comparing the results of the morphological, molecular and cytogenetic analyses, we were able to observe significant chromosomal variation between the samples from the Amazon and the Atlantic Forest. There seems to have been absence of hybridization and introgression between both populations, meaning an absence of gene flow. While most morphological characteristics do not show consistent differences in the two populations, four craniodental measurements differ significantly between them (S3 Table). The molecular analysis (Fig 5) shows that the medium genetic divergence between populations is only 1.7% (SD = 0.6%), a value lower that those observed between species within the same genus (8.5% to 12.6%) [15]. These chromosomal differences described here, in association with the absence of consistent differentiation between the two populations in our molecular analysis, suggest that these populations could be cryptic species, in which craniodental differences could indicate the beginnings of morphological differentiation.

Conclusion
The present paper confirms the presence of O. catherinae in the Amazonian biome, enlarging its geographic distribution in South America. These populations also contain a new karyotype for the species, with 2n = 62/FNa = 62. Comparative genomic mapping shows that this karyotype differs from all others already published, by fission-fusion and translocation. However, these disjunct populations from the Amazon and Atlantic Forest may represent two cryptic species, with chromosomal data and geographic distance indicating the absence of gene flow between them, although morphological and molecular data do not show significant differences. Our data indicate that further studies and more fieldwork for increasing samples are welcome to improve knowledge about the diversity of this small mammal in the Amazon and Atlantic Forest.
Supporting information S1 File. Methods and results of the morphological analysis of the sample. (DOCX) S1 Table. List of specimens included in the present study. For each specimen, voucher and/ or field number, GenBank accession number, locality, and type of analysis are provided (Cyt = Cytogenetics; Mol = Molecular; Morph = Morphology). In bold, sequences produced in the present study. States are ES (Espírito Santo), MG (Minas Gerais), MT (Mato Grosso), PA (Pará), RJ (Rio de Janeiro), and SP (São Paulo). Localities are plotted on map (Fig 1). (DOCX) S2 Table. Estimates of evolutionary divergence over sequence pairs between groups. The numbers of base substitutions per site from averaging over all sequence pairs between groups are shown. Standard error estimate(s) are shown above the diagonal. Analyses were conducted using the Kimura 2-parameter model [1]. The rate variation among sites was modeled with a gamma distribution (shape parameter = 0.627). The analysis involved 38 nucleotide sequences. Codon positions included were 1st+2nd+3rd+Noncoding. All positions containing gaps and missing data were eliminated. There were a total of 352 positions in the final dataset. Evolutionary analyses were conducted in MEGA6 [2].