Harvard Forest

Data Archive

HF296

Widespread Sampling Biases in Herbaria Revealed from Large-Scale Digitization 1656-2016

Related Publications

Data

• hf296-01: New England (preview)
• hf296-02: South Africa (preview)
• hf296-03: Australia (preview)
• hf296-04: rscripts

Overview

• Lead: Barnabas Daru, Charles Davis, Daniel Park, Aaron Ellison
• Investigators:
• Contact: Information Manager
• Start date: 1656
• End date: 2016
• Status: complete
• Location: New England, Australia, South Africa
• Latitude: -43.64806 to 47.45333 degrees
• Longitude: -73.72524 to 153.6242 degrees
• Elevation: 0 to 3254 meter
• Datum: WGS84
• Taxa:
• Release date: 2023
• Language: English
• EML file: knb-lter-hfr.296.4
• DOI: digital object identifier
• EDI: data package
• DataONE: data package
• Related links:
• Study type: historical
• Research topic: biodiversity studies; historical and retrospective studies; regional studies
• LTER core area: population studies
• Keywords: spatial variability, species abundance, species composition, species diversity
• Abstract:

  Non-random collecting practices may bias conclusions drawn from analyses of herbarium records. Recent efforts to fully digitize and mobilize regional floras offer a timely opportunity to assess commonalities and differences in herbarium sampling biases.

  We determined spatial, temporal, trait, phylogenetic, and collector biases in ~5 million herbarium records, representing three of the most complete digitized floras of the world: Australia (AU), South Africa (SA), and New England, USA (NE)

  We identified numerous shared and unique biases among these regions. Shared biases included specimens i) collected close to roads and herbaria; ii) collected more frequently during spring; iii) of threatened species collected less frequently; and iv) of close relatives collected in similar numbers. Regional differences included i) over-representation of graminoids in SA and AU and of annuals in AU; and ii) peak collection during the 1910s in NE, 1980s in SA, and 1990s in AU. Finally, in all regions, a disproportionately large percentage of specimens were collected by a few individuals. These mega-collectors, and their associated preferences and idiosyncrasies, may have shaped patterns of collection bias via ‘founder effects’.

  Studies using herbarium collections should account for sampling biases and future collecting efforts should avoid compounding these biases.

