# Evolutionarily Distinct and Globally Endangered (EDGE) bird species

The idea behind this post was born after a conversation with Simon Valdez-Juarez, a fellow graduate student at SFU. He mentioned an idea he had for creating a Shiny web app that allows people to see the evolutionary distinctiveness of bird species in their eBird checklists. Given that I have played around with wrangling eBird checklists to display summary stats of sightings by location and year (the myebird package, see a previous blog post about it here), we decided to join forces to create this website. The functions introduced below will be the backbone of such a project.

## What is ED?

Evolutionary Distinctiveness (ED) is a measure of species’ relative contributions to phylogenetic diversity. It is calculated by dividing the total phylogenetic diversity of a clade among all its members. David Redding, a former PhD student at SFU, came up with the idea of ED for his M.Sc. thesis at the University of East Anglia. Simply put, the ED score of a species is the sum of all the branch lengths to the root of the tree, each divided by the number of species each branch section is shared among. I’ll illustrate this using four hypothetical species (A, B, C, and D). For each branch, the number above represents the length of the branch (in million years), while the number below the branch shows the number of species that share that branch, so ED is the sum of all these fractions from the root of the tree to each tip (Figure is adapted from Isaac et al. [2007]):

ED is a pretty straightforward metric that can help identify individual (perhaps flagship?) species of high evolutionary distinctiveness.This metric can be combined with a threat score based on the Red List status compiled by the International Union for the Conservation of Nature (IUCN). The Red List status of each species can be converted to a quantitative measure by converting qualitative status into ordered integers, in this case denoted as GE: thus a species assess as “Least Concern” would have a GE score of 0, a “Near Threatened” species would be scored as a 1, and so on up to a species assessed as “Critically Endangered”, which would have a score of 4. Thus, by using Evolutionary Distinctiveness (ED) and Global IUCN Red List status (GE) you can create a composite score of Evolutionary Distinctive and Globally Endangered (EDGE) species. The one proposed by Isaac et al. (2007) and subsequently used by Jetz et al. (2014) on birds is as follows:

$EDGE = ln(ED + 1) + GE * ln(2)$

The data from Jetz et al. (2014) include a bird phylogeny covering almost 10,000 species (available here) and their associated EDGE scores and ranks (available here); these data are the basis for these R functions and I would not have been able to create them without their work.

## The myedge and mypd functions

To being playing with bird trees and your eBird data, you first have to install the myebird package from GitHub using devtools::install_github and load it in R:

devtools::install_github("sebpardo/myebird")
library(dplyr) # for piping and data wrangling
library(myebird)
library(ggplot2)

There are three functions in this package that I’ve already written about here (ebirdclean, myebird, and myebirdcumul). I have added two new functions: myedge and mypd, as well as a dataset named edge.

To use these functions, we first need to download all our data from eBird: an email with a download link will be sent to you after pressing “Sumbit” in the following eBird link (make sure you’re logged in to eBird!): http://ebird.org/ebird/downloadMyData. We then clean this csv file using the ebirdclean function:

mybirds <- ebirdclean("MyEBirdData.csv")

To obtain ED scores as well as EDGE scores and ranks of the species included in your eBird checklist, you can simply use the myedge function, whose first argument is the output provided by ebirdclean, and the second argument is edge.cutoff, which is the cutoff for species to be shown in terms of ranking (defaults to returning all species):

myedgebirds <-myedge(mybirds, edge.cutoff = 500)
myedgebirds %>% select(-sciName)
## # A tibble: 12 × 5
##                     comName  ED.Score EDGE.Score EDGE.Rank
##                       <chr>     <dbl>      <dbl>     <int>
## 1    Peruvian Diving-Petrel 28.418480   5.461065        29
## 2  Australian Painted-Snipe 26.901620   5.408126        32
## 3             Andean Condor 43.702267   4.493171       228
## 4          Marbled Murrelet  8.815285   4.363382       277
## 5       Hutton's Shearwater  8.341314   4.313888       293
## 6           Andean Flamingo 15.887414   4.212863       339
## 7              Magpie Goose 62.675689   4.153803       372
## 8                    Osprey 56.063247   4.044160       431
## 9                Great Knot 12.885775   4.017159       451
## 10         Beach Thick-knee 26.631221   4.012094       459
## 11   Yellow-nosed Albatross  5.730249   3.986054       471
## 12        Providence Petrel 11.894030   3.943059       500
## # ... with 1 more variables: sciName.edge <chr>

