Enhancing georeferenced biodiversity inventories: automated information extraction from literature records reveal the gaps

Biogeography

Global biogeographical and macroecological studies require data on aggregate entities, such as location-specific occurrences of taxa and regional species assemblages, in order to understand emergent patterns at global and/or temporal scales (McGill, 2019). Assembly of such detailed yet broad-scale data is highly labor-intensive; the sampling effort required for a specific research question can be daunting for any one researcher or single research team. This is one reason why collaborative and often public databases have gained traction (Klein et al., 2019; Heberling et al., 2021). For instance, empirical global biogeographic analyses (Costello et al., 2017; Rivadeneira & Poore, 2020; Chaudhary et al., 2021; Hughes et al., 2021) are increasingly based on public databases of georeferenced taxonomic occurrences, such as the Ocean Biodiversity Information System (OBIS, www.obis.org) and the Global Biodiversity Information Facility (GBIF, www.gbif.org). Analyzing such georeferenced databases with tools that partially alleviate incomplete or biased sampling (Edgar et al., 2017; Kusumoto et al., 2020; Zizka et al., 2020; Grenié et al., 2022), allows us to address questions on large-scale distributions of clades, especially those that are well-represented in such databases. For less well-studied clades, however, prospects for obtaining large amounts of such data are lower. Answering pattern-based questions such as ‘how many species of clade z are found in location y’ and more process-oriented questions such as ‘how did the current latitudinal diversity gradient form’ both require location-specific taxonomic data in substantial volume. In addition, generalized biogeographic hypotheses have the potential to be supported more robustly if they include a greater diversity of clades.

Automated information retrieval

Automated information retrieval (Hirschberg & Manning, 2015) is one recent approach to complement the time-consuming manual activity of data compilation from the scientific literature. We seek to advocate for the enormous interdisciplinary potential of automated information retrieval, while acknowledging that it is a rapidly developing area at an early stage of development whose capabilities are just beginning to be realized. Automated text-mining is well-established in the biomedical realm (Percha, Garten & Altman, 2012; Christopoulou et al., 2020), but has only recently been adopted for biodiversity studies (Peters, Husson & Wilcots, 2017; Kopperud, Lidgard & Liow, 2019). As far as we are aware, automated text-mining has never been applied to the literature for extraction of taxon occurrences in given locations for the purpose of understanding biogeography (but see Page, 2019). Natural language processing is primarily concerned with understanding human language. This can involve a range of simple tasks, for example sentence delimitation and tagging part-of-speech (identifying nouns, verbs, adjectives etc.), to more complex tasks such as language translation, text generation, or text-to-speech. We use natural language processing tools (De Marneffe et al., 2014; Bojanowski et al., 2017) to facilitate a task known as relation extraction (Wang et al., 2022), which involves understanding relations among words in a sentence. We use these extracted relations to compile cheilostome text-mined occurrence data (TMO), which we compare with data from DB.

Cheilostomes

Cheilostome bryozoans, though less well-studied than several metazoan clades of similar size, are ubiquitous in benthic marine habitats. They are the most diverse order of Bryozoa with a conservatively estimated 4,921 extant described species (Bock & Gordon, 2013), or 83% of all Bryozoa. Bryozoans are ecologically important habitat builders (Wood et al., 2013) and are vital components of the marine food chain (Lidgard, 2008). Despite important analyses of regional species distributions (Clarke & Lidgard, 2000; López Gappa, 2000; Barnes & Griffiths, 2008; Hirose, 2017, Boonzaaier-Davids, Florence & Gibbons, 2020; Denisenko, 2020), their global species richness distribution has never been quantified. We argue that even with concerns about the incompleteness of public records for the purpose of inferring regional to global diversity patterns (e.g., Klein et al., 2019; Chollett & Robertson, 2020, Moudrý & Devillers, 2020; Hughes et al., 2021), it is worth exploring cheilostome data in such public databases. We do so in order to identify spatial gaps in sampling but also to ask if automated information retrieval can enhance the species occurrence data available in public databases, specifically OBIS and GBIF (henceforth shortened as DB).

Taxon occurrence

