Introduction

Two of the most important naturalists from the turn of the 20^th Century were Edward William Nelson and Edward Alphonso Goldman. They contributed greatly to our knowledge, understanding, and documentation of the biota in the United States and México (^{López-Medellin and Medellin 2016}, https://sova.si.edu/record/SIA.FARU7364). The scientific material collected by both naturalists continues to be used as a rich resource in the systematic revision of many groups of birds and mammals (^{López-Medellin and Medellin 2016}). Nelson and Goldman’s biological surveys encompassed all of the states in México and lasted 14 years (1892 to 1906). In 1896, Nelson and Goldman conducted field work in Coahuila, México where they collected three individuals, later recognized as Peromyscus hooperi. These specimens were deposited and remain housed at the Smithsonian Institution’s National Museum of Natural History in Washington DC.

Peromyscus hooperi is a monotypic species, endemic to México and only known from portions of the states of Coahuila, Zacatecas, and San Luis Potosí (^{Álvarez-Castañeda 2002}). This species is sympatric with P. eremicus, P. melanophrys, and P. pectoralis in the states of Coahuila and Zacatecas (^{Schmidly et al. 1985}). Its preferred habitat is the grassland transition zone, a mixture of desert scrub and grassland vegetation (^{Schmidly et al. 1985}; ^{Lee and Schmidly 1977}). Its present fragmented and restricted distribution is considered a relict of a much larger historical distribution (^{Schmidly et al. 1985}).

Peromyscus hooperi is poorly represented in mammal collections and little is known about its current status in their restricted distribution; however, it is not protected by the Mexican government (Norma Oficial Mexicana - 059 - 2020, ^{Secretaría de Medio Ambiente y Recursos Naturales 201}0) and is classified as Least Concern by the International Union for Conservation of Nature - IUCN - (accessed on August 2022, ^{Álvarez-Castañeda 2016}). The species resembles P. eremicus, P. merriami, and P. pectoralis in cranial and external characters but differs in the karyotype (^{Lee and Schmidly 1977}; ^{Schmidly et al. 1985}). ^{Fuller et al. (1984}) and ^{Schmidly et al. (1985)} found that the karyotype of P. hooperi is very similar to P. crinitus, P. simulus, Osgoodomys banderanus and northern populations of P. boylii. However, P. hooperi has been described as a medium size mouse for the genus, with a long and bicolored tail (light grayish brown above and whitish below) with short hair. The upper parts, including face and top of head, are grayish with faint to moderate wash brown; lateral line is faint and near light buff; underparts are cream; and hind feet and lower legs are whitish. The skull contains large auditory bullae, and the first two upper and lower molars lack mesolophs. The glans penis is small but wide with a long protractile tip, and the baculum is long and slender with a cartilaginous tip ^{(Lee and Schmidly 1977}). The karyotype (2n = 48, FN = 52) comprises three pairs of biarmed autosomes and 20 pairs of acrocentric acrosomes (Lee and Schmidly 1977; ^{Schmidly et al. 1985}).

The taxonomic affinity of Hooper’s deer mouse has been problematic (^{Carleton 1989}). Based on a series of morphological characters (i. e., cranial characteristics, accessory lophs, and styles of the anterior molars, structure of the hyoid, and number and placement of the mammae) it was suggested to be closely related to the subgenus Haplomylomys (^{Lee and Schmidly 1977}). However, based on the anatomy of the phallus, it was more similar to species representing the subgenus Peromyscus (Lee and Schmidly 1977; ^{Schmidly et al. 1985}). Therefore, P. hooperi was characterized as an intermediate form between these two subgenera (^{Lee and Schmidly 1977}; ^{Fuller et al. 1984}; Schmidly et al. 1985). Peromyscus hooperi currently is recognized as the sole member of the Peromyscus hooperi species group (Schmidly et al. 1985; Carleton 1989), based on morphological characters, karyotypes, allozymes, and mtDNA - cytochrome b (Cytb; Carleton 1989; ^{Musser and Carleton 1993}, ²⁰⁰⁵; ^{Hogan et al. 1993}; ^{Dawson 2005}; ^{Bradley et al. 2007}).

^{Bradley et al. (2007}) used Cytb sequence data to conduct a phylogenetic analysis of the genus Peromyscus. They recovered strong nodal support for a sister group relationship between P. hooperi and P. crinitus with Maximum Parsimony (MP), however, using Maximum Likelihood (ML) and Bayesian Inference (BI) they did not resolve this relationship. In turn, this clade was sister to a clade including P. melanotis, P. polionotus, P. maniculatus, P. keeni, P. leucopus, P. gossypinus, P. eremicus, P. californicus, and Osgoodomys banderanus. ^{Platt et al. (2015}), included Cytb and three nuclear genes - Adh1-I2, Fgb-I7 and Rbp3, and concluded that the phylogenetic position of P. hooperi remains uncertain due to a lack of support values and the different placement between ML and BI analyses.

An additional problem for the systematic classification of the species within Peromyscus is the very definition of the genus. Several revisions and classifications have recognized subgenera - sensu lato - (^{Osgood 1909}; ^{Hooper and Musser 1964}; ^{Hooper 1968}) and genera - sensu stricto - (^{Carleton 1980}; ^{Carleton 1989}; ^{Musser and Carleton 2005}) within Peromyscus. However, the current resolution of this group does not fully adhere to either of those classifications. In addition, genetic and genomic studies have demonstrated the paraphyly of Peromyscus (^{Bradley et al. 2007}; ^{Miller and Engstrom 2008}; ^{Platt et al. 2015}; ^{Sullivan et al. 2017}; ^{Castañeda-Rico et al. 2022}). While clarifying the definition of Peromyscus is beyond the scope and objective of this manuscript, it is important to point out that whether we align to the sensu lato or sensu stricto classification of the genus, the phylogenetic placement of P. hooperi has not been well-resolved. However, hereafter, we recognized the genus Peromyscus as paraphyletic, including Habromys, Megadontomys, Neotomodon, Osgoodomys, and Podomys at the generic level (sensu stricto).

Uncertainty of the phylogenetic position of P. hooperi based on previous studies necessitates a revaluation using additional sequence data. To accomplish this, we used genome-wide data, including several thousand nuclear ultraconserved elements and whole mitochondrial genomes from a museum voucher specimen of P. hooperi collected by Nelson and Goldman combined with data from previous studies. These data provide new evidence about the phylogenetic placement of P. hooperi and its time of divergence from other peromyscines.

Materials and methods

Sample collection and laboratory methods. We used a museum specimen sample of Peromyscus hooperi - USNM 79619 - (ca. 2 mm² of dry skin) deposited at the Smithsonian Institution’s National Museum of Natural History; and collected by E. W. Nelson and E. A. Goldman on August 14, 1896 from Carneros, Coahuila, México. We followed strict protocols to avoid contamination during sampling, as described in ^{McDonough et al. (2018}) and ^{Castañeda-Rico et al. (2020)}. All pre-PCR steps were performed in a laboratory dedicated to processing ancient and historical DNA at the Center for Conservation Genomics, Smithsonian National Zoo and Conservation Biology Institute, Washington, DC. A silica column extraction protocol (^{McDonough et al. 2018}) was used to extract DNA. We quantified DNA with a Qubit 4 fluorometer (Thermo Fisher, Waltham, MA) using a 1x dsDNA HS assay and visualized DNA with a TapeStation 4200 System (Agilent Technologies, Santa Clara, CA) using High Sensitivity D1000 reagents. A dual-indexed library was prepared using the SRSLY PicoPlus NGS library prep kit (Claret Bioscience, LLC), according to the manufacturer’s protocol. We performed dual indexing PCR with TruSeq-style indices (^{Meyer and Kircher 2010}) using Kapa HiFi HotStart Uracil+ (Roche Sequencing), following the manufacturer’s protocol. This library was amplified with 12 cycles of PCR. We then pooled three PCRs from the same library before cleaning to increase DNA fragment representation. We cleaned the indexed library using 1.6X solid-phased reversible immobilization (SPRI) magnetic beads (^{Rohland and Reich 2012}), quantified concentration using a Qubit 4 fluorometer, and inspected size ranges and quality with a TapeStation 4200 System (conditions as mentioned above). Target-enrichment was performed to capture ultraconserved elements (UCE) and mitogenomes using the myBaits® UCE Tetrapods 5Kv1 kit (Daicel Arbor Biosciences) following the myBaits protocol v3, and the myBaits® Mito kit (Daicel Arbor Biosciences) for the house mouse Mus musculus panel, following the myBaits protocol v4. We amplified post-enrichment UCE and mitogenome libraries with 18 cycles of PCR using Kapa HiFi HotStart Ready Mix (Roche Sequencing), following the manufacturer’s protocol. A 1.6X SPRI magnetic bead clean-up was performed subsequently. We again quantified and visualized the enriched libraries using a Qubit 4 fluorometer and a TapeStation 4200 System, respectively (conditions as mentioned above). Finally, captured libraries were sequenced on a partial lane of a NovaSeq 6000 SP PE 2 x 150 base pairs (bp; Illumina, Inc., San Diego, CA, US) at the Oklahoma Medical Research Foundation, Oklahoma City (combined with samples from unrelated projects).

In addition to the data generated in this study, we also reanalyzed previously published data including the following: UCEs and full mitogenomes from ^{Castañeda-Rico et al. (2020}, ²⁰²²), as well as Cytb gene sequences from ^{Bradley et al. (2007}), ^{Platt et al. (2015}), and ^{Sullivan et al. (2017}; Table 1 and Appendix 1).

Ultraconserved elements. We analyzed the raw data following the PHYLUCE v1.6.7 pipeline with the default parameters (^{Faircloth 2016} https://github.com/faircloth-lab/phyluce). Illumiprocessor 2.10 (^{Faircloth 2013}) and Trim Galore 0.6.5 (https://github.com/FelixKrueger/TrimGalore) were used to trim adapters, barcode regions and low-quality bases. Reads were assembled into contigs using Trinity 2.8.5 (^{Grabherr et al. 2011}), and identified contigs matching UCE loci in the 5K UCE probe set (https://github.com/faircloth-lab/uce-probe-sets). A monolithic FASTA file was produced to extract sequences from each sample. We aligned FASTA sequences using MAFFT 7.4 (^{Katoh and Standley 2013}; ^{Nakamura et al. 2018}) and performed edge trimming. The resulting alignments were filtered to test them for various degrees of missing data (matrix completeness): 65 % matrix (35 % of the taxa missing for each UCE locus), 75 % matrix (25 % of taxa missing), 85 % matrix (15 % of taxa missing), 90 % matrix (10% of taxa missing), and 95 % matrix (5 % of taxa missing). Samples included in this dataset are shown in Table 1. We quantified informative sites with the PHYLUCE script phyluce_align_get_informative_sites.py. All of these analyses were performed on the Smithsonian Institution High Performance Computing Cluster (Smithsonian Institution, https://doi.org/10.25572/SIHPC).

We conducted a Maximum Likelihood (ML) analysis using RAxML 8.12 (^{Stamatakis 2014}) with a GTRGAMMA site rate substitution model and 20 ML searches for the phylogenetic tree for each of the aforementioned data matrices (i. e., 65 % to 95 % matrices). Nonparametric bootstrap replicates were generated using the -N autoMRE option which runs until convergence was reached. We reconciled the best fitting ML tree with the bootstrap replicate to obtain the final phylogenetic tree with support values using the -f b command.

We estimated the best evolutionary model of nucleotide substitution in jModelTest 2.1.1 (^{Guindon and Gascuel 2003}; ^{Darriba et al. 2012}) using the Akaike Information Criterion (AIC). The TVW+G model was selected as the best fitting model with the following parameters: base frequencies A = 0.2988, C = 0.2013, G = 0.2026, T = 0.2972; nst = 6; and gamma shape = 0.1070. A Bayesian Inference analysis (BI) using MrBayes 3.2.6 (^{Huelsenbeck and Ronquist 2001}; ^{Ronquist and Huelsenbeck 2003}) was performed on the 90 % matrix. The BI analyses comprised two independent runs with 50 million generations, sampling trees and parameters every 1,000 generations with four Markov-chains Monte Carlo (MCMC), three heated and one cold. We visualized output parameters using Tracer v1.7.1 (^{Rambaut et al. 2018}) to check for convergence between runs and we discarded the first 25 % of the trees as burn-in.

A species tree analysis under the multispecies coalescent (MSC) model with ASTRAL-III v.5.7.8 (^{Zhang et al. 2018}) was performed on the 90 % matrix. The local posterior probability - LPP - (^{Sayyari and Mirarab 2016}) was used as branching support. We used the uce2speciestree pipeline script (^{Campana 2019} https://github.com/campanam/uce2speciestree) to generate input files for ASTRAL. This script uses RAxML to infer individual gene trees under the GTRGAMMA substitution model, and 100 bootstrap replicates.

Mitogenomes. FASTQ files were analyzed using FastQC v0.11.5 (^{Andrews 2010}, www.bioinformatics.babraham.ac.uk/projects/fastqc). Adapter sequences and low-quality reads were removed using the default parameters (Phred:20, mean min-len:20) in Trim Galore 0.6.5 (https://github.com/FelixKrueger/TrimGalore). Exact duplicates were removed (-derep1,4) using Prinseq-lite v0.20.4 (^{Schmieder and Edwards 2011}). We mapped the resulting high-quality reads to the closest available reference genome (Peromyscus crinitus - GenBank accession number KY707308), using the Geneious algorithm in Geneious Prime® 2021.2.2 (https://www.geneious.com) with default parameters (Medium-Low sensitivity, Maximum mismatches = 20 %, Maximum gaps = 10 %). A consensus sequence was generated with Geneious Prime® 2021.2.2 (https://www.geneious.com), using 4X as the lowest coverage to call a base, and aligned them using MAFFT 7.45 plug-in (^{Katoh and Standley 2013}). We transferred annotations from the reference to rule out the presence of nuclear copies of mitochondrial genes (NUMTs), and translated all protein-coding genes to check for frame shifts or stop codons.

We aligned sequences with MAFFT 7.45 plug-in (^{Katoh and Standley 2013}) in Geneious Prime® 2021.2.2 (https://www.geneious.com). Samples included in this dataset are listed in Table 1. A BI analysis was conducted on a partitioned dataset using MrBayes 3.2.6 (^{Huelsenbeck and Ronquist 2001}; ^{Ronquist and Huelsenbeck 2003}). The best model and partition scheme were estimated using PartitionFinder 2.1.1 (^{Lanfear et al. 2016}). Our search was limited to the models available in MrBayes, with linked, corrected Akaike Information Criterion (AICc) and greedy parameters. The data block was defined by gene, tRNA, rRNA and D-loop selection. We conducted two independent runs with 50 million generations, sampling trees and parameters every 1,000 generations with four MCMC and parameters as mentioned above, to perform the BI analysis. Convergence between runs was checked using Tracer v1.7.1 (^{Rambaut et al. 2018}), and the first 25 % of the trees was discarded as burn-in.

We performed a ML analysis using the concatenated dataset in RAxML 8.12 (^{Stamatakis 2014}) with a GTRGAMMA site rate substitution model. Clade support was assessed by bootstrapping with the -N autoMRE option for a bootstrap convergence criterion. The -f b option was used to reconcile the best fitting ML tree with the bootstrap replicate to obtain the final phylogenetic tree (as mentioned above). DNA damage patterns were evaluated for the P. hooperi sample with mapDamage2.0 (^{Jónsson et al. 2013}) using -- rescale option.

Cytochrome b. We analyzed Cytb sequences extracted from the mitogenome that was generated in this study and from mitogenomes published by ^{Sullivan et al. (2017}) and ^{Castañeda-Rico et al. (2020}, ²⁰²²). We also used the Cytb sequences published by ^{Bradley et al. (2007}) and ^{Platt et al. (2015}) in order to compare the phylogenetic position of P. hooperi using genome-wide data as well as a single mitochondrial gene. The Cytb dataset allowed us to include more species within the genus Peromyscus and representatives of the genera Habromys, Megadontomys, Neotomodon, Osgoodomys, Podomys, Isthmomys, Onychomys, Reithrodonyomys, Neotoma, Ochrotomys, Baiomys, Ototylomys, Tylomys, Nyctomys, Oryzomys and Sigmodon. Samples included in this dataset are shown in Table 1 and Appendix 1.

The Cytb dataset was analyzed as follows: we performed alignment using MAFFT 7.45 plug-in (^{Katoh and Standley 2013}) in Geneious Prime® 2021.2.2 (https://www.geneious.com). We estimated the best evolutionary model of nucleotide substitution in jModelTest 2.1.1 (^{Guindon and Gascuel 2003}; ^{Darriba et al. 2012}) using the AIC method. The TPM3uf+I+G model was selected as the best fitting model with the following parameters: base frequencies A = 0.3896, C = 0.3336, G = 0.0500, T = 0.2267; nst = 6; proportion of invariable sites = 0.4080; and gamma shape = 0.6220. A BI analysis was run using MrBayes 3.2.6 (^{Huelsenbeck and Ronquist 2001}; ^{Ronquist and Huelsenbeck 2003}) as mentioned above for UCE and mitogenome datasets. We used Tracer v1.7.1 (^{Rambaut et al. 2018}) to check for convergence between runs, and the first 25 % of the trees was discarded as burn-in.

Divergence times estimation. Molecular dates of divergence were estimated in BEAST2 v2.6.6 (^{Bouckaert et al. 2019}) using the mitogenome dataset. First, we obtained the best model and partition scheme in PartitionFinder 2.1.1 (^{Lanfear et al. 2016}). The search was limited to the models available in BEAST, linked branch lengths, AICc model selection, and greedy schemes search. The data block was defined by codon position, tRNA, rRNA and D-loop selection, and the result was incorporated in the dating analysis. The BEAST analysis was performed under an uncorrelated lognormal relaxed molecular clock model. The calibrated Yule speciation processes model (^{Heled and Drummond 2012}) with a randomly generated starting tree were set up as priors. We used the same three calibration points with a lognormal distribution from ^{Castañeda-Rico et al. (2022)}. Calibrations were based on fossil records of 1) Reithrodontomys (mean = 1.8 million years ago [mya], stdev = 1.076, offset = 1.63), as used by ^{Steppan and Schenk (2017}); 2) Onychomys (mean = 4.9, stdev = 1.169, offset = 4.753), as used by Steppan and Schenk (2017); and 3) the most recent common ancestor of P. attwateri (mean = 2.7, stdev = 0.9, offset = 2.4 [^{Dalquest 1962}; ^{Karow et al.1996}; ^{Wright et al. 2020}]). Two independent runs of 50 million iterations were performed, each was sampled every 1,000 iterations. We checked convergence statistics for Effective Sample Sizes (ESS) using Tracer v1.7.1 (^{Rambaut et al. 2018}) and a 25 % of burn-in was performed on each run. We used LogCombiner v2.6.6 to combine trees and TreeAnnotator v2.6.2 to get a consensus tree with node height distribution (both packages available in BEAST). All phylogenetic and ultrametric dated trees from the UCE, mitogenome and Cytb datasets were visualized in FigTree 1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/). All analyses were performed on the Smithsonian Institution High Performance Computing Cluster (Smithsonian Institution https://doi.org/10.25572/SIHPC).

Results

Following the PHYLUCE v1.6.7 pipeline, we recovered 1,087 UCE loci (raw data are available in GenBank under BioProject PRJNA880321), and a complete mitogenome of 16,288 bp (GenBank accession number OP432689) from the P. hooperi sample. The average number of paired-end reads and fragment size after trimming were 13,075,112 reads and 67 bp long, respectively. The lowest-quality bases were detected at the end of the reads. We also recovered between 1,353 and 3,859 UCE loci from the reanalyzed samples. The average number of paired-end reads and fragment size after trimming for those samples ranged from 1,811,856 to 21,093,430 reads, and from 94 to 222 bp, respectively.

Ultraconserved Elements phylogenies. We recovered 9,840 contigs for P. hooperi after Trinity assemblies. The mean, minimum, and maximum length for contigs were 242, 201, and 3,784 bp, respectively. The incomplete matrix contained 4,406 UCE loci (n = 18, average = 3,136, min = 1,087, max = 3,859). A total of 1,087 UCE loci were obtained for P. hooperi with a mean, minimum, and maximum length of 235, 201, and 636 bp, respectively. The 65 % matrix contained 3,681 UCE loci (UL) with an average of 13.80 informative sites per locus (IS), the 75 % matrix (UL = 2,974, IS = 14.18, the 85 % matrix (UL = 1,514, IS = 14.29), the 90 % matrix (UL = 677, IS = 14.07), and the 95 % matrix (UL = 168, IS = 14.30).

The datasets representing various levels of matrix completeness yielded the same ML topology with high support values for all branches (Figure 1, phylogenetic trees obtained from the 65 %, 75 %, 85 %, and 95 % matrices are not shown). The BI tree topology, based on the 90 % matrix, showed the same topology with high posterior probability values for all branches (Figure 1). Both, ML and BI trees placed P. hooperi as sister to the clade containing Podomys floridanus, Neotomodon alstoni, P. mexicanus, P. megalops, P. melanophrys, P. perfulvus, P. aztecus, Habromys simulatus, P. attwateri, and P. pectoralis. The species tree supported, with high LPP values, the same phylogenic position of P. hooperi (Figure 1, based on the 90 % matrix). The only difference among the species tree and the concatenated ML and BI trees, was the relationship between P. mexicanus and P. megalops. These two species are sisters in the ML and BI trees but not in the species tree, where P. megalops is sister to the clade containing P. mexicanus, P. melanophrys, P. perfulvus, P. aztecus, Habromys simulatus, P. attwateri, and P. pectoralis.

Figure 1 (a) Bayesian Inference and Maximum Likelihood phylogenies based on a 90 % matrix UCE with 677 loci. Nodal support is provided with posterior probability/bootstrap values. (b) Species tree based on a 90 % matrix UCE with 677 loci. Nodal support is provided with local posterior probability values. The blue block highlights the phylogenetic position of Peromyscus hooperi. 

Mitogenome phylogenies. The final alignment included 21 taxa and was 16,272 bp in length. BI and ML analyses, with six partitions, provided slightly different topologies (Figure 2). However, both analyses supported (pp = 1, bootstrap = 76) the placement of P. hooperi as sister to the clade including Podomys floridanus, Neotomodon alstoni, P. mexicanus, P. megalops, P. melanophrys, P. perfulvus, P. boylii, P. aztecus, Habromys ixtlani, P. attwateri, and P. pectoralis. The phylogenetic position of P. crinitus changed across phylogenies (Figure 2), as did the position of the clade containing Podomys floridanus + Neotomodon alstoni. However, the BI tree yielded higher support values. The DNA damage analysis showed a weak signal of damage, typical of historical DNA (Appendix 2).

Cytochrome b phylogeny. The alignment included 64 taxa, 154 sequences, and was 1,143 bp in length. The BI analysis placed P. hooperi sister to the clade containing P. maniculatus, P. polionotus, P. keeni, P. melanotis, P. leucopus, and P. gossypinus (Appendix 3). However, the branch support value for this phylogenetic position was low (pp = 0.53). The two samples of P. hooperi, one sequenced in this study (USNM 79619) and the other by ^{Bradley et al. (2007}; TTU 104425, GenBank accession number DQ973103) clustered together with high support (pp = 1).

Divergence time estimation of Peromyscus hooperi. The mitochondrial divergence dating analysis, with six data partitions, estimated a Pliocene divergence time for P. hooperi around 3.98 mya (95 % HPD: 3.57 to 4.47 mya; Figure 3). The divergence of P. crinitus was dated ca. 4.31 mya (95 % HPD: 3.80 to 4.70 mya), the split of the clade including P. leucopus + (P. polionotus + P. maniculatus) at ca. 4.49 mya (95 % HPD: 4.03 - 5.02 mya), and the divergence of the genus Peromyscus was dated ca. 5.21 mya (95 % HPD: 4.79 - 5.71 mya).

Figure 2 Mitogenome phylogenies based on Bayesian Inference (a) and Maximum Likelihood (b). Nodal support is provided with posterior probability and bootstrap values, respectively. The blue block highlights the phylogenetic position of Peromyscus hooperi. 

Discussion

The biological expeditions undertaken by Nelson and Goldman in México were arguably among the most important ever achieved by two naturalists for a single country (^{López-Medellin and Medellin 2016}; ^{Guevara 2021}; https://sova.si.edu/record/SIA.FARU7364). To our knowledge, this is one of a few studies in which genome-wide data were obtained and analyzed from a specimen collected by these two naturalists (see ^{McDonough et al. 2022}). Our results not only provide new evidence about the phylogenetic position of P. hooperi but also joins a short list of mammal studies within the blooming field of Museomics (see ^{Card et al. 2021} for a review) that have successfully analyzed specimens collected before 1900 within a phylogeny (e. g.,^{Abreu-Jr et al. 2020}; ^{Sacks et al. 2021}; ^{Roycroft et al. 2021}, ²⁰²²; ^{Castañeda-Rico et al. 2022}; ^{McDonough et al. 2022}; ^{Tavares et al. 2022}).

Our nuclear DNA results strongly support P. hooperi as sister to a clade containing Podomys floridanus, Neotomodon alstoni, Habromys simulatus, P. mexicanus, P. megalops, P. melanophrys, P. perfulvus, P. aztecus, P. attwateri, and P. pectoralis (all Peromyscus species within the subgenus Peromyscus). In the mitogenome analyses, P. boylii (subgenus Peromyscus) and H. ixtlani joined the sister group of P. hooperi (Figure 1, 2). However, our results do not agree with those of ^{Bradley et al. (2007}), who found low support for P. hooperi as sister to P. crinitus (subgenus Peromyscus, Peromyscus crinitus species group), and both species sister to a clade including P. melanotis, P. polionotus, P. maniculatus, P. keeni, and P. leucopus (subgenus Peromyscus, Peromyscus leucopus and maniculatus species groups), P. gossypinus, P. eremicus, and P. californicus (subgenus Haplomylomys, Peromyscus californicus and eremicus species groups), and Osgoodomys banderanus. ^{Platt et al. (2015}) showed that P. hooperi could be related with the same species suggested by ^{Bradley et al. (2007)}, although P. polionotus and P. keeni were not included in their study. However, the phylogenetic position of P. hooperi remained uncertain due to lack of strong nodal support in both of these previous studies.

Figure 3 Dated whole mitochondrial genome phylogeny. Dates are provided in millions of years. The horizontal bars and numbers below the branches show the 95 % Highest Posterior Density. The blue block highlights the phylogenetic position of Peromyscus hooperi. 

Our phylogenomic analyses strongly support the placement of P. hooperi with the Peromyscus mexicanus, megalops, aztecus, melanophrys, and truei species groups (all within the subgenus Peromyscus). We did include three out of the five species groups studied by ^{Bradley et al. (2007}). We analyzed the only member of the Peromyscus crinitus species group in the nuclear and mitogenome dataset, and members of the Peromyscus maniculatus and leucopus species group only in the mitogenome dataset; but we did not find that P. hooperi is closely related to any of those groups as previously suggested. Despite the novel data generated here, denser taxon sampling is still required to better confirm and/or determine the closest relative of P. hooperi. For example, phylogenetic relationships between P. hooperi and members of the subgenus Haplomylomys still require further testing. However, despite this limitation, here we have provided strong nodal support for P. hooperi for the first time.

The Cytb analysis (Appendix 3) confirmed the identity of the P. hooperi specimen used in this study (USNM 79619), placing it in the same clade with the only other P. hooperi Cytb sequence available (^{Bradley et al. 2007}, TTU 104425 and GenBank accession number DQ973103). In addition, the phylogenetic position of the species in this analysis is similar to ^{Bradley et al. (2007)} and ^{Platt et al. (2015}). We found that P. hooperi is most closely related to the Peromyscus leucopus and maniculatus species groups but with a low support (pp = 0.53); therefore, its phylogenetic position is not resolved. In conclusion, we confirmed that the phylogenetic position of the Hooper’s deer mouse cannot be resolved using only Cytb sequences or a few genes, as Platt et al. (2015) documented. Our results demonstrate that genome-wide data allow a better resolution of the phylogenetic relationships of phylogenetically problematic species.

Our divergence times estimations indicated that the crown of Peromyscus was estimated ca. 5.21 mya (95 % HPD: 4.79 to 5.71 mya), and the diversification of the genus occurred ca. 4.49 mya (95 % HPD: 4.03 to 5.02 mya), both events during the Pliocene. We dated the split of P. hooperi during the Pliocene at ca. 3.98 mya (95 % HPD: 3.57 to 4.47 mya), following the split from P. crinitus at ca. 4.31 mya (95 % HPD: 3.80 to 4.70 mya). These dates coincide with previously dated phylogenies obtained from genome-wide data of peromyscines (e. g.,^{Castañeda-Rico et al. 2022}). They estimated the crown of the genus Peromyscus during the Pliocene at ca. 5.32 mya (95 % HPD: 4.85 to 5.98 mya), and the origin of P. crinitus at ca. 4.62 mya (95 % HPD: 4.05 to 5.28 mya), using mitogenomes. Our results also show that the Peromyscus hooperi, crinitus, maniculatus, and leucopus species groups were among the first to diverge within the genus Peromyscus (Figure 3), followed by the Peromyscus megalops, mexicanus, melanophrys, boylii, aztecus, and truei species groups, together with Neotomodon, Podomys, and Habromys. Based on our results and those of previous studies (e. g.,^{Hibbard 1968}; ^{Riddle et al. 2000}; ^{Dawson 2005}; ^{Castañeda-Rico et al. 2014}, ²⁰²²; ^{Platt et al. 2015}; ^{Sawyer et al. 2017}; ^{León-Tapia et al. 2021}), we suggest the Pliocene and Pleistocene as the time when speciation and diversification events took place within peromyscines, potentially associated with climatic cycles related to numerous vicariant and dispersal events.

Previous phylogenetic studies of the genus Peromyscus that analyzed single or few genes, provided older divergence times estimations (e. g.,^{Castañeda-Rico et al. 2014}; ^{Platt et al. 2015}; ^{Cornejo-Latorre et al. 2017}; ^{Bradley et al. 2019}). For example, ^{Platt et al. (2015)}, using Cytb, estimated the origin of Peromyscus and its diversification, during the Miocene, at approximately 8 mya and 5.71 mya, respectively; and the divergence of P. hooperi at ca. 5.2 mya, during the early Pliocene. However, estimates of the time to the most recent common ancestor (TMRCA) calculated from individual or few genes can be overestimated (^{Duchêne et al. 2011}).

The evolutionary uniqueness of P. hooperi is supported by our results and previous studies by ^{Fuller et al. (1984}) ^{and Schmidly et al. (1985)} who found that this species does not fit well with either of the subgenera Haplomylomys or Peromyscus. We hypothesize that P. hooperi will remain the sole member of the Peromyscus hooperi group as first proposed by ^{Schmidly et al. (1985)} and later supported by ^{Carleton (1989}) based on the morphological, karyotypic, and allozyme evidence.

The genetic and morphological uniqueness of P. hooperi, as well as its restricted distribution to grassland transition zones should make this a species of special concern for conservation. In addition, ^{Schmidly et al. (1985)} stated that P. hooperi is a relictual, monotypic species without close living relatives, and its survival is jeopardized/threatened by the fragile conditions of its habitat in central Coahuila as a result of overgrazing. During the last 21 years, habitat shifts from native grasslands to crops zones have increased with agricultural intensification, grain-fed cattle feedlots, and new land use policies in the Mexican states of Durango, Sinaloa, Chihuahua, Nuevo León, and particularly Coahuila where P. hooperi is mostly distributed (^{Galván-Miyoshi et al. 2015}; ^{Bonilla-Moheno and Aide 2020}). We recommend that future studies conduct population genetic analyses to determine the genetic diversity and structure of the different populations of P. hooperi. This species remains poorly known and potentially threatened by habitat loss, therefore new information is needed to determine an appropriate conservation strategy and category.