We can also easily subset to see the ED and EDGE ranks and scores for a single country or even region:

mybirds %>%
filter(Country == "Chile") %>%
myedge(edge.cutoff = 1000) %>%
select(-sciName)
## # A tibble: 15 × 5
##                         comName  ED.Score EDGE.Score EDGE.Rank
##                           <chr>     <dbl>      <dbl>     <int>
## 1        Peruvian Diving-Petrel 28.418480   5.461065        29
## 2                 Andean Condor 43.702267   4.493171       228
## 3               Andean Flamingo 15.887414   4.212863       339
## 4        Black-browed Albatross  5.340944   3.926469       509
## 5     Diademed Sandpiper-Plover 23.787857   3.903501       525
## 6              Peruvian Pelican 22.821709   3.863744       553
## 7               Westland Petrel 10.252474   3.806882       599
## 8              Chilean Flamingo 20.040577   3.739600       657
## 9          White-chinned Petrel  8.965589   3.685432       704
## 10             Humboldt Penguin  8.335894   3.620161       754
## 11            Pied-billed Grebe 35.996441   3.610822       763
## 12 South American Painted-Snipe 34.544290   3.570780       805
## 13                Black Vulture 33.390317   3.537775       858
## 14             James's Flamingo 15.887414   3.519716       875
## 15              Royal Albatross  7.167664   3.486477       931
## # ... with 1 more variables: sciName.edge <chr>
mybirds %>%
filter(County == "Metro Vancouver") %>%
myedge(edge.cutoff = 1000) %>%
select(-sciName)
## # A tibble: 6 × 5
##                  comName  ED.Score EDGE.Score EDGE.Rank
##                    <chr>     <dbl>      <dbl>     <int>
## 1       Marbled Murrelet  8.815285   4.363382       277
## 2                 Osprey 56.063247   4.044160       431
## 3      Red-throated Loon 49.484481   3.921666       511
## 4      Pied-billed Grebe 35.996441   3.610822       763
## 5       Long-tailed Duck  8.053587   3.589455       786
## 6 Golden-crowned Kinglet 32.860078   3.522237       873
## # ... with 1 more variables: sciName.edge <chr>

myedge also can provide a separate data frame for each country when used in conjunction with dplyr::do:

mybirds %>%
group_by(Country) %>%
do(myedge(., edge.cutoff = 500)) %>%
select(-sciName)
## Source: local data frame [17 x 6]
## Groups: Country [6]
##
##                     Country                  comName  ED.Score EDGE.Score
##                       <chr>                    <chr>     <dbl>      <dbl>
## 1                 Australia Australian Painted-Snipe 26.901620   5.408126
## 2                 Australia      Hutton's Shearwater  8.341314   4.313888
## 3                 Australia             Magpie Goose 62.675689   4.153803
## 4                 Australia                   Osprey 56.063247   4.044160
## 5                 Australia               Great Knot 12.885775   4.017159
## 6                 Australia         Beach Thick-knee 26.631221   4.012094
## 7                 Australia   Yellow-nosed Albatross  5.730249   3.986054
## 8                 Australia        Providence Petrel 11.894030   3.943059
## 9                    Canada         Marbled Murrelet  8.815285   4.363382
## 10                   Canada                   Osprey 56.063247   4.044160
## 11                    Chile   Peruvian Diving-Petrel 28.418480   5.461065
## 12                    Chile            Andean Condor 43.702267   4.493171
## 13                    Chile          Andean Flamingo 15.887414   4.212863
## 14               Costa Rica                   Osprey 56.063247   4.044160
## 15                Nicaragua                   Osprey 56.063247   4.044160
## 16 United States of America         Marbled Murrelet  8.815285   4.363382
## 17 United States of America                   Osprey 56.063247   4.044160
## # ... with 2 more variables: EDGE.Rank <int>, sciName.edge <chr>

## Problems with ED

As it has been previously pointed out, there are some issues with composite indices of evolutionary distinctiveness. One important problem with ED estimates is that they can be very sensitive to splits or lumps at the species level. Let’s take the hypothetical example presented in the tree above, but this time using a taxonomy that considers species A to actually be three separate species (A, E, and F) which split around 0.2 million years (MY) before present:

Using this taxonomy, the ED score of species A drastically drops, even though it is still nested in a very long branch. One way around this is to change the focus from scoring individual species (which is sexy) and instead look at the phylogenetic diversity of a community of species (not as sexy).

