Basic diversity dynamics of scleractinian corals with divDyn

Author
Affiliation

Ádám T. Kocsis

GeoZentrum Nordbayern, Friedrich-Alexander-Universität Erlangen-Nürnberg, Loewenichstr. 28, D-91054, Erlangen, Germany

Goals and introduction

The purpose of this vignette is to guide new users through data acquisition and analysis, as well as to demonstrate the basic capabilities of the divDyn R (Kocsis et al. 2019) package. This tutorial will use Scleractininan corals as an example group which was selected for data quality reasons - and as we have readily available additional information about the coral taxa with which to customize the analyses.

Occurrence data download

To run analyses, first the occurrence data have to be downloaded and processed. You can download all PBDB occurrences using the chronosphere R package. The chronosphere is an initiative we have started at GeoZentrum Nordbayern in Erlangen, to facilitate the use, tracability, portability and version controlling of data in scientific analyses.

The package is available from the CRAN repositories and can be installed with the following line of code (you will probably have to select your preferred download mirror).

install.packages("chronosphere")

After the package is installed, it is ready for use, you can attach it with:

library(chronosphere)
Chronosphere - Evolving Earth System Variables
Important: never fetch data as a superuser / with admin. privileges!

Note that the package was split for efficient maintenance and development:
 - Plate tectonic calculations -> package 'rgplates'
 - Arrays of raster and vector spatials -> package 'via'
# omit 'ver="20240506"' for most up-to-date data!
dat <- fetch(src="pbdb", ver="20240506")

------------------------------------------------------------
Accessing chronosphere registry tables.
------------------------------------------------------------
Downloading registry tables.

------------------------------------------------------------
Item no. 2445, src:pbdb, ser: occs3, ver: 20240506, res: species.
------------------------------------------------------------

------------------------------------------------------------
Downloading import code.
------------------------------------------------------------

------------------------------------------------------------
Downloading data file.
------------------------------------------------------------
MD5 checksum passed.

