Molecular Markers Reveal Limited Population Genetic Structure in a North American Corvid, Clark’s Nutcracker (Nucifraga columbiana)

The genetic impact of barriers and Pleistocene glaciations on high latitude resident species has not been widely investigated. The Clark’s nutcracker is an endemic North American corvid closely associated with Pinus-dominated forests. The nutcracker’s encompasses known barriers to dispersal for other species, and glaciated and unglaciated areas. Clark’s nutcrackers also irruptively disperse long distances in search of pine seed crops, creating the potential for gene flow among populations. Using the highly variable mitochondrial DNA control region, seven microsatellite loci, and species distribution modeling, we examined the effects of glaciations and dispersal barriers on population genetic patterns and population structure of nutcrackers. We sequenced 900 bp of mitochondrial control region for 169 individuals from 15 populations and analysed seven polymorphic microsatellite loci for 13 populations across the Clark’s nutcracker range. We used species distribution modeling and a range of phylogeographic analyses to examine evolutionary history. Clark’s nutcracker populations are not highly differentiated throughout their range, suggesting high levels of gene flow among populations, though we did find some evidence of isolation by distance and peripheral isolation. Our analyses suggested expansion from a single refugium after the last glacial maximum, but patterns of genetic diversity and paleodistribution modeling of suitable habitat were inconclusive as to the location of this refugium. Potential barriers to dispersal (e.g. mountain ranges) do not appear to restrict gene flow in Clark’s nutcracker, and postglacial expansion likely occurred quickly from a single refugium located south of the ice sheets.


Introduction
Phylogeography is the study of processes underlying spatial and temporal dimensions of genetic variation [1]. Patterns of genetic variation attained through phylogeographic methods can provide insights into population-level response to disturbance and the processes responsible for historical dispersal and colonization [2,3]. Examining patterns of genetic variation can help us to evaluate the roles of gene flow, bottlenecks, and historical or geological barriers in explaining contemporary patterns of diversity across geographical regions [4]. Contemporary genetic patterns are strongly influenced by postglacial expansion from refugia [5,6], historical and contemporary barriers to dispersal [7,8], and dispersal potential [9].
The repeated glaciations and climate oscillations of the Pleistocene epoch provide a natural tool to address population responses to large-scale landscape changes. Multiple expansions and contractions of ice sheets created a dynamic landscape that repeatedly fragmented historical populations, alternately creating barriers to dispersal and creating new habitats for colonization [10,11]. During the Pleistocene glacial maxima in North America, many populations were confined to refugia (ice-free areas) [10]. Plant and animal species expanded from several known refugia following the retreat of the ice sheets, including Beringia (present-day Bering Sea and parts of Alaska) and three main areas south of the ice sheets (Pacific Coast, Rockies, and Taiga) [10,12]. Additional areas along the periphery of the ice sheets are contested to have been ice-free, such as an Atlantic shelf near present-day Newfoundland [5,10,13].
High latitude resident bird species provide a unique opportunity to investigate patterns of postglacial and barriermediated dispersal; historical events shaping current population structure should be particularly evident in these species. Non-migratory species have the potential to retain patterns of genetic variation longer due to limited dispersal, allowing researchers to make inferences about past historic events [14][15][16]. For example, species of trees show distinct patterns of population genetic structure and have retained information on historic environmental changes [13,15]. As the number of studies on resident species increase, similar patterns are emerging in vertebrate taxa [16][17][18]. However, other resident species show limited differentiation despite potential barriers to dispersal [19,20].
Clark's nutcracker (Nucifraga columbiana) is a high latitude resident corvid species (Family Corvidae) that inhabits North American coniferous forests. Much of the published work on this species has focused on spatial cognition exhibited by their extensive food caching behaviours [21]. Additional research has concentrated on Clark's nutcracker's role as essential seed dispersers for many conifers, particularly Pinus species [22][23][24]. To date, no published phylogeographic work exists for this species, though genetic work has been identified as a priority [21]. Clark's nutcracker prefers higher altitude montane forests dominated by one or more Pinus species. Their close association with pine-dominated forests is due to a specialised diet of pine seeds, which influences dispersal patterns of both nutcrackers and seeds [21,25]. After nestlings fledge, Clark's nutcrackers have been known to undergo seasonal altitudinal movement to nearby Pinus-dominated subalpine areas [21]. In times of food shortage, and occasionally during the nonbreeding season, large numbers of birds have been shown to irruptively disperse, leaving home ranges, and travelling more than 100 km in search of large crops of pine seeds [21,26]. Range-wide morphological and plumage variation is well documented for some high latitude resident species [27,28], but there is little evidence of morphological variation in Clark's nutcracker [21]. Some suggest that this lack of variation is a reflection of high gene flow between populations due to irruptive dispersal, but this has not been investigated to date [21,29].
To address these gaps in phylogeographic knowledge, our goals for this study were to investigate patterns of dispersal and postglacial colonization in Clark's nutcracker as reflected in genetic structure. We predict that contemporary Clark's nutcracker populations will exhibit low levels of genetic differentiation between populations due to irruptive food-related dispersal patterns. Mountain ranges and unsuitable habitats will not act as strong barriers to gene flow due to this species' preference for high altitude sub-alpine habitats [21] where barriers are limited. Rather, dispersal will be restricted by availability of pine-dominated coniferous forests and seed crops. We expect that postglacial colonization likely occurred from a single refugium south of the ice sheets, similar to patterns of postglacial colonization in several North America Pinus species [13].

