Migratory Connectivity of a Neotropical Migratory Songbird Revealed by Archival Light-Level Geolocators

Michael T. Hallworth, T. Scott Sillett, Steven L. Van Wilgenburg, Keith A. Hobson and Peter P. Marra

accepted for publication in Ecological Applications

available here

Abstract

Understanding migratory connectivity is critical for interpreting population dynamics, seasonal interactions and for the implementation of conservation strategies of migratory species. We evaluated the migratory connectivity of a Neotropical migratory songbird, the Ovenbird (Seiurus aurocapilla) using archival light-level geolocators deployed at two breeding and four non-breeding locations while incorporating Ovenbird abundance as prior information using Bayes' Rule. We also included band recoveries submitted to the United States Geological Survey's Bird Banding Laboratory to assess connectivity of areas where geolocators were not deployed. We created a probabilistic map of origin for each capture site and mapped spring migration routes between non-breeding and breeding locations. We found a complete separation of eastern and western populations of Ovenbirds throughout the annual cycle. Breeding Ovenbirds from western Canada spent the non-breeding season throughout Central America and migrated through central North America during spring migration. Birds breeding in northeastern United States were distributed throughout the central Greater Antilles in the Caribbean and migrated through eastern North America during spring migration. Fall migration routes were not included because the timing of migration coincided with fall equinox when latitudinal estimates are unreliable. However, longitudinal estimates suggest no overlap between eastern and western populations during fall migration. Ovenbirds with geolocators attached in Jamaica bred in the northeastern United States with the highest posterior probability of origin found in Massachusetts, while Ovenbirds captured in Florida and Puerto Rico bred primarily in the mid-Atlantic. Incorporating Ovenbird abundance as a prior into geolocator estimates decreased the area of origin by 90.37% ± 1.05% for the breeding season and 62.30% ± 1.69% for the non-breeding season, compared to geolocator estimates alone. Ovenbirds exhibited strong migratory connectivity between breeding and non-breeding season which has important implications for various aspects of the ecology, evolution and conservation.


The following text provides sample code from our manuscript

relative abundance during the non-breeind season from eBird data

The abundance of Ovenbirds was estimated from Breeding Bird Survey data (Sauer et al. 2012) for the breeding season. Currently, no range wide monitoring data are available for the non-breeding distribution; therefore, we used eBird checklists (Avian Knowledge Network accessed 16 March 2013) reported from the Ovenbirds non-breeding range (see Ridgely et al. 2003) to determine non-breeding abundance.Non-breeding season abundance was generated using eBird checklists submitted during the non-breeding season between 1 November and 31 March 2010-2013.

Checklists submitted between 2010-2013 were pooled to increase the likelihood that sites were sampled multiple times between 1 November and 31 March. To maximize the amount of coverage across the non-breeding distribution we assumed no change in abundance between 2010 and 2013. Encounter histories were structured by month between November and March resulting in five sampling occasions at 17,103 sites.

eBirdData<-read.table("RoyleNichols_Abundance_covariates_NonSampled_removed.csv",header=TRUE,sep=",")
##       Site Nov_Oven Dec_Oven Jan_Oven Feb_Oven Mar_Oven Nov_num Dec_num
## 1 L1000342        0       NA        0       NA       NA       1      NA
## 2 L1000397       NA       NA       NA        0       NA      NA      NA
## 3 L1000680       NA       NA        0       NA       NA      NA      NA
## 4 L1000753        0       NA        0       NA        0       4      NA
## 5 L1001040       NA       NA       NA       NA        0      NA      NA
## 6 L1001219       NA       NA       NA        0       NA      NA      NA
##   Jan_num Feb_num Mar_num Nov_mins Dec_mins Jan_mins Feb_mins Mar_mins
## 1       2      NA      NA       20       NA       45       NA       NA
## 2      NA       2      NA       NA       NA       NA        0       NA
## 3       1      NA      NA       NA       NA       65       NA       NA
## 4       3      NA       2      105       NA       90       NA        0
## 5      NA      NA       1       NA       NA       NA       NA       20
## 6      NA       1      NA       NA       NA       NA        0       NA
##   LATITUDE LONGITUDE Elev NDVI_diff
## 1    18.17    -66.99  506  -0.05413
## 2    26.92    -82.27    5  -0.09055
## 3    26.57    -81.80   12  -0.04134
## 4    26.50    -81.96    6  -0.06594
## 5    26.12    -80.14   12  -0.05610
## 6    10.52    -84.98  697   0.03150
library(unmarked)
library(raster)
library(rgdal)
library(RColorBrewer)