• Methods:

  Spatial biases

  First, we evaluated the density of sampling localities across the focal regions using Delaunay triangulation polygons, which measure the land area covered by each sampling locale (Fortune, 1992). Larger triangles indicate sparser collecting effort, whereas smaller triangles indicate more concentrated effort. Second, we examined infrastructure bias by calculating the minimum distance of each collection locality to the nearest major road (GADM, 2015) and herbarium (following Thiers, 2016). Our dataset of roads derives from the publicly available Digital Chart of the World, which was compiled by the US Defense Mapping Agency from 1:1,000,000 scale paper maps (ESRI, 1992). All roads appearing at this scale were included in our analyses. Although this dataset includes only larger roads and has not been updated since 1992, it likely represents the most comprehensive digital record of roads around the world. We then compared these distances to those generated by a null model (1000 iterations) in which the same number of sample points was randomly (Poisson) distributed across each geographic region. Third, we mapped geographic biases in sampling density, defined as areas of excessive (hotspots) or insufficient (coldspots) collection (Hijmans et al., 2000). Hotspots and coldspots were determined at a spatial grain of 0.25° × 0.25° based on the number of specimens per grid cell, and identified using the 2.5% threshold (Ceballos and Ehrlich, 2006; Orme et al., 2005; Daru et al., 2015), based, respectively, on the 97.5th and 2.5th percentile values in the number of specimens collected per grid cell. Spatial distance calculations were done with the functions dist2Line and spDists in the R packages sp (Bivand et al., 2013) and geosphere (Hijmans, 2015), respectively. In our final predictive model of sampling density, we also included human population density (CIESIN, 2016), sampling localities, infrastructure (distance to herbaria and roads), number of specimens collected, and elevation or topographic relief.

  Temporal bias

  For each regional flora, we explored bias at several temporal scales. Collection dates ranged from 20 May 1664 to 9 January 2016 (AU), 15 November 1656 to 6 June 2016 (SA), and 28 July 1687 to 4 May 2016 (NE). We hypothesized that collectors tended to avoid fieldwork during unfavorable conditions (e.g., weekdays, winter, wartime). To test for temporal bias, we first re-coded collection dates as days of the week (Sunday = 1, Monday = 2, etc.), and day of year (DOY; where January 1 = 1 DOY and December 31 = 365 DOY, etc.). We then used a Rayleigh test of directional statistics in the R package circular (Agostinelli and Lund, 2013) to test whether each of these collection dates were randomly distributed against all dates spanning the entire duration of plant collection. If P less than α = 0.05, we rejected the null hypothesis of temporal uniformity at scales of weeks, days of the year, or decades.

  Trait bias

  We used customized R scripts to harvest information on growth duration (annual vs. perennial), growth form (woody vs. herbaceous), and height for each species from online regional databases (all accessed in June 2016), including: New South Wales Flora Online (http://plantnet.rbgsyd.nsw.gov.au); JSTOR Global Plants (https://plants.jstor.org); Atlas of Living Australia (http://bie.ala.org.au); Plants of Southwestern Australia (http://keys.lucidcentral.org); the African Plant Database (http://www.ville-ge.ch); Plants of Southern Africa (http://www.plantzafrica.com); Plant Resources of Tropical Africa (http://www.prota4u.org); Flora of North America (http://www.efloras.org); and the USDA Plants Database (http://plants.usda.gov). We then manually checked these data for inconsistencies in terminologies for defining certain traits. For example, ‘vines’ vs. ‘lianas’ for climbers, ‘forbs’ vs. ‘herbs’ for herbaceous life forms, ‘biennial’ for perennial growth duration. Extinction risk assessments for each species were retrieved from the IUCN Red List database (www.iucnredlist.org, accessed August 2016), which uses the following categories: Data Deficient (DD), Least Concern (LC), Lower Risk/Conservation Dependent (LR/CD), Near Threatened (NT), Vulnerable (VU), Endangered (EN), Critically Endangered (CR), and Extinct (EX). We grouped these narrow categories into two broader threat categories, threatened (EX+CR+EN+VU) or not threatened (LR/CD+NT+LC), following Yessoufou et al. (2012).

  Trait bias was evaluated using a Chi-squared test to contrast the number of observed specimens collected per species with the abundance of a species if specimen collection was equal across all species for each trait category. Because of dramatically unequal sampling effort in some species – e.g., Senna artemisioides with 10,167 specimens vs. Eucalyptus cordieri with only one – and the low coverage of taxa with available trait data, we randomly sampled 50 specimens from each available species with trait data using 1000 randomizations. Species with less than 50 specimens were excluded from this analysis.

  Phylogenetic bias

  We assessed phylogenetic signal in collection frequency as a measure of phylogenetic bias using two different tests (Wolkovich et al., 2013). A strong phylogenetic signal – closely related species sharing similar collection frequency – would suggest phylogenetic bias in collections. We first assembled a phylogeny using Phylomatic (Webb and Donoghue, 2005), enforcing a topology that assumed the APG III backbone (tree R20120829). This phylogeny included all species in our analysis, but provided only an approximate degree of relatedness based on taxonomic hierarchy at family level; many relationships, especially within genera, were unresolved. This is problematic because recent theoretical and empirical studies have shown that a lack of resolution in a community phylogeny may mask significant patterns by reducing statistical power (Schaefer et al., 2011) or suggest significant phylogenetic patterns that are not supported by more completely resolved phylogenies (Davies et al., 2012).

  To alleviate these concerns, we also tested for phylogenetic bias by including only those species sampled in the dated molecular phylogeny inferred from seven genes for 32,223 plant species (Zanne et al., 2014). Although this phylogeny has been criticized (Edwards et al., 2015), it nonetheless represents the single largest phylogeny to date for flowering plants. The taxon sampling for testing phylogenetic bias included 5814 species from AU, 3568 from SA, and 4269 from NE.

  We estimated phylogenetic signal using three common metrics: Abouheif’s Cmean (Abouheif, 1999), Blomberg’s K (Blomberg et al., 2003), and Pagel’s lambda (λ) (Pagel, 1999). Significance was assessed by comparing observed values to a null distribution created by shuffling the trait values across the tips of the phylogeny 1000 times. Pagel's λ uses a maximum-likelihood method with branch-length transformation to estimate the best-fit of a trait against a Brownian model. Values of Pagel’s λ range from 0 (no phylogenetic signal) to 1 (strong phylogenetic signal). Both Blomberg’s K (a significant phylogenetic signal is indicated by a K value > 1) and Pagel’s λ were calculated using the R package phytools (Revell, 2012). Abouheif’s Cmean was calculated using adephylo (Jombart and Dray, 2008). We tested the sensitivity of our analysis by exploring phylogenetic signal in collecting effort across nine well-sampled NE clades: Asteraceae, Brassicaceae, Cyperaceae, Ericaceae, Fabaceae, Lamiaceae, Poaceae, Ranunculaceae, and Rosaceae.

  In addition to phylogenetic signal, we also used phylogenetic generalized least squares regressions (PGLS) in the R package caper (Orme et al., 2012) to model collecting effort per species in each region as a function of species evolutionary ages, evolutionary distinctiveness (ED), and ‘evolutionary distinctiveness and global endangerment’ (EDGE; Isaac et al., 2007). Species ages were measured as the length of terminal branches (BL) linking species on a phylogenetic tree. ED measures the degree of phylogenetic isolation of a species, whereas the EDGE metric was determined by calculating the ED score of each species (Isaac et al., 2007) and combining it with global endangerment (GE) from IUCN conservation categories: EDGE = ln(1 + ED) + GE × ln(2), where GE represents expected probability of species extinction over a 100-year period (Redding and Mooers, 2006) categorized as follows: least concern = 0.001, near Threatened and Conservation Dependent = 0.01, Vulnerable = 0.1, Endangered = 0.67, and Critically Endangered = 0.999.

  Last, we examined the phylogenetic structure of collecting efforts across decades to test for patterns of phylogenetic overdispersion and clustering through time. Temporal phylogenetic structure by decade (i.e., 1901-1910, 1911-1920, etc.) was evaluated using the net relatedness index (NRI) and nearest taxon index (NTI; Webb et al., 2002, 2008). NRI describes a tree-wide pattern of phylogenetic dispersion, whereas NTI evaluates phylogenetic structure towards the tips of the phylogeny. Negative values of NRI or NTI indicate phylogenetic overdispersion whereas positive values indicate phylogenetic clustering.

  Collector bias

  We determined collector bias by tabulating the number of specimens amassed by each collector in the three floras. We then examined Pearson’s product-moment correlation between the numbers of specimens collected per collector with the number of species collected per collector.

  Known data issues

  hf296-01-new-england.csv

  The variable “year” contains the year in which a record was collected. Some records may have zero instead of NA to indicate missing information. In analysis, users may need to combine information from "year", "month", and "day" fields to generate accurate date information.

  hf296-02-south-africa.csv

  The variable "full.number" contains unique identifiers for each record. Values may include a mix of letters and numbers or may be missing. Users may need to generate their own unique identifiers.

  hf296-03-australia.csv

  The variables "occurrence.date” and “year” contain the date or year, respectively, on which a record was collected. Some records may have zero instead of NA to indicate missing information. In analysis, users may need to combine information from "year", "month", and "day" fields to generate accurate date information.

• Organization: Harvard Forest. 324 North Main Street, Petersham, MA 01366, USA. Phone (978) 724-3302. Fax (978) 724-3595.

• Project: The Harvard Forest Long-Term Ecological Research (LTER) program examines ecological dynamics in the New England region resulting from natural disturbances, environmental change, and human impacts. (ROR).

• Funding: National Science Foundation LTER grants: DEB-8811764, DEB-9411975, DEB-0080592, DEB-0620443, DEB-1237491, DEB-1832210.

• Use: This dataset is released to the public under Creative Commons CC0 1.0 (No Rights Reserved). Please keep the dataset creators informed of any plans to use the dataset. Consultation with the original investigators is strongly encouraged. Publications and data products that make use of the dataset should include proper acknowledgement.

• License: Creative Commons Zero v1.0 Universal (CC0-1.0)

• Citation: Daru B, Davis C, Park D, Ellison A. 2023. Widespread Sampling Biases in Herbaria Revealed from Large-Scale Digitization 1656-2016. Harvard Forest Data Archive: HF296 (v.4). Environmental Data Initiative: https://doi.org/10.6073/pasta/3126b2ce0916181327e3a9d0d8804e50.

Detailed Metadata

hf296-01: New England

1. family: family taxonomic information
2. genus: genus taxonomic information
3. specific.epithet: second element in the Latin binomial (scientific) name of a species, following the genus name, used to distinguish a species from others within the same genus
4. decimal.longitude: geographic decimal longitude where a record was collected (unit: degree / missing value: NA)
5. decimal.latitude: geographic decimal latitude where a record was collected (unit: degree / missing value: NA)
6. year: the year in which the specimen was collected; missing years noted by NA or 0
7. month: month of the year in which a record was recorded, with 1=January, 2=February, and 12=December
8. day: day of the month on which a record was collected
9. species: the Latin binomial (scientific) name of the specimen, comprising the generic element followed by the second epithet
10. country: name of the country in which the occurrence was recorded
11. state.province: description of the next smaller administrative geographic division than a country (state, province, canton, department, region, etc.) in which the specimen was recorded
12. county: description of the next smaller administrative geographic division than a state or province in which the specimen was recorded
13. municipality: description of the next smallest administrative geographic division than a county in which the specimen occurs
14. locality: specific description of the place where a specimen was recorded

hf296-02: South Africa

1. principal.collector: main person who collected the records
2. additional.collectors: additional person(s) who collected the records
3. full.number: unique identifiers as serial numbers for each record
4. family: taxonomic classification of the record at the family level
5. species: Latin binomial (scientific) name of the specimen, comprising the generic element followed by the specific epithet
6. taxon: full Latin binomial name of the specimen (genus + species), followed by name of the author(s) who validly published the scientific name
7. country.name: name of the country in which the occurrence was recorded
8. major.region: description of the next smaller administrative geographic division than a country (province) in which the specimen was recorded
9. locality.text: specific description of the place where a specimen was recorded
10. qds: the quarter degree square where a record was collected. For example, a QDS of 3419AD = longitude 19.375 and latitude −34.375.
11. basis.of.record: basis upon which the record was derived
    • PreservedSpecimen: occurrence record describing a preserved specimen usually in a museum or herbarium
    • FossilSpecimen: occurrence record describing a fossilized specimen
    • LivingSpecimen: occurrence record describing a living specimen, e.g. cultivated plants in a garden
    • Observation: occurrence record describing an observation
    • Literature: occurrence record based on literature alone
12. collection.day: day of the month on which a record was collected
13. collection.month: month of the year in which a record was collected, with 1 = January, 2 = February, and 12 = December
14. collection.year: the 4 digit year in which the specimen was collected
15. latitude: geographic decimal latitude where a record was collected (unit: degree / missing value: NA)
16. longitude: geographic decimal longitude where a record was collected (unit: degree / missing value: NA)
17. iucn: conservation status of each record according to International Union for the Conservation of Nature
    • EX: extinct
    • EW: extinct in the wild
    • RE: regionally extinct
    • CR PE: critically endangered, possibly extinct
    • CR: critically endangered
    • EN: endangered
    • VU: vulnerable
    • NT: near threatened
    • LC: least concern
    • DDD: data deficient - insufficient information
    • DDT: data deficient - taxonomically problematic
    • NE: not evaluated

hf296-03: Australia

1. family: family taxonomic information
2. genus: genus taxonomic information
3. species: species information
4. latitude: geographic decimal latitude where a record was collected (unit: degree / missing value: NA)
5. longitude: geographic decimal longitude where a record was collected (unit: degree / missing value: NA)
6. state: subregion/state/province where the record was collected
7. occurrence.date: date on which the specimen was collected, in the format, e.g. "Sun Apr 10 10:00:00 EST 2016
8. year: the year in which the specimen was collected; missing years noted by NA or 0
9. month: month in which a record was collected
10. day: day of the month on which a record was collected
11. julian: day of year on which a record was collected, where January 1=1 (unit: nominalDay / missing value: NA)
12. basis.of.record: basis upon which the record was derived
    • PreservedSpecimen: occurrence record describing a preserved specimen usually in a museum or herbarium
    • FossilSpecimen: occurrence record describing a fossilized specimen
    • LivingSpecimen: occurrence record describing a living specimen, e.g. cultivated plants in a garden
    • Observation: occurrence record describing an observation
    • Literature: occurrence record based on literature alone

hf296-04: rscripts

• Compression: zip
• Format: R script
• Type: script