## Phylogenetic Diversity

One neat way to estimate phylogenetic diversity of a clade is to simply add all branches of a tree; this represents the total evolutionary time covered by the species in the clade/tree. Again with our hypothetical four species, if we wanted to calculate PD for species B and C:

Our estimate of PD for species B and C is 4 MY of evolution. Now, if we also include species D, our estimate of PD will increase only slightly as species C and D are closely related.

On the other hand, if we include species A instead of species D, the phylogenetic diversity of this group will be much higher, as including species A means a much longer branch of evolutionary time is now being added:

The function mypd does exactly this: it calculates the phylogenetic distance of the species in your eBird checklist based on the Jetz et al. phylogeny. It is very straightforward to run, however it averages values across 20 very large trees, and therefore is a bit slow to run (shouldn’t take more than a minute though):

mypdvalues <- mypd(mybirds)
mypdvalues
## # A tibble: 1 × 3
##    mean_pd median_pd    sd_pd
##      <dbl>     <dbl>    <dbl>
## 1 20917.93  20755.27 1247.857

This means that the combined phylogenetic distance of the species recorded in all my eBird checklists encompasses 20.8 billion years of evolution!

As with the myedge function, you can look at PD values by country:

mybirds %>%
group_by(Country) %>%
do(mypd(.))
## Source: local data frame [16 x 4]
## Groups: Country [16]
##
##                             Country   mean_pd median_pd     sd_pd
##                               <chr>     <dbl>     <dbl>     <dbl>
## 1                         Argentina 2426.8959 2403.6887 121.87432
## 2                         Australia 7679.9299 7596.6293 466.10298
## 3                            Canada 6226.8165 6172.2301 357.44688
## 4                             Chile 5746.2029 5813.5036 277.16189
## 5                             China 1615.5454 1605.0135 104.20097
## 6                        Costa Rica 5331.2497 5273.0615 323.00325
## 7                           Estonia  782.8760  776.5314  55.08342
## 8                             Italy 2097.5834 2085.1940 133.16929
## 9  Lao People's Democratic Republic  576.1632  573.7042  39.61243
## 10                         Malaysia 1473.7752 1458.2762  94.28255
## 11                      Netherlands 2495.7160 2466.9962 140.25366
## 12                        Nicaragua 1569.0076 1567.9522  96.66135
## 13                        Singapore 1615.6165 1592.0610 108.28689
## 14                         Thailand 3207.9384 3154.9717 242.42609
## 15         United States of America 6125.7561 6052.1918 336.01728
## 16                          Vanuatu 1282.8242 1281.9231  80.43118

and even by county/region within a state or province (in this example, regions within British Columbia):

mybirds %>%
filter(State.Province == "CA-BC") %>%
group_by(County) %>%
do(mypd(.))
## Source: local data frame [15 x 4]
## Groups: County [15]
##
##                    County   mean_pd median_pd     sd_pd
##                     <chr>     <dbl>     <dbl>     <dbl>
## 1       Alberni-Clayoquot 1151.9561 1140.2276  63.73931
## 2                 Capital 2981.2589 2970.1906 160.05197
## 3        Central Okanagan 2367.7571 2347.5893 142.14835
## 4        Columbia-Shuswap 1442.1585 1444.1457  70.47375
## 5        Comox-Strathcona 1068.2720 1066.9558  44.22091
## 6           Fraser Valley 2447.5208 2436.1617 134.31824
## 7       Kootenay Boundary 1797.7464 1768.6904 116.87098
## 8         Metro Vancouver 5505.4194 5458.0431 315.90021
## 9        Mount Waddington 1232.4101 1223.7951  71.03491
## 10                Nanaimo  768.0895  762.9666  43.55813
## 11   Okanagan-Similkameen 3934.0560 3907.4926 234.13345
## 12 Skeena-Queen Charlotte  354.8324  350.6889  19.17975
## 13      Squamish-Lillooet 1825.5306 1803.6692 114.35689
## 14         Sunshine Coast 1528.9859 1511.2885  77.00035
## 15        Thompson-Nicola 2503.5555 2483.0836 147.10740

## Are ED and PD correlated?

Arne Mooers, who is a co-author of the bird tree paper, mentioned to me that there is often a strong correlation between ED and PD across species lists. With these two functions, it is easy to look at the relationship between PD and sum of EDs for a bunch of different lists, for example, ED and PD values by country:

