The family Aplodontiidae contains a single, monotypic extant genus, Aplodontia (mountain beaver), which was first described by Rafinesque in 1817. Phylogenetic studies have shown that it is the sister lineage to squirrels. Aplodontia rufa is endemic to the Pacific Northwest and ranges from central California to British Columbia, Canada. Currently, 7 described subspecies are recognized based on morphological taxonomic studies. In this study, mitochondrial and nuclear genes were sequenced to infer molecular phylogenies of A. rufa. One of the goals of this study was to use molecular data to test the current taxonomic hypothesis based on morphology. Another goal was to incorporate geographic information to elucidate distributions of major clades. Our results support the previously held subspecies designations based on morphological taxonomy, with 1 main exception: we determined that within A. rufa, the subspecies A. rufa rainieri and A. rufa rufa north of the Columbia River represent a single lineage and should revert to the name A. rufa olympica. Although we revised geographic boundaries for some groups (A. r. rufa, A. r. olympica, and A. r. pacifica), only the conservation status and management of A. r. olympica (previously 2 subspecies) in Canada may be affected. Our findings support the continued conservation efforts for the isolated and endangered lineages present in coastal California.
The mountain beaver, Aplodontia rufa, the sole extant species from a monotypic genus within the family Aplodontiidae, is endemic to the Pacific Northwest region of the United States and into Canada and has a limited distribution from southern British Columbia, Canada, into northern California (Hall 1981; Fig. 1). A. rufa is a semifossorial rodent weighing approximately 1 kg that tends to occupy moist environments or areas near water sources (Godin 1964; Hacker and Coblentz 1993). Aplodontia has limited renal function due to inefficient kidneys; this “primitive” structure is considered to be a physiological characteristic that limits its geographical distribution and provides evidence that this is a member of an ancient lineage (Pfeiffer et al. 1960). In fact, aplodontioids date back to the late Oligocene (Shotwell 1958; Rensberger 1975; Vianey-Liaud 1985; Hopkins 2008) and recent molecular studies have shown modern Aplodontiidae to be sister to the Sciuridae (squirrel) clade (Huchon et al. 1999; Adkins et al. 2001; Herron et al. 2004).
Rafinesque (1817) first described Anisonyx rufa from accounts and specimens of the Lewis and Clark Expedition. Later this species was placed in the genus Aplodontia (Richardson 1829a) and subspecies were described as they were discovered (Richardson 1829b; Merriam 1886, 1899). Taylor (1918) revised the genus Aplodontia and designated a single species with 9 subspecies based on cranial and pelage characters. Geographic variation within A. rufa was investigated along the Columbia River system by Finley (1941); at that time he acknowledged pronounced skull and pelage differences between populations on either side of the river but did not revise subspecific designations. Dalquest and Scheffer (1945) followed up with a systematic study of A. rufa that focused on subspecies from Washington and reduced the number of subspecies there from 4 to 2. A later revision by Hall (1981) affirmed this classification, and thus there are currently 7 recognized subspecific lineages based on geographical distributions and morphological characters, principally pelage coloration and skull measurements (Dalquest and Scheffer 1945; Hall 1981; Fig. 1).
Several investigators of Aplodontia have remarked on variability in pelage and cranial features among populations (Taylor 1918; Finley 1941; Dalquest and Scheffer 1945). Pelage coloration has commonly been used as a taxonomic character but the variability of pelage within subspecies, populations, or even litters can make it an unreliable character to infer taxonomic differences (Lindenmayer et al. 1995; Koyabu et al. 2008). Furthermore, some skull measurements used to differentiate subspecies of A. rufa were so variable that Dalquest and Scheffer (1945) deemed them to be useless. Therefore, the morphological characters currently used to delimit all subspecies of A. rufa may not be robust characters for inference of evolutionary lineages and hence taxonomy. We propose that taxonomic designations within this genus would benefit greatly from a molecular assessment.
Four of the 7 subspecies, A. r. californica, A. r. rainieri, A. r. rufa, and A. r. pacifica, have sizeable ranges; whereas the other 3 subspecies, A. r. humboldtiana, A. r. nigra, and A. r. phaea, have very restricted and isolated distributions (Fig. 1). One of these, A. r. nigra, is listed as endangered by the United States Fish and Wildlife Service (United States Fish and Wildlife Service 1991). A. r. phaea is considered a “species of special concern” by the state of California, whereas A. r. humboldtiana is found in greater densities and therefore does not have a conservation status designation (Williams 1986). The broader-ranging A. r. rainieri and A. r. rufa are classified as “vulnerable” and “critically imperiled, or imperiled,” respectively, in Canada (Gyug 2000); hence, they are managed following separate plans. Conversely, these 2 subspecies are managed as pest species in Washington and Oregon, where they rapidly colonize timber areas that have been recently harvested or reforested, or both (Hacker and Coblentz 1993; Arjo et al. 2007). The preferred diet of these subspecies of newly planted seedlings and even sapling-aged trees negatively impacts forest production (Cafferata 1992; Arjo 2010).
For conservation and management efforts to be successful, taxonomy and distribution must be clearly understood. In particular we were concerned about the taxonomic designations of subspecies in California (A. r. humboldtiana, A. r. nigra, and A. r. phaea) as well as the most northern subspecies in Canada (A. r. rufa and A. r. rainieri) because of their conservation status. Therefore, our primary goal was to use a molecular phylogenetic approach to test the current taxonomy of A. rufa (Hall 1981) as a hypothesis of evolutionary relationships. We inferred phylogenetic relationships using 2 mitochondrial DNA (mtDNA) genes and a nuclear DNA (nDNA) gene; cytochrome b (Cytb), control region (CR), and exon 10 of the growth hormone receptor (GHR), respectively. Secondarily, this study aimed to test the geographical distributions of subspecies when possible through thorough sampling efforts. Finally, divergence times were estimated for lineages within the genus Aplodontia to provide insight into the historical processes that have influenced evolutionary relationships and current geographical distribution of this taxon.
Materials and Methods
Sampling
The majority of samples were obtained from National Wildlife Research Center research projects as contributions from United States federal agencies, state agencies, or academic institutions; and opportunistically from population-control efforts (Table 1). Samples came in various forms: tissue, hair, and blood. We solicited samples from a variety of sources to achieve as much coverage as possible of geographical distribution of each designated subspecies. Trapping (in most cases) primarily occurred along slopes and draws in recently harvested timber units or in open or disturbed areas near old timber units where mountain beaver activity was detected (Ransome 2003; Arjo et al. 2007). Tissue samples consisted of a small notch removed from the ear and stored in a buffer (20% dimethylsulfoxide, 0.25 M ethylenediaminetetraacetic acid, saturated with NaCl, pH 8.0 solution—Seutin et al. 1991) or 100% ethanol, and blood was absorbed by filter paper. Blood samples were air dried and placed in coin envelopes; hair samples also were placed in coin envelopes after collection. Samples were sent to the National Wildlife Research Center or the Genetic Data Centre at the University of British Columbia for molecular analysis. Animals trapped during the course of this study were captured and handled under a protocol approved by the Institutional Animal Care and Use Committee of the National Wildlife Research Center and in compliance with guidelines of the American Society of Mammalogists (Sikes et al. 2011).
Extraction of DNA, amplification, and sequencing
Two methods were used for the isolation of DNA depending on the laboratory (National Wildlife Research Center or Genetic Data Centre). At the National Wildlife Research Center, whole genomic DNA (gDNA) was extracted from tissue using the DNeasy Blood & Tissue Kit (Qiagen, Valencia, California) following the manufacturer's protocol and from hair and blood samples with the QIAmp DNA Micro Kit (Qiagen) using a 30-μl and a 75-μl elution, respectively. At the Genetic Data Centre, gDNA was extracted from tissue following standard proteinase K–phenol–chloroform procedures (Sambrook et al. 1989).
Three sets of primers were used to amplify and sequence portions of 2 genes at the National Wildlife Research Center. The mtDNA Cytb gene was amplified with 2 sets of primers: MVZ 05/04 (Smith and Patton 1993—for the first 450 base pairs [bp] and including a transfer RNA [tRNA]) and self-designed A. rufa cytbII F/R (for the remaining portion of the Cytb gene; A. rufa cytbIIF 5′ AACATGAAACATTGGCATTC 3′ and A. rufa cytbIIR 5′ ATAGAGGCCAGTTGACCAAT 3′). Aligning sequences from these 2 sets of primers resulted in Cytb gene sequences approximately 1,050 bp in length. Also, a portion of the nuclear DNA exon 10 of the GHR gene was amplified and sequenced using the GHR EXON10F/750R primers (Adkins et al. 2001). Details for reactions and thermal profiles are provided in Appendix I. Each reaction used polymerase chain reaction buffer C (Invitrogen, Carlsbad, California), Go Taq polymerase (Promega, Madison, Wisconsin), and deoxynucleoside triphosphates (Invitrogen). Cleanup of polymerase chain reaction product and cycle sequencing were performed following Piaggio and Perkins (2005). Cycle sequencing cleanup was accomplished with PrepEase (USB, Cleveland, Ohio) protocols and resulting sequences were visualized on an ABI3130xl genetic analyzer (Applied Biosystems, Foster City, California).
At the Genetic Data Centre, 2 primer sets targeting the mtDNA CR were designed using conserved flanking regions: Cytb as the source of a forward primer (CytbARF 5′ GCCAGTTGAATACCCATTTATTGC 3′) and 12SRNA and Phe tRNA as sources for the reverse primers (12SAR2R 5′ GCATTTTCACTGGGGCGAGGAGT 3′ and PheARR 5′ GCTTTGCTTTATTTAAGCTAC 3′) resulting in polymerase chain reaction products approximately 1,300 bp in size. Details for reactions and thermal profiles are provided in Appendix I. Each used polymerase chain reaction buffer for Paq5000 (Agilent Technologies Canada Inc., Toronto, Ontario, Canada), Paq5000 polymerase (Agilent Technologies Canada Inc.), and deoxynucleoside triphosphates (New England Biolabs, Ipswich, Massachusetts). Cycle sequencing was performed following Franklin et al. (2011) using Sequitherm Excel II (Interscience, Markham, Ontario, Canada) and resulting sequences were visualized on a LiCor 4200 automatic sequencer (LiCor Biotechnology, Lincoln, Nebraska). A single CR sequence (AR 397) was generously donated by W. Zielenski and K. Pilgrim for this study (Table 1).
Sequence analyses
All sequences were aligned and edited in Sequencher (ver. 4.8; Gene Codes Corp., Ann Arbor, Michigan). Polymorphisms within GHR were coded with IUPAC nucleotide ambiguity codes. Each marker data set was assessed with ALTER (Alignment Transformation EnviRonment—Glez-Pena et al. 2010) to remove redundant haplotypes and compiled for downstream parsimony and likelihood analyses. Haplotype tables were compiled using DnaDiffer (Ritland 2012); haplotype tables can be accessed at Dryad Digital Repository (2012). Model selection was accomplished using MODELTEST version 3.6 (Posada and Crandall 1998) with Akaike information criterion (Akaike 1974). Total number of characters and parsimony-informative sites were calculated per data set using PAUP*4.0b10 (Swofford 2002). Insertions and deletions within sequence data sets were determined using DnaSP (Librado and Rozas 2009).
Phylogenetic analyses and divergence estimates
Phylogenetic relationships were estimated using maximum parsimony and maximum likelihood in PAUP*4.0b10 (Swofford 2002) and Bayesian methods in BEAST version 1.6.1 (Drummond and Rambaut 2007). Members of the Sciuridae were chosen as outgroup taxa for Cytb and GHR analyses (Appendix II; Wettstein et al. 1995; Kruckenhauser et al. 1999; Adkins et al. 2001; Piaggio and Spicer 2001; Harrison et al. 2003; Álvarez-Castañeda 2007; Huchon et al. 2007; Oshida et al. 2009; Kerhoulas and Arbogast 2010) because these have been determined to be sister taxa to Aplodontia (Huchon et al. 1999; Adkins et al. 2001; Herron et al. 2004). Because the squirrel CR sequences were too divergent and alignment with outgroup taxa was difficult, we performed the CR phylogenetic analysis with both an outgroup of samples of A. r. pacifica (see results of Cytb analysis) and also midpoint rooting. Initially, sequences from each marker (CR, Cytb, and GHR) were analyzed separately then we concatenated mtDNA data for a total mtDNA analysis. We performed the partition-homogeneity test in PAUP*4.0b10 (Swofford 2002) to test the validity of concatenating CR and Cytb sequences. This test was accomplished with 50 repetitions of a heuristic search, 100 tree-bisection-reconnection branch-swapping addition sequences each, and a reconnection limit of 8 and 10,000,000 rearrangements per repetition. For each data set, maximum-parsimony analyses consisted of heuristic searches with branch swapping by stepwise addition, 100 repetitions of random addition sequences, and tree-bisection-reconnection branch swapping with a reconnection limit of 8 and 10,000,000 rearrangements per replicate. All characters were unordered and had equal weights. Maximum-likelihood analyses were performed with estimates from the chosen evolutionary model through the initial generation of a neighbor-joining tree with tree-bisection-reconnection branch swapping as in the maximum-parsimony analysis. Nodal support was estimated using the nonparametric bootstrap (Felsenstein 1985) under maximum likelihood with 100 replicates and tree-bisection-reconnection branch swapping. Nucleotide diversity (Nei 1987:equation 10.5) within lineages and nucleotide divergences (Nei 1987: equation 10.20) between lineages were calculated in DnaSP (Librado and Rozas 2009). This was accomplished using primarily Cytb sequences but in some cases where unique Cytb haplotypes were not obtained (i.e., coastal California lineages) they had to be calculated from CR sequences.
The program BEAST version 1.6.1 (Drummond and Rambaut 2007) was used for Bayesian lineage divergence estimation using the most complete mtDNA data set (450 bp Cytb) and nDNA (GHR) analyzed separately. Initially a molecular clock was tested; when the molecular clock was rejected then divergence estimates were conducted under a relaxed clock uncorrelated lognormal (Drummond et al. 2006), and the appropriate evolutionary model. A coalescent exponential growth parameter was applied because most of the analyses involved within species relationships. We used fossil dates of 34–38 million years ago (mya—Vianey-Liaud 1985) for the Aplodontia–sciurid time to most recent common ancestor (tmrca) prior with a lognormal distribution. This closely matches the molecular estimated date for this divergence from mtDNA (Montgelard et al. 2002). This date served as an informative prior for the root of the tree. Only a few fossils of modern Aplodontia have been found and they all date to the Pleistocene. We chose one of these in an effort to balance out the much older prior placed at the root of the tree and hopefully achieve more robust posterior probabilities for estimates of divergences within this genus (Ho and Phillips 2009). Therefore, a 2nd fossil prior of 8,000 years ago was applied for the tmrca of the coastal California subspecies (Wake 2006) with an exponential distribution. Two other priors were set to a uniform distribution; one for the standard deviation of the uncorrelated lognormal relaxed clock (lower = 0, upper = 5) and the other for the mean of the branch rates (lower = 1E−5, upper = 3). We have little information about the values of standard deviation or the mean of the branch rates but we know that the default (0–∞) is not appropriate. Therefore, a uniform prior (noninformative prior) distribution is intended to allow Bayesian inference for parameters about which not much is known beyond the data included in the analysis at hand. BEAST was run multiple times to determine the most appropriate chain length. Finally, 3 independent runs were accomplished with 200 million generations of the Markov chain Monte Carlo analysis with sampling at every 1,000 generations. Convergence, effective samples sizes, and divergence times with upper and lower 95% highest posterior density (HPD) bounds were assessed in Tracer, version 1.5 (Rambaut and Drummond 2007) with the 3 runs considered separately to confirm consistency among runs. When consistency was confirmed, trees from 1 Cytb run were resampled at a lower rate in LogCombiner (6,000 states) to reduce the number of trees to be analyzed in TreeAnnotator. This was not necessary for GHR results. TreeAnnotator was used for both analyses with 10% of the samples removed for burn-in to generate a tree that was then visualized and edited for publication in FigTree, version 1.3.1 (Rambaut 2009).
Results
Sequence and trees statistics
Initially 114 representative samples (including 2 outgroup samples) from different geographical localities were sequenced for the whole Cytb (1,050 bp) and phylogenetic analyses were conducted. We found that phylogenetic trees generated from sequences of 1,050 bp and the first 450 bp of Cytb had identical topologies for cladogenesis among subspecies (Fig. 2); however, 1,050 bp of Cytb provided some further phylogenetic resolution within 1 clade (see inset in Fig. 2). Because we were concerned primarily with resolution at the deeper branches for taxonomic designations we decided to sequence the rest of samples for only the first 450 bp. Therefore, 380 samples were sequenced for this fragment (Fig. 2; Table 1). These samples represented complete taxon sampling for most subspecies (except A. r. californica). Subsets of these samples were selected for GHR and CR sequencing after initial phylogenetic analyses of Cytb data. Representative samples were selected from each clade; however, subsets were not the same for each marker. Among the 89 samples sequenced for CR there were 17 that were not included in Cytb sequencing: 16 from British Columbia and 1 from California (Table 1). All sequences obtained for each DNA fragment were deposited in GenBank (accession numbers JX419435–JX 420115). Appropriate models of evolution for each data set were assessed and relevant parameters were calculated or estimated (Table 2). Parsimony-informative sites, number of characters, and insertions and deletions also were calculated for each data set (Table 2).
Phylogenetic relationships
To generate phylogenetic trees, maximum parsimony and maximum likelihood were performed on Cytb, CR, and GHR data sets with all sequences and with nonredundant haplotypes, whereas Bayesian trees were inferred from Cytb and GHR with total sequence data (for access to data files, see Dryad Digital Repository [2012]). One concatenated data set with CR and Cytb nonredundant haplotypes (partition-homogeneity test, P = 0.75) was used to generate maximum-parsimony and maximum-likelihood trees. MtDNA phylogenetic trees generated under all approaches inferred the same evolutionary relationships, hence maximum-likelihood trees with nonredundant haplotypes are presented with bootstrap support and posterior probabilities for ease of viewing (Figs. 2 and 3). The evolutionary relationships inferred among subspecies in the concatenated mtDNA tree were analogous to the CR tree. The CR tree is presented because this marker alone contributed resolution to relationships among the coastal California subspecies (Fig. 3). Likelihood analyses of CR with 2 different rooting schemes (midpoint and outgroup) were identical; midpoint rooted trees are presented. The maximum-parsimony and maximum-likelihood GHR trees exhibited poor bootstrap support outside of separating the squirrels from Aplodontia. However, the Bayesian tree had greater resolution with posterior probability support for lineages within Aplodontia and therefore is presented (Fig. 4).
Five distinctive and well-supported clades within Aplodontia were formed from analysis of Cytb sequences (Fig. 2) under maximum-parsimony and maximum-likelihood analyses. A single clade, made up of samples identified when collected as A. r. pacifica, was quite divergent from the remaining lineages and its basal position was statistically well supported (>95%). This relationship was supported in the CR (maximum likelihood), GHR (Bayesian), and total mtDNA trees (Figs. 3 and 4). This suggested this lineage was distinct and unique from the rest of Aplodontia.
The relationships among the remaining lineages were not well supported in any analysis; however, the monophyly of each subspecies had high statistical support (Figs. 2–4). One exception is that there was support for relationships among the coastal California subspecies in the CR tree (Fig. 3); with A. r. humboldtiana and A. r. nigra as sister subspecies and A. r. phaea sister to those 2 clades. Within both mtDNA trees all samples of A. r. californica grouped into 1 lineage. However, this was not tested in the nDNA tree because we had only 1 sample of A. r. californica in that analysis. Samples from within the designated range of A. r. rufa fell into 2 separate lineages and this is statistically well supported in both mtDNA and nDNA analyses. One clade had samples from Washington and British Columbia and another had samples from areas of Oregon and northern California. No distinct phylogenetic difference was identified between purported A. r. rufa and A. r. rainieri throughout Washington and British Columbia. In fact, there appeared to be a single clade in this region (Figs. 2–4). The southern limit of the Washington and British Columbia clade and the northern limit of the Oregon and northern California clade is the Columbia River. In the Cytb tree, the 3 coastal California subspecies belonged to a single clade without support as unique lineages (Fig. 2). In fact, samples of A. r. nigra had a single Cytb and GHR haplotype (see haplotype data at Dryad Digital Repository [2012]) that was identical to a haplotype of A. r. humboldtiana (Figs. 2 and 4). However, in the CR phylogeny each subspecies had a unique haplotype(s) and the 3 subspecies formed monophyletic groups (Fig. 3). Therefore, each recognized coastal California subspecies occupied its own clade in the CR phylogeny.
Uncorrected pairwise sequence divergences estimated from the Cytb sequences within Aplodontia clades range from 0.004 to 0.01. The average sequence divergence between A. r. pacifica and the clade containing other Aplodontia subspecies was 0.13. Sequence divergences between the coastal California subspecies were calculated from the CR data because this was the only data set where each subspecies had unique haplotypes. Average divergences within these clades estimated from CR sequences ranged from 0.003 to 0.03. Divergences between the coastal California lineages ranged from 0.01 to 0.02. Between these clades and other Aplodontia clades, CR sequence divergences ranged from 0.03 to 0.04 with the exception of A. r. pacifica, where the divergences averaged 0.06, as they did between A. r. pacifica and all other Aplodontia clades.
Divergence time estimation
A strict molecular clock was rejected for the Cytb and GHR data sets (complete .xml files at Dryad Digital Repository [2012]) because the coefficient of variation estimation by BEAST was not close to zero. Therefore, the resulting divergence time estimate obtained from each data set was accomplished with the application of a relaxed clock. Posterior probabilities of divergences were not in agreement between Cytb and GHR (Table 3); nonetheless, the 95% HPDs for both markers were broadly overlapping. In fact, results from each marker indicated that the genus Aplodontia diverged about 1–2 mya. Further, divergence date estimates within Aplodontia occurred within the last million years. A. r. pacifica split from a common ancestor shared with other Aplodontia around 960,000 years ago (Cytb; Table 3). Among other clades, A. r. olympica (see “Discussion”) exhibited the earliest within-subspecies divergence and the most recent divergences were within the coastal California subspecies. There was significant among-lineage nucleotide substitution rate heterogeneity in both genes, as measured by coefficient of variation (Cytb = 2.26, 95% HPD: [1.00, 3.91]; GHR = 3.04, 95% HPD: [0.30, 5.61]). Posterior exponential growth was slightly negative for both genes (Cytb = −0.26, 95% HPD: [−0.40, −0.13]; GHR = −0.15, 95% HPD: [−0.28, −0.03]), suggesting a reduced effective population size.
Discussion
Phylogeny and taxonomy of Aplodontia
The inferred mtDNA and nDNA phylogenetic trees of Aplodontia closely reflect the subspecific designations that Hall (1981) described, with a few exceptions (Fig. 5). The subspecies designated as A. r. pacifica by Hall (1981) occupies its own well-supported clade (>95%), sister to all other lineages of A. rufa in both mtDNA and nDNA trees (Figs. 2–4). Notably, our samples included individuals from near the type locality of A. r. pacifica at Newport, Oregon (Merriam 1899). Sequence divergence of Cytb data between A. r. pacifica and the other subspecies range from 11.7% to 13.4%. Considering that the average Cytb divergence for mammalian species is 11% (Bradley and Baker 2001), these divergences suggest A. r. pacifica may be a separate species. In addition, when examining 450 bp of Cytb we identified 3 amino acid changes (42 bp Ala to Ile/S; 102 bp Met to Leu; and 110 bp Tyr to Ser) that are unique to the A. r. pacifica lineage relative to all other Aplodontia lineages (see haplotype data at Dryad Digital Repository [2012]). The morphology of A. r. pacifica is distinct from that of other subspecies of A. rufa: smaller body size, narrower defined rostrum, and larger, rounder ears (Merriam 1899). Karotypes of A. r. pacifica, A. r. californica, and A. r. phaea have been documented (Carrasco and Humphrey 1968; McMillin and Sutton 1972) and all have diploid chromosome numbers of 46. However, A. r. pacifica was described as having 5 pairs of metacentric and 17 pairs of submetacentric autosomes (Carrasco and Humphrey 1968), whereas the other 2 subspecies had 6 metacentric autosomes (McMillin and Sutton 1972). Further, the Y-chromosome of A. r. pacifica was described as metacentric, whereas the Y-chromosome of the other 2 subspecies was described as submetacentric (Carrasco and Humphrey 1968; McMillin and Sutton 1972). This difference alone is considered sufficient to warrant taxonomic distinction up to the subgeneric level (Nadler 1966a, 1966b) and is yet another piece of evidence that A. r. pacifica and A. rufa have had separate evolutionary histories for some time. Elevation of A. r. pacifica to species-level status would reinstate the classification of Merriam (1899), yet we acknowledge that such a revision requires a thorough morphological assessment. Specifically, further study is needed to verify that all of the unique morphological and karyotypic characteristics attributed to A. pacifica are consistent with our redefined geographic range for the subspecies. Likewise, a more detailed molecular assessment of samples from between the southern limit of our sampling of A. r. pacifica and the northern limit of our samples of A. r. rufa is needed to confirm distinction and a lack of a hybrid zone. Therefore, we retain the subspecific status of A. r. pacifica.
Individuals collected within the range of A. r. rufa (Hall 1981) were found in 2 well-supported clades in both mtDNA and nDNA trees. One clade had a geographical range in Oregon and northern California (Figs. 2–5). The other clade contained individuals from within the ranges of both A. r. rufa and A. r. rainieri (Hall 1981) north of the Columbia River (Figs. 1–5). The finding that A. r. rainieri and A. r. rufa fell into a single clade with individuals from Washington and British Columbia is not surprising because many researchers have noted that these subspecies are difficult to distinguish using morphological characters (Dalquest and Scheffer 1945; Ransome 2003). We targeted samples from the type locality of A. r. rainieri (Merriam 1899), yet we found no apparent genetic distinction of samples from this site or its proximity. The Washington and British Columbia clade can be identified as A. r. olympica based on Merriam's (1899) designation of samples from this region and Taylor's (1918) revision of the genus. The type locality for Aplodontia olympica was from “Queniult Lake, Olympic Mts., Washington” (Merriam 1899:20), which Taylor (1918) spelled “Quiniault” and Carraway and Verts (1993) found to be in Grays Harbor County, Washington. A great number of our samples in this clade were trapped from various localities in this county. We also suggest synonymizing A. r. rainieri Merriam, 1899, into A. r. olympica Merriam, 1899 (Fig. 5). The result would be that any A. rufa north of the Columbia River would be identified as a single subspecies, A. r. olympica, including individuals from British Columbia, Canada, where there are currently 2 separate conservation plans for the 2 presumptive subspecies (Ransome 2003). The clade from Oregon and northern California should retain the name A. r. rufa (Fig. 5) because the type locality for A. r. rufa is on the southern side of the Columbia River in Oregon (Richardson 1829b; Merriam 1899). Furthermore, a greater effort should be put into sampling A. r. rufa throughout its distribution; our sampling identified an apparent divergence between Willamette Valley and coastal Oregon individuals in the Cytb analysis (see inset in Fig. 2).
The current designation of California subspecies, A. r. nigra, A. r. humboldtiana, and A. r. phaea, as distinct subspecies is concordant with only the CR phylogenetic results. The distinction of A. r. nigra and A. r. humboldtiana was not supported by Cytb (nor GHR) analyses but A. r. phaea was well supported (Fig. 2). However, there are clear clades that are statistically well supported corresponding to each of these subspecies in the CR phylogeny (Fig. 3). This is likely because CR evolves more rapidly than Cytb and GHR and therefore can detect more-recent divergences (Randi 2000). Further, although the GHR Bayesian tree did provide support for some other clades within Aplodontia there were only 3 parsimony-informative sites in this data set. Also within California, our results support A. r. californica as a valid subspecies; however, our sampling was only in a very limited portion of their distribution (Fig. 1), as described by Hall (1981). Further study is required to test the geographic distribution of this subspecies.
Geographical distributions and implications for management and conservation
A map showing an approximation of revised geographic boundaries within the genus Aplodontia is presented (Fig. 5). This map was compiled with data from extensive sampling within the ranges of most subspecies. However, in some cases our sampling did not cover the entire range of a subspecies (e.g., A. r. californica) and therefore we can only be sure this lineage is represented in the shaded area (= sampling area), but this does not mean it is not found in other areas of its distribution as designated by Hall (1981); hence, we have left Hall's boundaries. In the future, addition of samples from these areas may expand this range map.
Our phylogeographic assessment has reduced the geographic distribution of A. r. rufa as compared to that of Hall (1981; Fig. 5); however, this should not affect the management of this subspecies because they are still widely distributed. A similar situation is found for A. r. pacifica. Through targeted sampling along the Oregon coast we found that the range of A. r. pacifica did not extend as far south as previously described (Hall 1981; Figs. 1 and 5); some samples from Oregon that fell into the A. r. rufa clade were from areas thought to be part of the range of A. r. pacifica. Again this should not change the management practices given they are still common throughout the revised range.
The combination of A. r. rainieri and northern A. r. rufa into A. r. olympica in Washington and British Columbia may have an impact on how mountain beavers are managed in British Columbia, Canada. Previously considered “vulnerable” and “critically imperiled, or imperiled” as separate subspecies, respectively (British Columbia Conservation Data Centre 2012; Gyug 2000), the status of A. rufa may be reduced to “vulnerable” or removed from the species of concern list once recognized as a single subspecies, thereby increasing its numbers and range. Because mountain beavers, even as 2 separate subspecies, were not previously protected in Washington, the management of A. r. olympica would be status quo in that state.
Unfortunately, our sampling of the range of A. r. californica is insufficient to address the limits of its distribution. The endangered species status for the subspecies A. r. nigra and species of special concern for A. r. phaea will not be altered by our findings because our results support the distinction of these lineages along with A. r. humboldtiana at the subspecific level. Notably, we found low diversity among samples of A. r. nigra with only 1 Cytb haplotype and 2 CR haplotypes (see haplotype data at Dryad Digital Repository [2012]). It is also noteworthy that divergences within the A. r. humboldtiana clade are at least as large as those between A. r. humboldtiana, A. r. nigra, and A. r. phaea together. Further, haplotype diversity of A. r. humboldtiana is relatively high, suggesting possible cryptic diversity, thus warranting further investigation. We sampled A. r. nigra and A. r. phaea from their type localities but not A. r. humboldtiana (Merriam 1899; Taylor 1914, 1916a). Clearly, continuing research and conservation efforts for the isolated, disjunct coastal California populations are critical for the maintenance of these unique lineages.
An examination of the geographic distributions of subspecies that we documented (Fig. 5) will undoubtedly raise questions about what possible geographic barriers explain the proximate distributions of A. r. olympica, A. r. rufa, A. r. humboldtiana, and A. r. pacifica. Pleistocene influences may explain some of the current geographical boundaries, yet there may be other influences. The current distribution of Aplodontia correlates to the distribution of the temperate broadleaf and mixed forests ecoregion in the Pacific Northwest region of the United States as defined by Ricketts et al. (1999). This explains the southern limitations of the distribution but not the northern range in British Columbia, where the ecoregion is much more broadly distributed than is Aplodontia. The northern limit of the genus and the remaining barriers among continuously distributed Aplodontia may be explained to some degree by the presence of large rivers. However, specifics of river size that limit mountain beavers need to be further defined. Gyug (2000) conducted extensive field surveys and investigations of historical records to determine the limits of A. rufa in British Columbia. Gyug (2000) determined that the Fraser River limited northward dispersal even with suitable habitat on the northern side of this river; this distribution limit was confirmed by a later study (Ransome 2003). The Columbia River also clearly delimits A. r. rufa and A. r. olympica. Finley (1941) first recognized the Columbia River as a barrier for Aplodontia but he believed it was an effective barrier only in its lower reaches and that further upriver it was less effective. Steele (1989:15) stated that A. r. humboldtiana is “geographically isolated on three sides by the Smith River to the north, the coastline to the west, and a distance of approximately 160 km of unoccupied territory to the south.” We collected an A. rufa from the Jedediah Smith State Park in Del Norte County, California, which is on the other side of the Smith River (northwest) from the distribution of A. r. humboldtiana and our sample from this park is clearly an A. r. rufa (Table 1), confirming that the Smith River is a probable barrier. Without an obvious explanation for the phylogeographical split between A. r. pacifica and A. r. rufa, one must wonder if the Siuslaw River and its tributaries, the Alsea River, or the Yaquina River may serve as a boundary for A. r. pacifica. This could be investigated through the identification of the southernmost population of A. r. pacifica. Nonetheless, we conclude that rivers can be barriers to dispersal of Aplodontia because they delimit some current geographical boundaries among lineages within this genus.
Dates of divergences
Although powerful analytical tools have been developed recently for divergence date estimation and phylogenetic assessment (Drummond et al. 2006), there are still limitations to the inferences that can be made from some data sets. In our case, the estimates from mtDNA and nDNA are broadly overlapping, with mtDNA providing more resolution than nDNA; we regard these estimates as providing a general timeline for the divergences of Aplodontia. Therefore, we conclude from both lines of evidence that modern Aplodontia arose <5 mya (probably 1–2 mya) and the divergences within this genus occurred during the Pleistocene. The divergence estimate of 1–2 mya for modern Aplodontia closely matches the estimate for this divergence obtained from phylogenetic analysis of fossil data (Hopkins 2008).
There is a large gap in fossil data for aplodontiids from 6 mya to about 70,000 years ago (Shotwell 1958; Hopkins 2007). Interestingly, the fossils that date to about 70,000 years ago (Kurtén and Anderson 1980) are all from 3 northern California caves and this time period approximates our estimated divergences (Cytb) of A. r. rufa and A. r. californica (Table 3), both of which currently occur in California. Shafer et al. (2010) suggest the occurrence of multiple refugia in the Pacific Northwest during the Pleistocene climatic oscillations. The southern boundary of the ice sheets occurred in this region and it is clear that the climatic oscillations of this time period shaped the evolutionary relationships of the region's biota (Shafer et al. 2010), including Aplodontia, and our estimates of divergence times within this genus confirm this impact. Our results also suggest a reduction in overall demography through time, which could be due to anthropogenic or climatic alteration, or both, of habitat of Aplodontia, as reflected by the endangered and threatened California subspecies. This issue requires further investigation before final conclusions are reached.
The genus Aplodontia of the family Aplodontiidae was understood to be a monotypic genus with the species A. rufa. One objective of this study was to use a molecular phylogenetic approach to test the morphological taxonomy (Hall 1981) as a hypothesis of evolutionary relationships. Using this approach, this study has identified 2 possible species-level divergences within this genus, which may require the resurrection of a historically recognized species, pending further investigation. Also, by incorporating geographical information the current location and distribution of subspecies of A. rufa were elucidated. As a result, this study identified a single subspecies (A. r. olympica) in Washington and British Columbia, Canada, where previously 2 subspecies were designated. Further, distributions of subspecies in Oregon were revised and taxonomic designations and geographical distributions of coastal California subspecies of A. rufa were confirmed. Other than possible management changes for A. r. olympica, the conservation status for all other subspecies of Aplodontia should remain as status quo. Finally, our estimated divergence dates within Aplodontia demonstrate Pleistocene influences and thus provide a better understanding of evolution within this group.
Aplodontia Richardson, 1829b
Aplodontia Richardson, 1829b:334. Type species Aplodontia leporina Richardson, 1829b:335 (= Anisonyx rufa Rafinesque, 1817:45).
Aplodontia rufa Rafinesque, 1817
Anisonyx? rufa Rafinesque, 1817:45. Type locality “Neighbourhood of the Columbia River.”
Aplodontia rufa: Merriam, 1886:316. First use of current name combination.
Aplodontia rufa californica Peters, 1864
H[aplodon]. leporinus var. Californicus Peters, 1864:179. Type locality Sierra Nevada of California (Merriam 1886; Taylor 1918).
Aplodontia major Merriam, 1886:316. Type locality “Sierra Nevada Mountains, in Placer County, Cal.”
[Aplodontia rufa] californica: Trouessart, 1904:348. First use of current name combination.
Aplodontia rufa humboldtiana Taylor, 1916a
Aplodontia humboldtiana Taylor, 1916a:21. Type locality “Carlotta, Humboldt County, California.”
Aplodontia rufa humboldtiana: Taylor, 1918:470. First use of current name combination.
Aplodontia rufa pacifica Merriam, 1899
Aplodontia pacifica Merriam, 1899:19. Type locality “from Newport, mouth of Yaquina Bay, Oregon.”
Aplodontia rufa pacifica: Taylor, 1918:480. First use of current name combination.
Aplodontia rufa rufa Rafinesque, 1817
Anisonyx? rufa Rafinesque, 1817:45. Type locality “Neighbourhood of the Columbia River.”
Aplodontia leporina Richardson, 1829b:335. No type locality, but near Columbia River implied.
Aplodontia chryseola Kellogg, 1914:295. Type locality “Jackson Lake, Siskiyou County, California altitude 5900 feet.”
Aplodontia rufa grisea Taylor, 1916b:497. Type locality “Renton [near Seattle], Washington.”
Aplodontia rufa nigra Taylor, 1914
Aplodontia nigra Taylor, 1914:297. Type locality “Point Arena, Mendocino County, California.”
Aplodontia rufa nigra: Taylor, 1918:479. First use of current name combination.
Aplodontia rufa olympica Merriam, 1899
Aplodontia olympica Merriam, 1899:20. Type locality “Queniult Lake, Olympic Mts., Washington.”
Aplodontia major rainieri Merriam, 1899:21. Type locality “Paradise Creek, south side Mt. Rainier, Washington (alt., 5200 ft.).”
Aplodontia californica columbiana Taylor, 1916b:499. Type locality “Roab's Ranch, Hope, British Columbia.”
Aplodontia rufa olympica: Taylor, 1918:460. First use of current name combination.
Aplodontia rufa phaea Merriam, 1899
Aplodontia phaea Merriam, 1899:20. Type locality “Pt. Reyes, Marin Co., California.”
Aplodontia rufa phaea: Taylor, 1918:480. First use of current name combination.
Acknowledgments
We thank W. Zielinski, M. Schwartz, and K. Pilgrim of the United States Forest Service for their kindness, generosity, and collaboration. We also extend our gratitude to T. Best for helpful discussions on taxonomic designations. We also thank M. A. Neubaum and J. Figueroa; J. Jeffers of the Nevada Department of Wildlife; and D. Girman and D. Crocker from Sonoma State University for their assistance in this project, including, in some cases, donations of samples. We also acknowledge K. Hamm, B. Youst, D. Lancaster, F. Arnold, B. Covell, J. Carr, R. Best, E. Meister, J. Duvall, J. Yost, J. Schaberl, T. Kocket, T. Losli, P. Happe, G. Phelan, R. Gonzalez, M. Dykzeul, and E. Myers for their help in obtaining samples. Thanks to J. Fillon and K. Ritland for their invaluable assistances in haplotype submission and discernment, respectively. We appreciate the support of the sabbatical program of the United States Department of Agriculture, Wildlife Services, National Wildlife Research Center; the Washington Forest Protection Association; and Oregon Forest Industries Council. Finally, we acknowledge the helpful comments and guidance from an anonymous reviewer.
Literature Cited
Appendices
Appendix I
Appendix I
Polymerase chain reaction chemistry and conditions for each primer set targeting 3 DNA regions. dNTP = deoxynucleoside triphosphate; °C/TID = temperature in degrees Celsius/initial denaturation time in minutes; °C/TD = temperature in degrees Celsius/initial denaturation time in seconds; °C/TA = temperature in degrees Celsius/annealing time in seconds; °C/TE = temperature in degrees Celsius/extension time in minutes; °C/TFE = temperature in degrees Celsius/final extension time in minutes; Cytb = cytochrome b; GHR = exon 10 of the growth hormone receptor; CR = control region.
Appendices
Appendix II
Appendix II
Outgroup samples sequenced and analyzed in this study (Hall 1981), localities by state and county, collector identification, loci sequenced per sample, and GenBank accession numbers indicated. Cytb = cytochrome b; CR = control region; GHR = exon 10 of the growth hormone receptor.
Fig. 1
Map of distribution of Aplodontia rufa adapted from Hall (1981:map 233). Circles represent samples collected for this study. Circle size correlates to the number of samples collected. Some sampling locations are proximate to one another and therefore caused overlap of some circles or appear as a single circle in some cases.
Fig. 2
Maximum-likelihood phylogram inferred from nonredundant haplotypes from 380 samples sequenced for the first 450 bp of cytochrome b (Cytb). This tree is the same topology as a tree inferred from nonredundant haplotypes from 112 samples sequenced for 1,050 bp of Cytb with the exception of phylogeographic resolution within Aplodontia rufa rufa (see inset with bootstrap support). These 2 data sets have different samples representing the tips of the trees but infer the same clades and overall topology. Bootstrap support ≥ 70 is indicated with ○, whereas support ≥ 90 is indicated with •. The 1st symbol represents bootstrap support estimated from the first 450 bp of Cytb and the 2nd symbol is the bootstrap from the 1,050 bp of Cytb. The 3rd symbol represents the posterior probability from the Bayesian analysis ≥ 0.95. Outgroup taxa relationships are collapsed in this figure for ease of presentation. State(s) where samples were collected are designated in parentheses.
Fig. 3
Maximum-likelihood phylogram inferred from nonredundant haplotypes from 89 samples sequenced for the control region (CR). Bootstrap support ≥66 is indicated with ○, whereas support ≥90 is indicated with •. The 1st symbol is the bootstrap from maximum-likelihood analysis of the CR sequences and the 2nd symbol is from maximum-likelihood analysis of concatenated CR and cytochrome-b nonredundant haplotypes, which produced an identical topology. Revised subspecies designations are labeled on the tree. State(s) where samples were collected are designated in parentheses. We inferred the same topology with or without outgroup taxa so outgroup taxa were removed from this tree after analysis for ease of viewing and this tree is midpoint rooted.
Fig. 4
Nuclear DNA phylogeny inferred from Bayesian analysis of 35 growth hormone receptor exon 10 sequences. The symbol • represents the posterior probability from the Bayesian analysis ≥ 0.95. Revised subspecies designations are labeled on the tree. State(s) where samples were collected are designated in parentheses. The branches to the outgroup have been cut to reduce the overall size of the tree.
Fig. 5
Map adapted from Hall's (1981:map 233) distribution of Aplodontia rufa and revised based on the results of this study. Shaded areas indicate geographical areas within the range defined by Hall (1981) where we sampled and are confident that subspecies exist. Nonshaded areas remain as areas to be confirmed as a portion of the range and taxonomic association.
Table 1.
Samples used in this study per subspecies of Aplodontia rufa (Hall 1981), number per collection locality (n), sample identification (ID), collection location by state and county (Co.) or Canadian province, collector identification (CI), number of samples sequenced per locus with CB = 1,050 bp cytochrome-b (Cytb) gene, CB1 = 450 bp Cytb gene, CR = control region, and GHR = exon 10 of the growth hormone receptor.
Table 2.
Statistics and model parameters for sequence data. n = number of samples sequenced without outgroup taxa; nnrh = number of nonredundant haplotypes; bp = base pairs; pi = parsimony-informative sites without outgroup taxa; indels = insertions–deletions; ML = maximum likelihood; MP = maximum parsimony; CI/RI = confidence index/retention index; inv. sites = number of invariable sites; α = gamma shape parameter.
Table 3.
Bayesian posterior age estimates for Aplodontia using 2 different genes. TMRCA = time to most recent common ancestor; ka = thousand years ago; Cytb = cytochrome b; HPD = highest posterior density; GHR = exon 10 of the growth hormone receptor.