Set up data to create an unmarked dataframe

Create “capture history”

History<-eBirdData[,c("Nov_Oven","Dec_Oven","Jan_Oven","Feb_Oven","Mar_Oven")]
##   Nov_Oven Dec_Oven Jan_Oven Feb_Oven Mar_Oven
## 1        0       NA        0       NA       NA
## 2       NA       NA       NA        0       NA
## 3       NA       NA        0       NA       NA
## 4        0       NA        0       NA        0
## 5       NA       NA       NA       NA        0
## 6       NA       NA       NA        0       NA

Create Unmarked data frame

Site covariates include Elevation, difference in NDVI score from Nov-Mar
Observation Covariates include number of checklists, total minutes birding

eBird.umf<-unmarkedFramePCount(y=History, 
siteCovs=data.frame(
Elev=eBirdData[,"Elev"],
NDVI=eBirdData[,"NDVI_diff"],
Long=eBirdData[,"LONGITUDE"]),
obsCovs=list(Ncounts=eBirdData[,c("Nov_num","Dec_num","Jan_num","Feb_num","Mar_num")],
CountMins=eBirdData[,c("Nov_mins","Dec_mins","Jan_mins","Feb_mins","Mar_mins")]) )

### Summary Data ###
summary(eBird.umf)
## unmarkedFrame Object
## 
## 17103 sites
## Maximum number of observations per site: 5 
## Mean number of observations per site: 1.33 
## Sites with at least one detection: 563 
## 
## Tabulation of y observations:
##     0     1     2     3     4     5     6     7     8     9    10    11 
## 22120   350   161    69    41    10    16     8    15     6     2     2 
##    12    13    14    16    18    20    23    24    35    46  <NA> 
##     3     1     2     2     1     1     1     1     1     1 62701 
## 
## Site-level covariates:
##       Elev           NDVI              Long       
##  Min.   : -54   Min.   :-0.3878   Min.   :-106.7  
##  1st Qu.:   7   1st Qu.:-0.0846   1st Qu.: -84.3  
##  Median :  14   Median :-0.0591   Median : -81.8  
##  Mean   : 220   Mean   :-0.0619   Mean   : -82.7  
##  3rd Qu.:  60   3rd Qu.:-0.0394   3rd Qu.: -80.6  
##  Max.   :3962   Max.   : 0.4104   Max.   : -59.5  
## 
## Observation-level covariates:
##     Ncounts        CountMins    
##  Min.   :  1     Min.   :    0  
##  1st Qu.:  1     1st Qu.:    0  
##  Median :  1     Median :   60  
##  Mean   :  2     Mean   :  220  
##  3rd Qu.:  2     3rd Qu.:  186  
##  Max.   :172     Max.   :29602  
##  NA's   :62701   NA's   :62701
## Standardize covariates to optimize model fit during MLE ##
obsCovs(eBird.umf)<-scale(obsCovs(eBird.umf))
siteCovs(eBird.umf)<-scale(siteCovs(eBird.umf))

Determine which distribution fits the data

dist.P   <-pcount(~1 ~1, eBird.umf, mixture = "P")    # Poisson 
dist.NB  <-pcount(~1 ~1, eBird.umf, mixture = "NB")   # Negative Binomial 
dist.ZIP <-pcount(~1 ~1, eBird.umf, mixture = "ZIP")  # Zero Inflated Poisson

DistMods<-fitList(
"Poisson"                 = dist.P,
"Negative Binomial"       = dist.NB,
"Zero Inflated Poisson"   = dist.ZIP)

#Rank models by AIC #