If you use the data in publications, please cite its
reference(s), as well as that of the 'chronosphere' project.
Please remember to acknowledge the Paleobiology Database (http://paleobiodb.org) in your publication.
Alternatively, you can use: The Paleobiology Database Community. (year). Denormalized occurrence table of the Paleobiology Database. Zenodo. https://doi.org/10.5281/zenodo.7808071

This is a pre-processed version of the PBDB that was downloaded with this API call:

attributes(dat)$chronosphere$API
[1] "https://paleobiodb.org/data1.2/occs/list.csv?datainfo&rowcount&interval=Ediacaran,Holocene&show=class,classext,genus,subgenus,coll,coords,loc,paleoloc,strat,stratext,lith,env,ref,crmod,timebins,timecompare,refattr,entname,attr,ident,img,plant,abund,ecospace,taphonomy,etbasis,pres,prot,lithext,geo,methods,resgroup,refattr,ent"

You can copy and paste this URL to your browser and the PBDB would compile the current corresponding dataset for you. It normally takes about 15-30 minutes to get the dataset, another couple of minutes to download the .csv, which you still have to read in using read.table().

The PBDB API allows much more flexible download, but this denormalized version of the PBDB is enough for the overwhelming majority of the analyses. There are some cases, when the data you are interested in are not in this subset and you have to use the API yourself.

Data processing

Installing and loading divDyn

To process the stratigraphic information in the downloaded occurrence file, and do the analyses, you have to to install and attach the divDyn package.

The package is available from the CRAN repositories and can be installed with the following line of code (you will probably have to select your preferred download mirror).

install.packages("divDyn")

After the package is installed, it is ready for use, you can attach it with:

library(divDyn)
v0.8.2. See latest changes with 'news(package="divDyn")'. Default: 
Time flows forward with 'bin', and backward with 'age' and 'slice()'.

Initial filtering

The downloaded dataset needs to be subsetted to represent the order of scleractinian corals.

dat <- dat[dat$order=="Scleractinia",]

The analyses will be conducted at the level of genera.

The entries that have no valid genus names in the PaleoDB can be omitted at this stage. They will not provide us any information about the evolution of the group.

# omit non-genus entries
dat<- dat[!dat$genus=="", ]

The dataset should be cleaned as much as possible at this stage, for didactic purposes and for keeping things simple, this step is skipped in this tutorial.

Stratigraphic resolution

Every collection in the the PaleoDB has a user entered 'early_interval' and 'late_interval' field. These are fields are used to place the collection in stratigraphy/geological time. The 'early_interval' collection is mandatory for every entry while the 'late_interval' is only necessary if the taxon does not fit in the interval specified as 'early_interval'. Otherwise the 'late_interval' field is blank (empty quotes: "").

The divDyn package uses two frequently used, discrete time scales that we use as the finest stratigraphic resolution of our analyses: the stage level time scale (stages, 95 bins) and the ten million year time scale (tens, 49 bins). These can be attached with the data() function:

data(stages)
data(tens)

These will be useful for plotting as they contain the age estimates for the intervals, but they are not too useful for finding which interval the 'early_interval' and 'late_interval' fields point to. This can be found in the keys object.

data(keys)

This list is designed to be used with the categorize() function of the package, that combines the large number of entries and replaces them with the name of the group they belong to. This has to be run on both the 'early_interval' and the 'late_interval' fields.

# the 'stg' entries (lookup)
stgMin <- categorize(dat[ ,"early_interval"], keys$stgInt)
stgMax <- categorize(dat[ ,"late_interval"], keys$stgInt)

The output of this function is either an interval number that corresponds to the stages timescale, -1, indicating that the value was empty quotes, or NA, if the interval is either not in the lookup table or is too broad.

Because this scheme is designed for general use, the output of this function is a vector of character values. These have to be demoted to numeric values, otherwise their order would not be correct.

# convert to numeric
  stgMin <- as.numeric(stgMin)
  stgMax <- as.numeric(stgMax)

Now these two vectors have to be combined. There are more solutions to solve the problem of stratigraphic uncertainty, but for the sake of simplicity, let’s just omit every occurrence that cannot be assigned to a single stage in the stages timescale. For that, we will create an empty container, check whether an the interval entry for the occurrence satisfies the condition above, and if so, save the interval ID.

# empty container
dat$stg <- rep(NA, nrow(dat))

# select entries, where
stgCondition <- c(
# the early and late interval fields indicate the same stg
  which(stgMax==stgMin),
# or the late_intervar field is empty
  which(stgMax==-1))

# in these entries, use the stg indicated by the early_interval
  dat$stg[stgCondition] <- stgMin[stgCondition] 

Now every occurrence entry either has a single integer number in the stg column, or NA as its interval identifier. Note that the code above cannot be used for Cambrian and Ordovician occurrences due to stratigraphic problems. If you are interested in assigning these to stage-level intervals, we recommend to take a look at the vignette we wrote that describes global marine, Phanerozoic scale analyses:

https://github.com/divDyn/ddPhanero/blob/master/doc/1.0.1/dd_phanero.pdf

The assignments can be repeated for the 10-million year bins with the following commands:

# a. categorize interval names to bin numbers
  # categorize is the new function of the package
  binMin<-categorize(dat[,"early_interval"],keys$binInt)
  binMax<-categorize(dat[,"late_interval"],keys$binInt)
  # convert to numeric
  binMin<-as.numeric(binMin)
  binMax<-as.numeric(binMax)
# b. resolve max-min interval uncertainty
# final variable (empty)
  dat$bin <- rep(NA, nrow(dat))
# use entries, where
  binCondition <- c(
  # the early and late interval fields indicate the same bin
    which(binMax==binMin),
  # or the late_interval field is empty
    which(binMax==-1))
# in these entries, use the bin indicated by the early_interval
  dat$bin[binCondition] <- binMin[binCondition]

Sampling assessment

You can assess how well the stratigraphic resolution turned out by running the table() function on the interval-identifier vector:

table(dat$stg)

  48   49   54   55   56   57   58   59   60   61   62   63   64   65   66   67 
   1    5  114   72  448  910  602  101   86  248   80   87  493  383  231 1954 
  68   69   70   71   72   73   74   75   76   77   78   79   80   81   82   83 
1149  988  157  288  207  483 1023  390  536  128   71  158  298  742  396  289 
  84   85   86   87   88   89   90   91   92   93   94   95 
 296  295  466  652  519 1193 1772 1194 1465 1903 7829 1201 

Summing it tells us how many of the overall occurrences we can assign to stratigraphic stages.

sum(table(dat$stg))
[1] 31903
# which is a
sum(table(dat$stg))/nrow(dat)
[1] 0.8511779
# proportion of the total data.

As we cannot use unresolved occurrences in any way, coding will be easier if we omit them at this stage.

# omit unbinned
dats <- dat[!is.na(dat$stg),]

If you take a look at the stage object, you can notice that besides the few entries in stages 37, 49 and 50, Scleractinian become important fossils in the Anisian stage.

# omit Paleozoic
dats <- dats[dats$stg>52,]

You can get a more organized view of sampling parameters by running the binstat() function in divDyn, that calculates the occurrence, collection, and reference counts in a single line. This is the general use of the high-level function of the package: you state the occurrence data.frame as the first, and then the column names as additional arguments. The column tax stands for the taxon names, bin for the discrete time bins, coll for collection identifiers and ref for reference identifiers.

bsFull <- binstat(dats, tax="genus", bin="stg", 
  coll="collection_no", ref="reference_no")
The database contains duplicate occurrences (multiple species/genus).
bsFull$occs
 [1]   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
[16]   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
[31]   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
[46]   NA   NA   NA   NA   NA   NA   NA   NA  114   72  448  910  602  101   86
[61]  248   80   87  493  383  231 1954 1149  988  157  288  207  483 1023  390
[76]  536  128   71  158  298  742  396  289  296  295  466  652  519 1193 1772
[91] 1194 1465 1903 7829 1201

This output is organizes so that the index of the values in the vector match up the bin identifier (e.g. 60th value is for stg 60). The default setting of the function will output a message about duplicate occurrences. This warns us that there are collections with more than one genus entries in a collection (i.e. more than one species/genus). If you are interested in genus-level analyses, it is probably better to count these as one, which you can do with duplicates=FALSE option.

bs <- binstat(dats, tax="genus", bin="stg", 
  coll="collection_no", ref="reference_no", duplicates=FALSE)
bs$occs
 [1]   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
[16]   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
[31]   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA
[46]   NA   NA   NA   NA   NA   NA   NA   NA   99   66  357  769  444   77   72
[61]  228   77   74  413  292  206 1665  906  701  126  243  190  387  822  349
[76]  452  119   63  146  228  629  360  252  256  226  411  508  410  941 1192
[91]  964 1119 1474 6157  972

Plotting

Plotting these variables is probably better then just looking at the numbers. The package includes a powerful time-scale plotting function that lets you visualize the how time is broken down to discrete intervals. This highly-customizable function is built on the basic plot() function, and most of its arguments are inherited. The following function call will draw plot with in the past 250 million years, with series-level shading and system-level boxes at the bottom:

tsplot(stages, boxes="sys", boxes.col="systemCol", 
    shading="series", xlim=c(250, 0), ylim=c(0,2000))

For the sake of simplicity it is best to combine the output of binstat() wth that timescale object. Since the stg column indicates that same time slices in both objects, this step is easiest to accomplish with a simple merge() function call. This step should take care of any temporal mismatches between the two objects and will ensure that the correct values are plotted on th the geological timescale.

bsMerged <- merge(stages, bs, by="stg")

# inspect object
str(bsMerged)
'data.frame':   95 obs. of  25 variables:
 $ stg      : int  1 2 3 4 5 6 7 8 9 10 ...
 $ sys      : chr  "E" "E" "E" "Cm" ...
 $ system   : chr  "Ediacaran" "Ediacaran" "Ediacaran" "Cambrian" ...
 $ series   : chr  "Upper Ediacaran" "Upper Ediacaran" "Upper Ediacaran" "Terreneuvian" ...
 $ stage    : chr  "Avalon Assemblage" "White Sea Assemblage" "Nama Assemblage" "Fortunian" ...
 $ short    : chr  "Ava" "Whs" "Nam" "For" ...
 $ bottom   : num  580 560 550 539 529 ...
 $ mid      : num  570 555 544 534 525 ...
 $ top      : num  560 550 539 529 521 ...
 $ dur      : num  20 10 11.2 9.8 8 6.5 5.5 4.5 4 3.5 ...
 $ systemCol: chr  "#fddd88" "#fddd88" "#fddd88" "#94ac72" ...
 $ seriesCol: chr  "#fddd88" "#fddd88" "#fddd88" "#9db989" ...
 $ col      : chr  "#fcd589" "#fdd587" "#fed583" "#a9be93" ...
 $ occs     : num  NA NA NA NA NA NA NA NA NA NA ...
 $ colls    : num  NA NA NA NA NA NA NA NA NA NA ...
 $ refs     : num  NA NA NA NA NA NA NA NA NA NA ...
 $ SIBs     : num  NA NA NA NA NA NA NA NA NA NA ...
 $ occ1     : num  NA NA NA NA NA NA NA NA NA NA ...
 $ occ2     : num  NA NA NA NA NA NA NA NA NA NA ...
 $ ref1     : num  NA NA NA NA NA NA NA NA NA NA ...
 $ ref2     : num  NA NA NA NA NA NA NA NA NA NA ...
 $ u        : num  NA NA NA NA NA NA NA NA NA NA ...
 $ chao1occ : num  NA NA NA NA NA NA NA NA NA NA ...
 $ uPrime   : num  NA NA NA NA NA NA NA NA NA NA ...
 $ chao1ref : num  NA NA NA NA NA NA NA NA NA NA ...

Then you can draw the number of occurrences with lines().

tsplot(stages, boxes="sys", boxes.col="systemCol", 
  shading="series", xlim=c(250, 0), ylim=c(0,2000), ylab="Number occurrences")
lines(bsMerged$mid, bsMerged$occs)

We will mostly use the same basic function call for plotting, and typing all these arguments is time consuming. But we can spare some time, if we defined a wrapper function around the plotting that executes the plotting call the with the same arguments that we want to keep the same, while allowing us to change what we would like to change. This is fairly straightforward with the ... argument that allows you to pass arguments through function calls.

tp <- function(...) tsplot(stages, boxes="sys", boxes.col="systemCol", 
  shading="series", xlim=52:95, ...)

Now, every time you call the newly defined function tp() it will have the same colour, shading and x-interval, while you can change other arguments of the tsplot() function. For instance if you wanted see how the number of collections change over time, you can draw a similar plot as above, with:

tp(ylim=c(0,350), ylab="Number of collections") 
  lines(bsMerged$mid, bsMerged$colls)