Population size estimates for plant species listed as threatened and endangered under the Endangered Species Act (ESA) are largely speculative. The lack of quantitative methods for assessing species abundance has contributed to ambiguity when assessing the status of ESA-listed species, including recognizing when species may be at an increased risk of extinction or determining when recovery has been achieved. In this paper we describe a sampling-based procedure used to estimate the minimum size of the greater Grand Valley population of the federally threatened Colorado hookless cactus (Sclerocactus glaucus). The estimation procedure applies plant density estimates derived from sampled macroplots to known habitat areas to obtain an estimate of the minimum population size for the entire area of occupation of the taxon. We found that previously reported population size estimates for the species were much lower than those resulting from our sampling-based approach.
INTRODUCTION
Understanding the abundance of species that are listed as threatened or endangered under the Endangered Species Act (“listed species”) is vital to determining their status and making management decisions aimed at their conservation (He and Gaston 2000). Despite a broad need for accurate assessments of the abundance of listed species, such measures are largely lacking. This information gap is especially severe for listed plant species, which are chronically overlooked and understudied as the result of insufficient resources allocated to endangered species programs (Schemske et al. 1994). In many cases, the distributions of listed plants are well defined and readily known while other, equally, if not more, important types of biological information, including species abundance and population trends, fall largely in the realm of speculation (Tear et al. 1995). This paradigm confounds the already difficult and uncertain decision-making processes associated with species conservation, potentially putting species at greater risk of extinction, or contributing to ambiguity in determining when a species has achieved recovery (Menges and Gordon 1996; Kaye et al. 2019).
Previous research has found that estimates of total population size pertaining to listed plant species are routinely based on the cumulative results of surveys or simple guesswork (Tear et al. 1993). While surveys are imperative to defining the spatial extent and distribution of the species of concern, there are at least two reasons why population size estimates derived from surveys may lead to imprecise assessments of abundance: (1) complete systematic inventories (i.e., censuses) require a prohibitive amount of time and resources to complete due to the spatial extent of all but the most geographically restricted populations; and (2) detection rates of rare plant species are often poor due to patchy distributions, low densities, and interannual variability in the detectable number of plants, depending on cryptic life phases, responses to environmental variables, or performance in previous years (Schemske et al. 1994; Le Lay et al. 2010).
Furthering the need for accurate assessments of the abundance of listed species, the US Fish and Wildlife Service (USFWS) has often considered benchmarks of minimum population size in guiding ESA-related policy decisions. Frequently minimum population sizes are included in species-specific recovery plans as a criterion to consider their delisting (Foin et al. 1998; Hoekstra et al. 2002; Zeigler et al. 2013). More recently, the USFWS has adopted the analytical Species Status Assessment framework that incorporates available biological information, including population sizes, to clarify the current status of species and inform both listing and delisting decisions. Therefore, the reliability of population size estimates improves decision making related to species protection.
In this paper, we present a sampling-based procedure whereby estimates of plant density obtained from sampled macroplots are applied to spatial occurrence data in order to arrive at a minimum population estimate of the listed plant species Colorado hookless cactus (Sclerocactus glaucus [K. Schumann]) L.D. Benson. In this case, we estimate the minimum population size in order to provide the most conservative estimate of the total population size, in part as a precaution because smaller populations are at a greater risk of extinction (Rabinowitz et al. 1986; Soulé 1987; Primack 2006).
Study Species
Colorado hookless cactus is a small, cylindrical barrel cactus endemic to the high-elevation deserts of western Colorado. The species typically grows as a single perennial stem from a taproot (BLM 2019). Most plants possess characteristic straight “hookless” abaxial spines, although variation in spine morphology is relatively common and has proven to be an unreliable trait for field identification (Schwabe et al. 2015). Individuals of reproductive stage will typically produce striking pink flowers in April and May. Across its limited range, the species exhibits a somewhat generalist distribution, inhabiting a variety of upland desert habitats and plant communities, ranging from alluvial river benches to salt desert shale barrens, gravelly colluvial washes, and sparse pinyon-juniper woodlands (USFWS 2010). Individual cacti are long-lived, capable of persisting for more than 10 y (BLM 2019). Due to its long-lived nature, population abundance remains relatively stable from year to year.
Colorado hookless cactus received Federal protection as a threatened species under the Endangered Species Act in 1979 (USFWS 1979). At the time of its listing, all Sclerocactus occurrences possessing straight central spines in western Colorado and northeastern Utah were classified as Sclerocactus glaucus. Phylogenetic analysis has since led to a series of taxonomic revisions ultimately dividing the Sclerocactus glaucus complex into three distinct species (Heil and Porter 2004; USFWS 2009). Two of these species (Sclerocactus brevispinus K.D. Heil and J.M. Porter and Sclerocactus wetlandicus Hochstätter) are endemic to the Uintah Basin in northeastern Utah while S. glaucus is restricted to western Colorado.
More recently, detailed molecular work has defined two geographically delineated populations of Colorado hookless cactus in western Colorado (McGlaughlin and Neale 2017). One of these populations is a relatively large and genetically diverse group associated with the Gunnison River and Grand Valley, hereafter referred to as the greater Grand Valley (GGV) population. The second of these populations is a smaller, genetically differentiated group associated with the Colorado River to the northeast. The first of these groups, the GGV population, is the focus of our study.
Taxonomic revisions have complicated the ability of resource managers to quantify range-wide population sizes and relative abundances for the three species leading to uncertainty in previous agency estimates. Database records housed at the Colorado Natural Heritage Program (CNHP), the agency responsible for tracking and ranking rare Colorado taxa, estimate the size of the GGV population of Colorado hookless cactus at 16,800 individuals (CNHP 2017). We performed a detailed review of these records and found that a number of the occurrences had not been observed in 20 or more years, or were ranked as either “historical,” “failed to find,” or “no data.” Additionally, the population values reported in individual records varied between partial and total counts of individuals within a given area and population estimates. In several instances population values were not representative of a complete mapped occurrence. Adding to uncertainty, not all inventory data compiled over the past 5 y has been integrated into the database.
While these records constitute the best available spatial representation of occupied Colorado hookless cactus habitat, long intervals between observations and the lack of a standardized methodology used to estimate site-level species abundance have contributed to uncertainty in the population estimates associated with the individual occurrence records, and therefore cast doubt on the cumulative population size. Based on these factors and observations in the field, we believe there are likely many more GGV Colorado hookless cactus individuals than accounted for in official reported estimates.
METHODS
Study Design
To estimate the minimum population size of GGV Colorado hookless cactus we implemented a two-stage sampling design (Thompson 2012) whereby 16 habitat areas were selected for sampling using a spatially balanced random sample. Within each of these habitat areas, a rectangular macroplot was subjectively sited to represent the larger habitat area. These macroplots functioned as the primary sampling units in the two-stage design. The macroplots were then sampled using rectangular quadrats (the secondary sampling units) selected using a simple random approach whereby each quadrat had an equal chance of selection. Mean plant density and population totals were then estimated for each of the 16 macroplots.
To be conservative in our estimation of the number of GGV Colorado hookless cactus plants per habitat area, we assumed that the only plants within each habitat area sampled were those we estimated within the macroplot itself. Under this assumption, the density of each habitat area could be calculated by dividing the estimated total number of plants within the macroplot by the area of the habitat, in square meters. Since the density calculations for each habitat area only accounted for the number of plants within each macroplot, and not the overall habitat area it was chosen to represent, we can be reasonably certain that the resulting estimate is the minimum number of plants per habitat area. We then used a ratio estimator to estimate habitat area density and applied this estimate to the total area of the study region to obtain the minimum population size of GGV Colorado hookless cactus. Sampling occurred between mid-September and the end of October in 2017.
Estimating the Area of Occupied Habitat
Understanding the spatial extent of the species of interest is foundational when calculating a population estimate using density, since this area calculation is the multiplier to which the resulting density estimate is ultimately applied. In our case, the spatial extent of GGV Colorado hookless cactus occupation was defined using inventory and mapping data that have been collected during the period since the species was listed in 1979. The resulting spatial dataset is the most comprehensive and complete representation of Colorado hookless cactus locations in Colorado (Figure 1).
To produce the area calculation, all spatial occurrence data from CNHP occurrence records were combined with BLM survey data from each of the field offices where the species occurs. The resulting spatial representation consisted of both point data and polygons. In situations where polygon data were concurrent between the two datasets, we used the union geoprocessing tool in ArcGIS Desktop (ESRI 2018) to combine the polygons to ensure that any given area was only counted once. Any additional overlapping areas were manually eliminated in ArcGIS with priority given to the dataset that was most current or complete.
In many cases, point data were clustered in areas of occupation, but because point data are nondimensional, we were not readily able to apply our plant density value to these areas. To obtain a conservative area calculation from the clustered points we converted these clusters to polygons using the point density tool in the spatial analyst extension for ArcGIS. The point density analysis used the merged point dataset to perform a density calculation: spatial data points that fell within a search area of 20 m2 were summed, then divided by a search radius of 50 m to derive a point density value for each 20 m2 cell (Silverman 1986). The point density output raster was then reclassified for density values >0.0005 points/m2 to define the clustered population areas (Figure 2). The resulting raster output was then converted to a polygon layer and unioned with the existing polygon dataset.
The resulting output from the GIS exercise was a constellation of 1744 individual polygons, each representing a discrete habitat area.
Figure 2.
Spatial sampling hierarchy of quadrat, macroplot, and habitat area at sampling Site 12. Point clusters encompassed by hatched polygons represent GGV Colorado hookless cactus point data that were converted to habitat area polygons—note data points that fell outside of our selection criteria and were excluded from our final area calculation.
Sample Macroplots
Sixteen habitat areas were ultimately selected for sampling based on a random, spatially balanced approach implemented using the Shiny spatially balanced sampling tool in the Landscape Toolbox (Karl 2007). Within each selected habitat area, a rectangular macroplot was subjectively defined to represent it. In most cases the macroplots were many times smaller than the habitat area they were chosen to represent. The macroplots ranged in size from 216 m2 to 2000 m2, with a mean of 749 m2 and a standard deviation of 448 m2 (Table 1).
Each macroplot was then sampled via a series of equally sized rectangular quadrats nested within the macroplot. Rectangular quadrats were chosen instead of square or circular quadrats because they give much more precise density estimates for plants that are clumped and/or irregularly spaced, a distribution pattern displayed by most plant species (Elzinga et al. 1998). Within each macroplot the necessary number of quadrats were selected for sampling using a simple random process. The total number of quadrats per macroplot (N) and the number of quadrats ultimately sampled per macroplot (n) varied by site (Table 1).
Each sampled quadrat was surveyed, and the total number of Colorado hookless cactus individuals was recorded. For each macroplot, density in terms of the mean number of plants/m2 was then calculated by dividing the estimated macroplot total by the area of the macroplot, expressed in m2. A power analysis was completed for each sample to ensure enough quadrats were sampled to ensure estimates of mean density and population size for each macroplot met a minimum criterion of being at least 80% confident of being within 30% of the estimated true value (Elzinga et al. 1998).
Table 1.
Summary statistics for each of the 16 samples.
Estimating Plant Density of the Habitat Areas
Using the estimated macroplot plant totals, density estimates were then produced for each of the 16 habitat areas. This was done by dividing the estimated macroplot total by the area of the given habitat area it represents, also expressed in m2. This procedure assumes that the only plants that occur within each habitat area were those estimated within its associated macroplot. This assumption is certainly not true and ultimately underestimates the true density of the habitat areas since they are many times larger than the macroplots and almost certainly contain an indeterminable number of Colorado hookless cactus plants not captured by the macroplot itself. In this way the density estimates for each of the habitat areas represent the minimum plant density for those areas.
Although this approach is conservative in the sense that it ignores any plants in the habitat area that are outside the macroplot (and there are certainly many such plants), it is not conservative in the sense that this method assumes the macroplot totals are precise. To address this issue, we used the estimated number of plants corresponding to the lower 90% confidence limit for each macroplot instead of the estimated macroplot total itself. Using the values corresponding to the lower confidence limit adds additional insurance that the estimated minimum population size of GGV Colorado hookless cactus is not biased high.
Ratio Estimation
It is tempting to calculate the mean density and variance for the collection of 16 habitat areas by averaging the 16 individual density values, although such an approach would only be warranted if the habitat areas were all the same size. Because of the difference in size of the 16 habitat areas, a ratio estimator is required to estimate the total population size and variance (Stehman and Salzer 2000). The ratio estimator employed is:
where 
is the estimated density, ȳ = sample mean number of plants per habitat area, and ā = sample mean area per habitat area sampled (Stehman and Salzer 2000; Thompson 2012:94).
The variance of this ratio estimator is:
Where 
number of plants in sampling unit (habitat area) u, au = area of sampling unit (habitat area) u, and n = number of habitat areas sampled (Stehman and Salzer 2000; Thompson 2012:94). Note that the variance formulas given in the two references cited above multiply the calculated variance by the finite population correction factor (fpc), (N – n)/N. The fpc is not included here because only a small fraction of the total habitat area was sampled and when n is very small compared to N, the fpc reduces to 1.
A confidence interval around the estimated density is constructed by first calculating the standard error of the estimate:
then applying it to the formula:
where t* is the t value from the t-distribution that corresponds to the desired confidence level and n – 1 degrees of freedom.
The estimated density and confidence interval derived from the sample of habitat areas, given in the number of plants/m2, is then converted into an estimate of minimum population size by multiplying this value by the total area of the habitat, also expressed in m2:
Where 
= estimated total number of plants,
= confidence interval around the total number of plants, and A= total area of occupied GGV Colorado hookless cactus habitat. Ratio estimation analysis and confidence interval generation was completed using the survey package (Lumley 2004 and 2019) in R 3.6.1 (R Core Team 2019).
RESULTS
Macroplot Densities
Estimated population totals and estimated population densities were calculated for each of the 16 macroplots sampled. Densities ranged from 0.063 plants/m2 to 0.623 plants/m2 based on an average study site of 749 m2. The mean density of the collection of 16 macroplots is 0.234 plants/m2 (90% CI [0.156, 0.312]). As previously discussed, this average density applies only to the total area inside the collection of macroplots and is not representative of the larger landscape where individual plants are often dispersed at low densities over large areas. Estimated macroplot totals ranged from 40 to 617 individual plants based on the same series of variously sized macroplots (Figure 3).
Table 2.
Output statistics resulting from the ratio estimation procedure using the lower 90% confidence level of the 16 macroplot totals.
Ratio Estimation
The estimation procedure calculates a ratio estimator of plant density using the value corresponding to the lower limit of the 90% confidence interval from each of the 16 estimated macroplot totals shown in Figure 3 and the habitat areas within which the estimates were obtained. This procedure assumes that the only plants within each mapped habitat area are those estimated within the macroplot, resulting in very conservative estimates of density for each habitat area and, ultimately, of the total GGV population size.
Determining a Minimum Population Estimate (
)
Based on our sample the estimated density is = 0.0061. Again, this estimate is based on the conservative assumption that the only plants that exist within each mapped habitat area we sampled were captured within its associated macroplot and corresponds to the lower 90% confidence limit of each of the 16 samples. The standard error of = 0:0012. was applied to the total area of habitat, A = 16,996,891 m2 in order to derive the estimated minimum population total = 103,086 plants ± 34,966 plants (90% CI [68,120; 138,053]; Table 2).
DISCUSSION
Our findings suggest that the size of the GGV population of Colorado hookless cactus is likely many times larger than previously reported. Based on our estimation procedure, approximately 68,000 plants should be considered the minimum size of the GGV population of Colorado hookless cactus. Despite the conservative approach taken to reach this minimum estimate, the result of our sampling-based estimation procedure is higher than previously reported estimates of its total population size. The CNHP estimates the size of the GGV population at 16,800 plants and estimates of the total range-wide number of Colorado hookless cactus (the GGV and the northern population) vary between 19,000 and 22,000 plants (CNHP 2017).
Throughout our process we have been vigilant to avoid biases that could result in overstating the total number of plants in the GGV population. One important measure taken to ensure a very conservative estimate was the application of the estimated number of plants corresponding to the lower 90% confidence limit for each macroplot in our calculations, rather than the estimated macroplot totals. If we were instead to apply the estimated macroplot totals the resulting estimate would range from approximately 115,000 to 231,000 plants. At the upper end, if we were to apply the values corresponding to the upper 90% confidence limit, the estimated minimum size would range from approximately 162,000 to 324,000 plants.
Figure 3.
Estimated macroplot totals for n = 16 Colorado hookless cactus macroplots. Error bars are 90% confidence intervals. The values corresponding with the lower limit of the confidence intervals were used in the ratio estimation procedure.
The primary factor inhibiting our ability to extrapolate our macroplot density estimates directly to the collection of habitat areas was the fact that the macroplots were subjectively sited. A completely random site selection process could have prevented the bias resulting from sampling targeted populations. Although such a completely random design would be preferable, it would be much more costly and time intensive to implement than the design employed here. Inevitably, a totally random selection process would require many more samples in order to overcome the variability between sampling units (many sampling units may fall in areas with few or no plants, while others may fall in areas with many plants).
While we are able to estimate the minimum size of the GGV population of Colorado hookless cactus, we are limited in our ability to approximate the true total population size, although we can conclude that the total population size is very likely higher than the minimum population size estimated using the procedure documented here. While a precise or exact estimate of population size would be ideal, understanding the minimum number of individuals in a species or population is still a useful metric and allows for conservation planning to occur under precautionary premises. This approach could be broadly applied to rare plant species whose range of occupation is well defined, especially long-lived perennial species that exhibit stable population trends, including both Sclerocactus brevispinus and S. wetlandicus and the northern population of Colorado hookless cactus.
Despite representing only a single case, this study highlights the potential disparity that can exist between subjective, ad hoc assessments of species abundance and those based on sampling. Resource managers should approach nonsampling estimates of population size with a measure of caution, realizing that such estimates may be highly inaccurate. In many cases, these inaccuracies may be biologically significant, as was the case here. While this study demonstrates that the total size of the GGV population of Colorado hookless cactus is larger than previously described, such a situation should be considered a best-case scenario. While one species' status might be less precarious than previously believed based on its population size, there are very likely other species and populations where the opposite reality is true.
Since abundance estimates remain a primary component of the relative status of a species that in turn may play an outsized role in the amount of resources allocated to its management—including whether or not it will receive federal protections—the threat posed by either underestimating or overestimating the abundance of a rare species may have major implications on its management and, ultimately, on its conservation outcome.
ACKNOWLEDGMENTS
The authors graciously acknowledge the contributions of David Sinton and Marian Sanone from the Bureau of Land Management Uncompahgre Field Office for their support in developing and fine-tuning the GIS process. This effort would not have happened without all those who assisted with field data collection, especially Sheila Cloud, Jenny Fausey, Robyn Oster, and Marian Sanone, in addition to Taryn Contento, Anna Lincoln, and Brooke Palmer. We also extend our gratitude to Rooster Barnhardt for his river navigation expertise and assistance in accessing remote reaches of the study area, and to three anonymous reviewers for their thoughtful comments and suggestions, which helped strengthen this manuscript.