(DistMods_AIC<-modSel(DistMods))
##                       nPars      AIC   delta  AICwt cumltvWt
## Negative Binomial         3  8676.59    0.00  1e+00     1.00
## Zero Inflated Poisson     3  9817.31 1140.72 2e-248     1.00
## Poisson                   2 14021.30 5344.71  0e+00     1.00

Construct and fit models

fm1<- pcount(~1 ~1, eBird.umf, mixture = "NB")
fm2<- pcount(~CountMins ~1, eBird.umf, mixture = "NB")
fm3<- pcount(~Ncounts ~1, eBird.umf, mixture = "NB")
fm4<- pcount(~CountMins ~Elev, eBird.umf, mixture = "NB")
fm5<- pcount(~CountMins ~NDVI, eBird.umf, mixture = "NB")
fm6<- pcount(~1 ~Elev+NDVI, eBird.umf, mixture = "NB")
fm7<- pcount(~1 ~I(Elev^2)+NDVI, eBird.umf, mixture = "NB")
fm8<- pcount(~1 ~I(Long^2), eBird.umf, mixture = "NB")
fm9<- pcount(~1 ~I(Long^2)+Elev, eBird.umf, mixture = "NB")
fm10<- pcount(~1 ~NDVI, eBird.umf, mixture = "NB")
fm11<- pcount(~CountMins ~I(Elev^2)+NDVI, eBird.umf, mixture = "NB")
fm12<- pcount(~NDVI ~I(Elev^2)+NDVI, eBird.umf, mixture = "NB")

Show results

(OvenModsRank<-modSel(OvenMods))

#                             nPars     AIC  delta   AICwt cumltvWt
#lam(Elev^2+NDVI)p(CountMins)     6 8242.74   0.00 8.6e-01     0.86
#lam(NDVI)p(CountMins)            5 8246.37   3.63 1.4e-01     1.00
#lam(.)p(CountMins)               4 8269.57  26.83 1.3e-06     1.00
#lam(Elev)p(CountMins)            5 8271.54  28.80 4.8e-07     1.00
#lam(.)p(Ncounts)                 4 8380.39 137.65 1.1e-30     1.00
#lam(Elev^2+NDVI)p(NDVI)          6 8618.22 375.47 2.5e-82     1.00
#lam(Elev^2+NDVI)p(.)             5 8653.27 410.53 6.1e-90     1.00
#lam(NDVI)p(.)                    4 8657.42 414.68 7.7e-91     1.00
#lam(Elev+NDVI)p(.)               5 8659.27 416.53 3.1e-91     1.00
#lam(.)p(.)                       3 8676.59 433.85 5.3e-95     1.00
#lam(Long^2)p(.)                  4 8678.41 435.67 2.1e-95     1.00
#lam(Long^2 + Elev)p(.)           5 8680.33 437.59 8.2e-96     1.00

Create new data with the predicted values

# Save the mean and sd of variables
mean.Elev<-mean(eBirdData$Elev)
sd.Elev<-sd(eBirdData$Elev)
mean.NDVI<-mean(eBirdData$NDVI)
sd.NDVI<-sd(eBirdData$NDVI)

# Create new data with predicted values
new.Elev<- data.frame(Elev=seq(-0.536,7.32,,100),NDVI=seq(-3.859,5.595,,100))
pred.Elev<-predict(fm11,type="state",newdata=new.Elev,appendData=TRUE)