Ethics statement
All animals captured in the field were handled following animal welfare protocols (#0614 and #1028) approved by the University of Lethbridge Animal Welfare Committee using guidelines set by the Canadian Council on Animal Care (CCAC). Banding in Canada was performed under Canadian Wildlife Service banding permit #10425W (2007) and #10804 (2008 onward) and in the US under US Fish and Wildlife banding permit #23522. All netting and blood sampling was conducted under the appropriate state, provincial, federal, and institutional authorities; detailed permit information is available upon request.

Sample collection and preparation
From 2009-2012, we captured up to 30 individual Clark's nutcrackers at each sampling site (Table 1) using standard mistnetting techniques with call playback. For call playback, we used a medley of group interaction vocalizations sourced from multiple tracks found on xeno-canto.org. Sample collection locations for each defined population (i.e. sampling site) were limited to within a 50 km of radius and were not separated by any obvious barriers to dispersal. We chose sites on either side of possible barriers to dispersal (e.g., Rocky Mountains, Great Basin) and from areas that were previously glaciated and unglaciated during the last glacial maximum (e.g., Central Alberta and New Mexico, respectively; Figure 1). We collected less than 100 μL of blood from each bird by pricking the brachial vein with a sterile needle and collecting blood in a capillary tube. Blood was stored in 95% ethanol and archived at -80°C upon return to the University of Lethbridge. Each bird was banded with a US Fish & Wildlife aluminum band, aged and sexed (if possible), and mass, tarsus, and other morphological measurements taken. Additional tissue and feather samples were obtained from museum collections taken from birds during the breeding season within the past 20 years to ensure all samples were from contemporary populations (Table S1).
Laboratory procedures DNA extraction. DNA was extracted from blood, feather, and tissue samples using a modified Chelex protocol [30]. Once extracted, DNA was stored at -20°C in low TE buffer.
Polymerase chain reactions (PCR) were performed on a Master gradient thermocycler (Eppendorf: Hauppauge, NY) in a 25 μL volume containing 1x goTaq Flexi buffer (Promega: Madison, WI), 2.5 mM MgCl 2, 200 μM dNTP, 0.4 μM of each primer, and 0.5 units goTaq Flexi taq polymerase (Promega) under the following conditions: one cycle of 94°C for 120 s, 52°C for 45 s, and 72°C for 60 s, 37 cycles of 94°C for 30 s, 52°C for 45 s and 72°C for 60 s and one cycle of 72°C for five min. PCR products were run on a 0.8% agarose gel to confirm DNA amplification.
DNA sequencing was performed at McGill University and Génome Québec Innovation Centre on a 3730xl DNA Analyzer (Applied Biosystems: Carlsbad, CA) or at the University of Lethbridge on a 3130 DNA Analyzer (Applied Biosystems). For in-house sequencing we used a shrimp alkaline phosphataseexonuclease clean up followed by sequencing and sodium acetate precipitation [28] prior to gel electrophoresis.
Microsatellite DNA. Three individuals from geographically distant populations (Montana, southern Oregon, and Colorado) were initially screened for 32 microsatellite primer pairs developed for and used in other corvids. Thirteen additional nutcracker-specific microsatellite loci have since been developed by another independent research group [32]. If screened loci were suspected to be polymorphic and/or amplified clearly, two additional individuals from each of the initial test populations plus three individuals from another population potentially separated by dispersal barriers were screened. Eight loci were polymorphic ( Table 2). All forward primers were modified by adding an M13 sequence (5' -CAC   (Table 2). PCR products were mixed with a stop solution (95% formamide, 20 mM EDTA and bromophenol blue), denatured for 3 min at 94°C, and run on a 6% polyacrylamide gel using a LI-COR 4300 DNA Analyzer (LI-COR Inc.: Lincoln, NE). Alleles were scored via visual inspection, and genotypes were independently confirmed by a second person. Three controls of known allele sizes (pre-screened individuals) plus a size standard were included on each gel to ensure consistent scoring along with a negative control to ensure no contamination was present.

Analyses of genetic structure
Mitochondrial DNA. Chromatograms were edited and sequences aligned manually in MEGA v5.0 [33]. Shared haplotypes were determined using DnaSP v5.10 [34]. We assessed population structure and relationships among haplotypes using a maximum parsimony network constructed using Network v4.611, which is designed to construct the shortest, least complex network [35]. The median-joining (MJ) tree algorithm was used as this allows for multi-state data (e.g., four nucleotide base options at each variable site), with settings Clark's nutcracker population sampling locations. Current range overlaid in yellow [76]. Darker shades imply higher elevation areas. See Table 1 for coordinates, sample sizes, and abbreviations per population. doi: 10.1371/journal.pone.0079621.g001 as follows: weight = 10 for each site, epsilon = 0, and the MJ squared option turned off. After the tree was constructed, the MP (maximum parsimony) option was used to remove superfluous links and unnecessary median vectors from the network, which can be produced during the initial network calculations [36].
Genetic variation within populations was measured by calculating haplotype (H d ) and nucleotide (π) diversity using DnaSP v5.10 [34]. We calculated pairwise Φ ST values (an analogue of Wright's fixation index (F ST ) using Arlequin v3.5.1.3 [37] to examine population structure and test for genetic differentiation among populations and haplogroups. Significance values were corrected using a Benjamini-Hochberg correction [38], to control for false discovery rate (FDR), the expected proportion of falsely rejected hypotheses in multiple significance testing situations where multiplicity might occur.
A Mantel's test was performed to examine the correlation between genetic and geographic distances in GenAlEx v6.0 [39]. We calculated geographic distances using weighted central coordinates for each population (Table 1)   significance was tested using 9,999 permutations. A mismatch distribution was used to help visualize signatures of demographic expansion, and to test the null hypothesis of historic population growth and expansion [41]. Genetic variation allocated within and among populations was examined using an analysis of molecular variance (AMOVA) in Arlequin v3.5.1.3 [37]. We investigated the potential impact of barriers on genetic variation by performing a spatial analysis of molecular variance (SAMOVA; K = 2-12; 100 iterations) [42]; however, the process cannot test for rangewide panmixia as K cannot be set to K = 1. Alternatively, BAPS v6.0 [43], a Bayesian method for analysing population structure, can detect range-wide panmixia. A preliminary population mixture analysis for spatial clustering of individuals in BAPS used no a priori population information and K max = 13. A subsequent BAPS run was then conducted using a priori population information.
Microsatellite DNA. Allele frequencies, deviation from Hardy-Weinberg equilibrium (observed (Ho) versus expected (He) heterozygosity), and pairwise F ST values [44] were calculated with 999 permutations using GenAlEx v6.0 [39] and P values were corrected for multiple tests using a Benjamini-Hochberg correction for FDR [38]. Allelic richness was calculated in FSTAT v2.9.3 [45]. A Mantel's test was performed to test for isolation by distance using the same procedure as for mitochondrial DNA.
Bayesian clustering analyses were conducted using the programs STRUCTURE v2.3.3 [46,47] and BAPS v6.0 [48]. STRUCTURE was run with the following settings: K = 1-15, a burn-in of 100,000 followed by 500,000 runs, admixture assumed, correlated allele frequencies and including a priori population information 10 replicates were performed for each value of K. In STRUCTURE, it can be difficult to decide when K captures major structure in the data due to similar lnP(X|K) values, thus two additional analyses were conducted using STRUCTURE results: Bayes factor calculations [46] and ΔK [49] using STRUCTURE HARVESTER [50]. One limitation with using STRUCTURE HARVESTER is that it cannot detect when K = 1 is the true number of clusters. However, Bayes factors and BAPS are able to detect if K = 1, so the addition of these methods is useful in species where panmixia is a possibility. A preliminary population mixture analysis for spatial clustering of individuals in BAPS used no a priori population information and K max = 13. A subsequent BAPS run was then conducted using a priori population information.

Species distribution modeling
Sample points. In addition to our sample (field and museum) locations, geo-referenced Clark's nutcracker locations were obtained from the Global Biodiversity Information Facility (GBIF; http://data.gbif.org/, accessed on 22 January 2013). We excluded occurrences outside of North America, without geo-references, or recorded before 1950 from modeling and further inspected accuracy by plotting points using ArcMap v10.

Genetic analyses
Mitochondrial DNA. Mitochondrial DNA control region sequences (n = 169; Table S1) from 15 populations range-wide revealed 48 variable sites in a 900 base pair (bp) region. In total, 68 haplotypes were recognized in the median-joining analysis (GenBank accession nos. KF687612-KF687679; Figure 2; Table 3). Eighteen haplotypes accounted for 119 individuals; 50 singleton haplotypes were found (Figure 2; Table 3). The largest haplotype group contained a total of 50 individuals from 14 of the 15 populations analysed. Most haplotypes were a single mutational step removed from each other. Overall, limited geographic clustering was observed in the maximum parsimony network (Figure 2).
Haplotype diversity (H d ) values ranged from 0.583 (southern Alberta (SAB)) to 0.964 (Montana (MT), New Mexico (NM), and central California (CCA)) and nucleotide diversity (π) ranged from 0.00127 (SAB) to 0.00402 (northeastern Oregon) for populations with greater than five samples (Table 1). When all samples were combined, overall H d was 0.854, and overall π was 0.00324 indicating high levels of genetic diversity in Clark's nutcracker. In 78 population pairwise comparisons of genetic differentiation, Φ ST values ranged from -0.050 (P corrected = 0.779) for Utah (UT) and northeast Washington (NEWA) to  Table 3). Small black circles represent inferred nodes. See Table 1 for population abbreviations. , and NEWA). The mismatch distribution did not deviate significantly from the expected signature, so we cannot reject the null hypothesis of historical population growth and expansion [41]. Geographic distance in 78 paired populations ranged from 339 km for WY and UT to 2074 km for NM and central Alberta (CAB). A low, but significant correlation between genetic and geographic distance was detected for mitochondrial DNA (R 2 = 0.082, P = 0.040; Figure 3). When the NM population was removed from the analyses, the relationship between genetic and geographic distance still trended towards isolation by distance, though it was no longer significantly correlated (R 2 = 0.064, P = 0.058; not shown).
Both AMOVA and SAMOVA showed that most genetic variation is contained within rather than among population groups. The highest value of between-group variation was observed when K = 1 for AMOVA analyses (3.4% between groups, 96.6% within groups) and for K = 2 in SAMOVA (17.8% between groups, 80.7% within groups), with one group consisting of NM and the other of all other populations pooled. The preliminary (no a priori population information) and secondary (using a priori population information) BAPS analyses for spatial clustering of individuals found no significant differentiation between individuals or populations (i.e., K = 1). (Table 5) were used for analyses; an eighth locus (ApCo37) was removed due to limited amplification success and subsequent missing data (~50%). We used the thirteen populations with greater than five samples for microsatellite analyses (Table 1; Table S1). The total number of alleles detected ranged from four in PnuA106W to 10 in ApCo40 with an average of 6.43. Only 13 of the 91 loci-population comparisons showed significant deviations from Hardy-Weinberg and all showed lower heterozygosity than expected. ApCo41 has a significant heterozygote deficit in five of the 13 populations sampled and while two of these (NEWA and CAB) had fewer samples (n ≤ 8), the others had a larger number of samples (n ≥ 14; Table 5).

Microsatellite DNA. A total of seven polymorphic loci
Limited population differentiation was detected with a global F ST value of 0.070. Paired F ST values ranged from 0.000 (14/78 comparisons) to 0.123 (CCA and central Alberta (CAB); Table  6). After Benjamini-Hochberg FDR correction, 21 of 78 paired comparisons were significant (Table 6). A weak, but significant, correlation was found between genetic and geographic distance using Mantel's test (R 2 = 0.078, P = 0.029; Figure 4). STRUCTURE results suggested that the optimal number of groups was K = 1 ( Figure 5A). Further analyses produced peaks at K = 1, K = 4, and K = 12 ( Figure 5B), though K = 1 produced the highest peak and highest posterior probabilities (mean (Pr (K=1)) = -2625.45) and is supported by Bayes factor values (Pr (K=1) = 1) and BAPS results.
Species distribution modeling. A total of 1563 presence records were used for training and 174 for testing the distribution models. After removing correlated bioclimatic layers, 10 layers were used for model construction (bio1-4, 8,  Haplotype codes correspond to those found in Figure 2. Population abbreviations and locations are found in Table 1. doi: 10.1371/journal.pone.0079621.t003 12,[14][15][18][19]. Mean area under the receiver operating curve (AUC) was 0.907 (range = 0.906 -0.908), suggesting that the current model was highly suitable for backcasting to the paleodistribution model. Isothermality (mean diurnal temperature range/mean annual temperature range (bio3)), temperature seasonality (bio4), and mean annual temperature (bio1) were the most important variables contributing to the model at 44.3%, 37.7%, and 11.1% respectively. The current distribution predicted by the Maxent model ( Figure 6A) closely approximated the known Clark's nutcracker range in North America (Figure 1). The model returned a maximum 66% probability of suitable habitat throughout the range. The paleodistribution model returned a maximum 25% probability of predicted suitable habitat, mostly concentrated in the western mountainous regions and along the southeast edge of the ice sheets during the last glacial maximum ( Figure  6B).

Population genetic structure
Using a combination of measures of population differentiation and clustering analyses, our study suggests high levels of population connectivity between most Clark's nutcracker populations in North America. We found limited differentiation in mitochondrial DNA (Table 4), with only the New Mexico population showing significant pairwise differences from four other populations. Microsatellite data indicate higher levels of population differentiation between Clark's nutcracker populations, though this predominantly occurs between populations at the edges of the species' range (e.g. CAB and SCA; Table 6). Geographic distance rather than barriers appear to influence mitochondrial and nuclear gene flow, as there is a weak, but significant effect of isolation by distance on nutcracker populations for both molecular markers ( Figures 4 and 5). However, for mitochondrial DNA, this pattern seems to be driven by populations in New Mexico, suggesting limited peripheral isolation for this marker. For microsatellite data, additional peripheral populations (CAB, SAB, and SCA) are significantly differentiated in multiple comparisons ( Table  6), suggesting peripheral isolation in addition to isolation by distance. Clustering analyses do not support strong population differentiation for either nuclear or mitochondrial markers, with a single population being the best fit for the data in both cases. In the case of microsatellite data, these results should be interpreted with caution as both STRUCTURE and BAPS may struggle with situations where differentiation is low (F ST < 0.030) [56], though global values (F ST = 0.070) suggest this is likely not the case.
High levels of haplotype, nucleotide, and allelic diversity throughout most nutcracker populations suggest large population sizes and few founder effects or population bottlenecks. The exception to this is the southern Alberta population, which exhibits reduced haplotype and nucleotide diversity for mitochondrial markers relative to the rest of the populations, though microsatellite loci have high levels of allelic richness and heterozygosity (Table 5). This difference between marker genetic diversity may be due to a historical founder effect in previously glaciated southern Alberta (mtDNA), with nuclear gene flow being less limited (microsatellites).
Patterns of potential panmixia and limited population differentiation are not unique to Clark's nutcracker. Many other species that undergo irruptive and seasonal dispersal also exhibit limited population differentiation despite the presence of proposed barriers to dispersal. Studies of Pinus seed specialists, such as pygmy nuthatch (Sitta pygmaea) [57], crossbill species (Loxia spp.) [19], and Clark's nutcracker's sister species, Eurasian nutcracker (Nucifraga caryocatactes) [58], found no significant genetic differentiation among populations despite potential barriers to dispersal. Limited    differentiation has been found in other species with irruptive or seasonal dispersal, such as boreal owl (Aegolius funereus) [59], and European redpoll species (Carduelis spp.) [60], yet Steller's jay (Cyanocitta stelleri) [14] and white-breasted nuthatch (Sitta carolinensis) [61] exhibit significant differentiation between populations and limited gene flow corresponding with barriers to dispersal. In other taxa, peripheral populations are more likely to be isolated due to reduced gene flow [16,62,63], while maintaining panmixia in the majority of the range [64]. Overall, Clark's nutcracker population connectivity seems to be primarily limited by distance with some peripheral isolation at the southern and northern edges of the range, which is reasonable for a high elevation species with irruptive and seasonal altitudinal dispersal.

Postglacial colonization
After the LGM, species that rapidly expanded from a single refugium with high levels of gene flow are expected to exhibit three characteristics: 1) a star-like phylogeny with many low Overall A n = 10 A n = 5 A n = 10 A n = 7 A n = 7 A n = 7 A n = 4 Only populations with greater than five samples were used; n = number of samples used in genotyping and analyses; A n = number of alleles; A r = allelic richness; Ho = observed and He = expected heterozygosity; P = departures from Hardy-Weinberg equilibrium (ns = not significant, *P < 0.05, **P < 0.01, ***P < 0.001. See frequency single haplotypes separated from high frequency central ancestral haplotypes by few mutational steps [65]; 2) low levels of genetic subdivision between and within populations [4,66]; and 3) a mismatch distribution of pairwise differences among haplotypes indicating a sudden increase in expansion from a single population [67]. These signals of expansion from a single refugium have been reported in several North American bird species (e.g., multiple species reviewed in 68, downy woodpecker (Picoides pubescens) [20]) as well as many plant species, including several Pinus species that Clark's nutcracker specializes upon (reviewed in 13). Clark's nutcracker mitochondrial DNA exhibits a star-like phylogeny with one major shared haplotype, and short connections between this haplotype and others (Figure 2), in addition to a modeled paleodistribution that suggests that this species expanded from a single refugium after the LGM. This pattern is less complex than that found in whitebark pine postglacial expansion: colonization is thought to have occurred from three refugia (Oregon, Idaho, and Colorado) [25]. Given that Clark's nutcrackers range extends farther south than that of whitebark pine, nutcrackers may have been isolated in a more southern refugium and expanded into whitebark pine refugia from there, subsequently assisting with seed distribution northward. In addition to the aforementioned characteristics, temperate species that expanded north from a southern refugium have been shown to exhibit genetic differentiation between northern and southern populations corresponding to previously unglaciated and glaciated areas [12]. Southern populations were less impacted by glaciation and should show higher levels of nucleotide diversity and phylogeographic structure compared to northern populations [69]. This pattern occurs to varying degrees in Pinus species on which Clark's nutcracker commonly feed: whitebark pine (Pinus albicaulis) [25], white pine (P. monticola) [70], limber pine (P. flexilis) [71,72], and Ponderosa pine (P. ponderosa) [73,74]. However, this pattern does not distinctly appear in Clark's nutcracker. All southern and northern populations sampled exhibit high levels of  heterozygosity, nucleotide and allelic diversity, and limited differentiation from other populations, making it very challenging to pinpoint a single refugium for this species. In addition, while species distribution models for current suitable nutcracker habitat performed very well (high AUC values), paleodistribution models returned reduced probabilities (< 30%) of suitable LGM habitat throughout unglaciated areas. Taken together, these results point to a single refugium or multiple small, connected refugia, the location of which is not obvious through genetic signatures or paleodistribution modeling.

Conclusions and future research
Throughout its range, Clark's nutcracker exhibits low levels of population differentiation and does not appear to be limited by potential barriers to dispersal, though peripheral populations may be slightly isolated. Future research could integrate additional mitochondrial markers to elucidate the location of a glacial refugium or refugia for this species. Mitochondrial markers with lower rates of mutation in corvids (e.g. cytochrome B) [31] may retain historic genetic signals longer and thus provide additional historic information to this study. Samples at the northern and isolated southern extremes of the range may also increase our understanding of peripheral isolation and postglacial colonization for this species. Given the importance of Clark's nutcracker to dispersal of Pinus species in North America [21,22] and current conservation and management concerns for many Pinus species [75], further understanding of nutcracker dispersal may also assist in future forest management and recovery plans.