Taxon occurrence data from DB and TMO are not expected to be the same. DB records are comprised of diverse sources, but the bulk of contributions come from research institutes that are willing to digitize and organize field notes and collections (Saran et al., 2022). The bulk of TMO data are, on the other hand, based on published taxonomic literature. We ask if they could, separately or in combination, shed light on a long-standing biogeographic hypothesis in the bryozoological literature. Many different groups of organisms show the canonical LDG, a species richness peak in tropical regions and decreasing species richness towards the temperate and polar zones (Hillebrand, 2004; Menegotto, Kurtz & Lana, 2019). Despite being common across marine and terrestrial realms, and among diverse eukaryote clades, the LDG is not universal (Chaudhary et al., 2021). Marine extratropical bimodal species richness peaks have been observed, for example in brittle stars (Woolley et al., 2016), polychaetes (Pamungkas, Glasby & Costello, 2021), crustaceans (Rivadeneira & Poore, 2020), fishes (Lin et al., 2021), and brown macroalgae (Fragkopoulou et al., 2022), among other groups. Bimodality has also been suggested for cheilostome bryozoans (Schopf, 1970; Clarke & Lidgard, 2000; Barnes & Griffiths, 2008).

The TMO and DB data in combination support the view that the latitudinal diversity pattern of living cheilostomes is bimodal. These data reveal highest levels of estimated species richness in temperate latitudes, but TMO species richness has a peak in the northern hemisphere while DB has a peak in the temperate south. Moreover, two datasets differ significantly in the geographic richness patterns in Atlantic vs. Pacific ocean basins (Schopf, 1970; Barnes & Griffiths, 2008). We discuss the pros and cons of TMO and public databases such as OBIS and GBIF and how their differences can help us understand the uncertainties of the retrieved spatial diversity patterns, beyond what is estimated within the confines of each dataset.

DB data retrieval

We used the R-package robis (Provoost & Bosch, 2021) to access OBIS (2022), and the web interface of GBIF to retrieve latitude/longitude occurrence records of cheilostomes (GBIF dataset; DOI 10.15468/dl.58pd9h). Both databases were accessed on 21.02.2022. We removed records without species epithets. For taxonomic ambiguities such as cf., aff., we disregard the uncertainty; for instance, Microporella cf. ciliata becomes Microporella ciliata. Records with genus names that are not accepted according to either the Working List of Genera and Subgenera for the Treatise on Invertebrate Paleontology (DP Gordon, 2019, personal communication), World Register of Marine Species (WoRMS Editorial Board, 2022) or www.bryozoa.net (Bock, 2022) were also removed. For all unaccepted species names that are found in WoRMS, we translated the species name to the accepted species name according to WoRMS. We also dropped all Linnean binomials that were not found in WoRMS. The result is 831 unique genus names and 4,549 unique genus-species combinations (henceforth simply species) in 149,042 retained OBIS and GBIF records.

TMO (Text-Mined Occurrence) data retrieval