new.Mins<-data.frame(CountMins=seq(-0.316,42.239,,17103))
pred.Mins<-predict(fm11,type="det",newdata=new.Mins,appendData=TRUE)
##   doing row 1000 of 17103 
##   doing row 2000 of 17103 
##   doing row 3000 of 17103 
##   doing row 4000 of 17103 
##   doing row 5000 of 17103 
##   doing row 6000 of 17103 
##   doing row 7000 of 17103 
##   doing row 8000 of 17103 
##   doing row 9000 of 17103 
##   doing row 10000 of 17103 
##   doing row 11000 of 17103 
##   doing row 12000 of 17103 
##   doing row 13000 of 17103 
##   doing row 14000 of 17103 
##   doing row 15000 of 17103 
##   doing row 16000 of 17103 
##   doing row 17000 of 17103
pred.Mins$Minutes<-(pred.Mins$CountMins*1391.692)+293.4254
pred.Mins$Hours<-pred.Mins$Minutes/60
pred.Elev$Elevation<-(pred.Elev$Elev*sd.Elev)+mean.Elev
pred.Elev$NDVI_diff<-(pred.Elev$NDVI*sd.NDVI)+mean.NDVI
head(pred.Elev)
##   Predicted     SE lower upper    Elev   NDVI Elevation NDVI_diff
## 1     4.191 1.0437 2.572 6.828 -0.5360 -3.859    -53.99   -0.3878
## 2     4.084 0.9945 2.534 6.582 -0.4566 -3.764    -13.43   -0.3797
## 3     3.977 0.9465 2.494 6.341 -0.3773 -3.668     27.13   -0.3716
## 4     3.871 0.8999 2.454 6.105 -0.2979 -3.573     67.69   -0.3636
## 5     3.765 0.8547 2.413 5.875 -0.2186 -3.477    108.25   -0.3555
## 6     3.660 0.8109 2.371 5.650 -0.1392 -3.382    148.81   -0.3474

Generate figure to display expected results of model

plot of chunk unnamed-chunk-14

Create abundance map using Beta estimates of best model (fm11)

plot of chunk unnamed-chunk-17

Geolocator Assignments

See text for information regarding transition events and sun-elevation angles. The following text begins with kernel density estimates of the stationary periods of the annual cycle (breeding/non-breeding)

Load the required packages

library(raster)
library(rgdal)
library(RColorBrewer)

read in prior distribution Ovenbird relative abundance during the non-breeding season
This file was made using eBird checklists and the unmarked package (see above)

eBird_Abundance<-raster("nmix_predict.txt")
NB_dist<-shapefile("OVEN_NonBreeding_Dist_WGS.shp")        #Non-breeding distribution shapefile

Clip the abundance map to the size of the Non-breeding distribution

eBird_Abun_dist<-mask(eBird_Abundance,NB_dist)
## Found 1 region(s) and 786 polygon(s)

A few cells in the raster were incorrect, replaced cell with NA

eBird_Abun_dist[(eBird_Abun_dist==Inf)]<-NA

Create a probability surface for the prior surface by dividing each cell in the surface by the sum of the surface

cellStats(eBird_Abun_dist,sum)
## [1] 19938
eBird_Abun_Prob<-calc(eBird_Abun_dist, function(x){x/cellStats(eBird_Abun_dist,sum)})

Read in the kernel density estimates for the Non-breeding season and make a raster of the files

KDE_names<-list.files("NonBreeding_Clipped_KDE/SnapRaster_WGS", pattern="*_clip.txt",full.names=TRUE)
# read in file names #
KDE_listnames<-list.files("NonBreeding_Clipped_KDE/SnapRaster_WGS", pattern="*_clip.txt",full.names=FALSE)
# Make raster layer from folder #
KDE_layer<-lapply(KDE_names, raster)

Extend each KDE layer to the relative abundance layer

KDE_extend<-list(
kde25044<-extend(KDE_layer[[1]],eBird_Abun_dist,value=0),
kde25060<-extend(KDE_layer[[2]],eBird_Abun_dist,value=0),
kde25061<-extend(KDE_layer[[3]],eBird_Abun_dist,value=0),
...)

Crop the KDE layers to the non-breeding distribution

KDE_crop<-lapply(KDE_extend,function(x){crop(x,eBird_Abun_dist)})

Bayesian Assignment using abundance model as prior information

Bayes_Numerator_Nmix<-lapply(KDE_crop,function(x){overlay(x,eBird_Abun_dist, fun=function(x,y){return(x*y)})})

for(i in 1:45){
Bayes_assignment_Nmix<-lapply(Bayes_Numerator_Nmix,function(x){calc(x,fun=function(x) {x/cellStats(Bayes_Numerator_Nmix[[i]],sum)})})
}

Plot non-breeding assignments using non-breeding abundance as prior information

plot of chunk unnamed-chunk-30