# calculating ED by country
edcountry <- mybirds %>%
group_by(Country) %>%
do(myedge(.)) %>%
select(-sciName) %>%
summarise(EDs = sum(ED.Score), n = n())

# calculating PD by country
pdcountry <- mybirds %>%
group_by(Country) %>%
do(mypd(.))

edpdcountry <- left_join(edcountry, pdcountry, by = "Country")

# fitting linear model and extracting coefficients to print on plot
m <- lm(EDs ~ median_pd, edpdcountry)
rout <- list(paste('Fitted model: ', round(coef(m)[1], 3), ' + ',
round(coef(m)[2], 3), '* PD', sep = ''),
paste('R^2 == ', round(summary(m)[['r.squared']], 3),
sep = ''))

ggplot(data = edpdcountry, aes(median_pd, EDs, size = n)) +
geom_point(color = "grey60") +
ggtitle("ED and PD by country") +
geom_smooth(method = "lm", show.legend = FALSE) +
geom_text(aes(label = Country), size = 4,
nudge_y = -100, check_overlap = TRUE) +
geom_label(aes(x = 500, y = 3000, label = rout[[1]]), hjust = 0,
size = 5,  inherit.aes = FALSE) +
geom_label(aes(x = 500, y = 2600, label = rout[[2]]), hjust = 0, size = 5,
inherit.aes = FALSE, parse = TRUE)

At least in the example above, PD and ED are strongly correlated, with an $$R^2$$ of 0.98 and a coefficient of 0.45, which indicates that the sum of ED scores for a group of species is roughly half of the PD covered by that same group.

Lastly, the edge object provides a data frame that is useful for the two previous functions as it links the Jetz et al. (2014) phylogeny with your eBird checklist. You can peruse it independently, as it shows updated scientific names, common names, and EDGE scores and ranks as reported by the Jetz et al. paper:

edge
## # A tibble: 9,993 × 6
##                    sciName                       comName ED.Score
##                      <chr>                         <chr>    <dbl>
## 1       Pseudibis gigantea                    Giant Ibis 56.62592
## 2        Aegotheles savesi New Caledonian Owlet-nightjar 40.26083
## 3  Gymnogyps californianus             California Condor 33.39032
## 4      Strigops habroptila                        Kakapo 32.37154
## 5      Rhynochetos jubatus                          Kagu 55.38186
## 6  Houbaropsis bengalensis               Bengal Florican 24.27388
## 7          Athene blewitti                  Forest Owlet 24.08449
## 8    Pithecophaga jefferyi              Philippine Eagle 23.63732
## 9         Fregata andrewsi  Christmas Island Frigatebird 22.67818
## 10     Carpococcyx viridis        Sumatran Ground-cuckoo 21.15445
## # ... with 9,983 more rows, and 3 more variables: EDGE.Score <dbl>,
## #   EDGE.Rank <int>, sciName.edge <chr>

The species with the highest EDGE score is the Giant Ibis, which is classified as Critically Endangered by the IUCN (GE = 4) and also has the highest ED score among all Critically Endangered bird species.

Important note: There are many discrepancies between the taxonomy used in the Jetz paper and that one used in eBird; the latter even gets updated yearly! I’ve wrangled both species lists to try and match as many species as possible (you can see the code here), however the coverage is not perfect as there are many species which have been recently split and can’t be matched either by common or scientific name. Unfortunately, any species in your eBird checklists that can’t be matched to the Jetz phylogeny are ignored by myedge and mypd, yet there isn’t many of them (~50 out of 1300+ species in my checklist) however, this shouldn’t affect (too much) the species with highest EDGE scores or considerably bias PD estimates.

## Bugs/Problems?

This package is still a work in progress, so if you find any issues/bugs/problems please let me know. There are a few more features I’d like to add (e.g., graphing functions, plotting location of sightings on a map, to name a few); but any ideas are welcome!

## References

Isaac, N. J. B., Turvey, S. T., Collen, B., Waterman, C. & Baillie, J. E. M. (2007) Mammals on the EDGE: Conservation Priorities Based on Threat and Phylogeny. PLOS ONE, 2, e296. http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0000296

Jetz, W., Thomas, G. H., Joy, J. B., Redding, D. W., Hartmann, K. & Mooers, A. O. (2014) Global Distribution and Conservation of Evolutionary Distinctness in Birds. Current Biology, 24, 919–930. http://www.cell.com/current-biology/fulltext/S0960-9822(14)00270-X