We follow a previously detailed text-mining procedure (Kopperud, Lidgard & Liow, 2019) with modifications. We extracted text from two collections of published works, our own corpus (3,233 pdf documents) and the GeoDeepDive archive (GDD, https://geodeepdive.org/), which contains full-text contents of journal articles. Only English language publications and those likely to feature extant bryozoans were used for information extraction (see Appendix S1). The natural language processing steps are detailed in the following sections, but we summarize them here. First, we perform a linguistic annotation task on the sentences of the text. Second, we extract so-called “candidates”, that is, pairs of location names and species names that co-occur in the same sentence. Third, we develop a location-name verifier to remove several mis-identified location names. Finally, we use a relational classifier that tells us whether the sentence actually said that the species was in the location.

Relation classifier

We manually annotated 4,938 unique candidates (species-location pairs) to form our training dataset. If the sentence explicitly stated or strongly implied that the taxon was found in the mentioned location, we labelled the candidate as ‘positive’. If not, we labelled the candidate as ‘negative’. These annotations were made by two persons, with intra- and inter-annotator accuracies of 91% (n = 200) and 85.6% (n = 211), respectively. This classifier is analogous to the location-name classifier (see previous section “Location-name verification”), except we used neither the indicator variables, nor the entire sentence. Instead, we used the sequence of tokens along shortest path in dependency grammar between the two spans in the sentence (see Xu et al., 2015; Kopperud, Lidgard & Liow, 2019). We trained our machine-classifier using the training set of these sequences to evaluate the classifier. This relation classifier achieved an F1 of 76.8%, accuracy of 73.1%, precision of 74.8%, recall of 78.9% and false positive rate of 34.3% (see Fig. S6b).

The estimated false positive rate of 34.3% is better than a random classifier baseline, but not as good as the false positive rates between human annotators at 14% and 16% (n = 200, assuming annotator A is correct, then evaluating annotator B, and vice versa). Similarly, the classifier accuracy of 73.1% is better than a weighted coin random classifier. If we use a random classifier that is as unbalanced as the labelled candidates (60% positives), this classifier gives an accuracy baseline of 52% (0.6² + 0.4² = 0.52). Yet, this is not as good as the intra- and inter-annotator accuracies at 91% and 85.6%, respectively. Incidentally, applying the relation classifier from Kopperud, Lidgard & Liow (2019) trained on relating taxon spans and geologic time interval spans, yielded virtually the same performance (evaluated on the same test data), despite not being trained on the specific task or having “seen” any location-name spans.

The relation classifier could potentially be improved by increasing the amount of training data, and/or using a more complex classifier, such as taking into account context-dependent word embeddings (Peters et al., 2018; Devlin et al., 2019). However, we believe that the major bottleneck is not lack of natural language understanding. Rather, the candidates themselves are not always linguistically sound, coherent and self-contained sentences. Specifically, much of the information that relays the relation between the two spans in the taxonomic literature of interest is coded in titles, sub-titles, variation in font type, font size, and spatial layout of the paragraphs. This introduces errors for the sentence splitting procedure in CoreNLP. The feature that we used to capture the relation (dependency grammar) is inherently limited since it is designed to work on relatively coherent, self-contained and complete sentences. Standard natural language processing tools are flexible and relatively easy to adopt, but for this particular problem it could be advantageous to use other, non-linguistic features to facilitate the information extraction, as has been suggested in the knowledge base creation literature (e.g., Schlichtkrull et al., 2018) but is outside the scope of our study.

We treat taxonomic ambiguities within TMO data in the same manner as OBIS and GBIF records (see section “DB data retrieval”).

From TMO location names to spatial data

Location names (e.g., New Zealand, Tasmania) are submitted to the Google geocoding service (https://developers.google.com/maps/documentation/geocoding/) to acquire a bounding box with four latitude-longitude coordinates and a centroid, based on Google’s defaults (Fig. S1). We remove species occurrences in locations represented by bounding boxes that are larger than about 2% of the Earth’s surface using area calculations assuming a spherical globe. The location names and their bounding boxes vary extensively in terms of their spatial resolution. The bounding boxes range from the scale of a few hundred meters across (for example “Fort Pierce Inlet”, in Florida, see Fig. S1), to major ocean basins like the “Mediterranean” or “Japan”. See Fig. S2 for how the bounding box sizes are distributed, and Fig. S3 for how the range-through diversity is affected when alternative cutoffs for the location size (lower or higher than ≈ 2%) are used.

Estimating latitudinal species richness

We initially evaluate species richness in thirty-six 5° latitudinal bands using two standard richness estimators that perform relatively well under a suite of conditions (Walther & Moore, 2005): Chao2 and Jackknife using the function specpool in the R package vegan (Oksanen et al., 2015). We treat these latitudinal bands as independent. We then repeat the procedure using thirty-six equal area bands, since areas represented within equal angle bands decrease poleward. To apply these estimators, we divide each (equal angle or area) latitudinal band into 5° longitudinal sampling units. We use the bias-corrected form of $Chao2=S_{obs}+Q_1^2(N-1)/(2N{Q}_2)$, and incidence-based Jackknife = S_obs + Q₁(N − 1)/N. Here, S_obs is the number of observed species in each band, N is the number of (longitudinal) sampling units, Q₁ is the number of species observed in only one sampling unit, and Q2 is the number observed in two sampling units. The quantities Q₁ and Q₂ are to an extent sensitive to our choice of dividing the latitudinal bands into 5° units. The larger the bands, the higher chance that there will be more singletons, since the discrete state space is merged or reduced (see Fig. S7). Note that we measure incidence as whether a species is either observed or not observed, for each geographical unit. Hence, any duplicate records in OBIS and GBIF do not inflate richness. Because terrestrial regions are not suitable habitats for marine cheilostomes, we mapped all landlocked longitudinal sampling bins (Fig. S4) based on a 1:10 m map of global coastlines (Patterson, 2019). We removed the landlocked bins prior to richness estimation. For DB data where spatial coordinates are points, it is trivial to assign data to sampling units. For TMO, we assume that a species occurs in all of the sampling units (latitudinal bands, and longitudinal bins) that intersect the bounding box associated with the location. TMO bounding boxes vary in size, but most are smaller in area than our sampling units (Fig. S2).

In addition to Chao2 and Jackknife estimators, we also determined range-through species richness. Here, we assume that a species spans its southernmost and northernmost occurrence record, regardless of whether it is observed in any intermediate latitudinal band. We did not split the incidence data by hemispheres. If one species was present in the southern and northern hemispheres, we would consider it present in all latitudinal bands inbetween. We acknowledge that all richness estimators, including the ones we chose, have different limitations (Gwinn et al., 2016). Confidence that the inferred patterns are real is improved by the extent that different estimators making different assumptions yield consistent results.

Capturing species richness: comparing DB and TMO

Applying the text-mining procedure to our corpora, we retrieved 1,408 species in 343 genera, and 1,400 unique location names among 7,204 TMO records. Only 23% of the species in the DB data that we retained were also in TMO. On the other hand, 68% of species in the TMO occurred in DB. 21% of the species richness is common to both (Fig. 1). In combination with DB data, we have species-location information from 4,910 species, almost tallying with the 4,921 described cheilostome species (Bock & Gordon, 2013).

Our machine-classifier achieved an accuracy of 73.1%, F1 of 76.8%, recall of 78.9%, FPR of 34.3% and precision of 74.8% as estimated with our test set (Fig. S6b). These results are substantially better than a random classifier baseline, but not as good as the human annotator repeatability. Specifically, the FPR among annotators is about 15% (n = 200). A random classifier that is as unbalanced as our training data (60% positive labels) would yield 60% false positives, but a random classifier equaling our classifier’s recall of 78.9% would have the same false positive rate of 78.9% (see Appendix S1).

Latitudinal species richness patterns

Combined TMO and DB data in plots of range-through species richness show a bimodal pattern with species richness peaks in both hemispheres surrounding 40° and −40° (Fig. 1). Inferred species richness in both of these peaks is about double that in the tropics (Fig. 2). The two data sources contribute different latitudinal constituents, as suggested by the limited overlap in their species composition (Fig. 1 inset).

Chao2 and Jackknife estimated species richness from DB shows two peaks between −20° and −45° that are more than double the next highest peak between 25° and 50° (Fig. 2A). In contrast, TMO estimated richness shows a highest peak between 30° and 45° (Fig. 2B). With minor exceptions in the Antarctic where spatial distortion is largest, equiangular and equi-areal bands yield nearly identical inferences (compare Figs. 2A, 2C). The latitudinal pattern appears smoother when using larger latitudinal band sizes (Fig. S7), but at the same time inflated, while retaining a qualitatively similar picture. Longitudinal sampling bins of varied sizes did not lead to notable variation for the Jackknife and Chao2 estimators (Fig. S8).

The northern hemisphere peak in richness (Fig. 1) and reflects TMO records from the Mediterranean and Japan, but also from the Atlantic Ocean (Figs. 3A, 3E), including the British Isles. Note that we did not include the Mediterranean as part of the Atlantic basin for Fig. 3. A portion of the TMO data are spatially imprecise, for example the location names “France”, “Spain” or “Morocco” may be associated with Mediterranean endemics, yet these records could contribute to the Atlantic richness counts in Fig. 3. The spatially precise DB data show a much lower peak in the Eastern Atlantic (Fig. 3E, orange line shifted slightly northward), reflecting data from the British Isles and northern Europe. Conversely, DB data mainly from Australia and New Zealand contribute disproportionately to the huge southern hemisphere peak. The richness captured by DB in Australia and New Zealand is not reflected by TMO species richness (Figs. 3B, 3D). The western Atlantic and eastern Pacific do not display such pronounced temperate zone peaks (Figs. 3C, 3F). Looking at individual ocean basins, TMO and DB are sometimes congruent and other times incongruent. For example, there is an absence of DB records in Japanese waters, and there are similarly few TMO and DB records in the Indian Ocean (Fig. 4).

Such varied regional species richness patterns are in part influenced by the geographic occurrence of samples. Figure 4 summarizes the relative distribution of species-location records for TMO and DB data as global heatmaps. For DB data, there are about one order of magnitude fewer records in tropical regions than for subtropical and temperate ones (Fig. S9a). While there are also fewer TMO records in tropical regions, the effect is not as pronounced (Fig. S9b). Northern and southern hemisphere species richness peaks in the two data sets (Figs. 3E, 3D) correspond with high regional densities of TMO and DB records, respectively (Figs. 3E, 3D).

Supporting information

• Appendix S1: Extended methods and supplementary figures.

• Appendix S2: Bibliographic references for TMO data.

The code and data required to reproduce the analyses and figures is available at www.zenodo.org/record/6770200 (DOI 10.5281/zenodo.6475658).