Michael T. Hallworth and Peter P. Marra
For the first time, we use a small archival global positioning system (GPS) tag to identify and characterize non-breeding territories, quantify migratory connectivity, and identify population boundaries of Ovenbirds (Seiurus aurocapilla), a small migratory songbird, captured at two widely separated breeding locations. We recovered 15 (31 %) GPS tags with data and located the non-breeding territories of breeding Ovenbirds from Maryland and New Hampshire, USA (0.50 ± 0.15 ha, mean ± SE). All non-breeding territories had similar environmental attributes despite being distributed across parts of Florida, Cuba and Hispaniola. New Hampshire and Maryland breeding populations had non-overlapping non-breeding population boundaries that encompassed 114,803 and 169,233 km2, respectively. Archival GPS tags provided unprecedented pinpoint locations and associated environmental information of tropical non-breeding territories. This technology is an important step forward in understanding seasonal interactions and ultimately population dynamics of populations throughout the annual cycle.
Hallworth, M. T. and P. P. Marra. 2015. Miniaturized GPS Tags Identify Nonbreeding Territories of a Small Breeding Migratory Songbird. Scientific Reports 5:11069 doi: 10.1038/srep11069
The following text provides code for the analyses in our manuscript
# Load required packages for anaylsis #
library(ks)
library(raster)
library(sp)
library(rgdal)
library(R2jags)
library(ade4)
library(sp)
library(RColorBrewer)
library(dismo)
library(rgeos)
Read data file with locations and environmental data from MoveBank.org. A few rows and columns of the data are displayed below.
# Location data and attributes #
GPSdata<-read.csv("Data/Ovenbird_GPS_MoveBank_HallworthMT_EnvData.csv")
# Spatial Layers #
worldmap<-shapefile("Spatial_Layers/TM_WORLD_BORDERS-0.3.shp")
states<-shapefile("Spatial_Layers/st99_d00.shp")
## Latitude Longitude Duration DOP Satellites BirdID Date Time
## 1 43.94176 -71.73830 31 3.0 5 2521_45043 7/1/2013 1:40:18
## 2 43.94037 -71.73757 37 2.2 5 2521_45043 7/29/2013 1:38:42
## 3 43.94147 -71.73782 36 3.0 4 2521_45043 8/26/2013 1:36:42
## 4 42.97363 -72.01157 34 1.4 8 2521_45043 9/23/2013 1:34:36
## 5 19.71088 -70.79377 42 2.2 5 2521_45043 10/21/2013 1:32:48
## 6 19.23941 -69.94896 35 3.0 4 2521_45043 11/18/2013 1:30:43
Create spatial points data frame in order to do spatial analyses. Location data were collected in WGS84 projection (or lack of a projection)
coords<-cbind(GPSdata[,2],GPSdata[,1])
GPSspdf<-SpatialPointsDataFrame(coords,GPSdata,proj4string=(CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0")))
## class : SpatialPointsDataFrame
## features : 145
## extent : -84.35636, -69.23454, 18.31352, 43.94719 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0
## variables : 27
## names : Latitude, Longitude, Duration, DOP, Satellites, BirdID, Date, Time, CaptureLocation, utm_easting, utm_northing, utm_zone, NASA.Distance.to.Coast, ASTER.ASTGTM2.Elevation, TRMM.Combined.Precipitation, ...
## min values : 18.31352, -69.23454, 20, 0.8, 3, 1831_95478, 1/13/2014, 1:01:12, HBEF, 256288.3, 2026253, 16N, 3.987785, 8.269577, 0.000000000, ...
## max values : 43.94719, -84.35636, 66, 34.4, 10, 2521_45068, 9/23/2013, 2:03:48, JB, 789563.1, 4869587, 19N, 297.486469, 1043.294549, 1.632081388, ...
# Combine latitude and longitude for each bird #
WinterTerrs<-cbind(GPSdata[c(6:10,16:20,25:30,35:36,41:43,50:55,61:65,69:75,
79:85,90:95,100:105,109:115,119:125,129:135,139:145),2],
GPSdata[c(6:10,16:20,25:30,35:36,41:43,50:55,61:65,69:75,
79:85,90:95,100:105,109:115,119:125,129:135,139:145),1],
GPSdata[c(6:10,16:20,25:30,35:36,41:43,50:55,61:65,69:75,
79:85,90:95,100:105,109:115,119:125,129:135,139:145),6],
GPSdata[c(6:10,16:20,25:30,35:36,41:43,50:55,61:65,69:75,
79:85,90:95,100:105,109:115,119:125,129:135,139:145),9])
# determine the dates #
winterdates<-GPSdata[c(6:10,16:20,25:30,35:36,41:43,50:55,61:65,69:75,
79:85,90:95,100:105,109:115,119:125,129:135,139:145),7]
# combine into single data frame #
birds<-GPSdata[c(6:10,16:20,25:30,35:36,41:43,50:55,61:65,69:75,
79:85,90:95,100:105,109:115,119:125,129:135,139:145),c(6,7,9)]
# make empty vector to store bandwidth estimator - least-square cross validation
HlscvBird<-rep(NA,17)
# loop through birds to determine bandwith
for(i in 3:15){
sel.rows<-WinterTerrs[,3]==i
HlscvBird[i]<-Hlscv(WinterTerrs[sel.rows,1:2])
}
HlscvBird[16:17]<-mean(HlscvBird,na.rm=TRUE)
# make an empty list to store territories
BirdTerr<-list()
# loop through birds to generate territories
for(i in 3:17){
sel.rows=WinterTerrs[,3]==i
BirdTerr[i]<-raster(kde(WinterTerrs[sel.rows,1:2],h=HlscvBird[i]))
}
## Error in raster(kde(WinterTerrs[sel.rows, 1:2], h = HlscvBird[i])): error in evaluating the argument 'x' in selecting a method for function 'raster': Error in chol.default(S) :
## the leading minor of order 2 is not positive definite
## Calls: kde ... pre.sphere -> matrix.sqrt -> chol2inv -> chol -> chol.default
# Scale kernel density estimate between 0-1 for each bird
TerrScale<-list()
for(i in 3:15){
TerrScale[i]<-calc(BirdTerr[[i]],fun=function(x){x/cellStats(BirdTerr[[i]],max)})
}
# Convert raster to polygon
TerrPoly<-list()
for(i in 3:15){
TerrPoly[i]<-rasterToPolygons(TerrScale[[i]],fun=function(x){x>0.049999})
}
# Simplify the polygon
TerrrelPolyD<-list()
for(i in 3:15){
TerrrelPolyD[i]<-gUnaryUnion(TerrPoly[[i]])
}
# project the polygon in Lambert Conic to get an area in hectares
PolyMeters<-list()
for(i in 3:15){
PolyMeters[i]<-spTransform(TerrrelPolyD[[i]],CRS("+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))
}
# determine area of individual territories
areaPoly<-rep(NA,15)
for(i in 3:15){
areaPoly[i]<-gArea(PolyMeters[[i]])/10000
}
# Generate the mean territory size
mean(areaPoly,na.rm=TRUE)
## [1] 0.496146
We used a Mantel test to determine the strength of migratory connectivity between breeding and non-breeding locations. Values close to 1 indicate a strong correlation between breeding and non-breeding distance matrices (i.e., individuals breeding in close proximity also winter in close proximity). We used the first location taken as the breeding ground location -(all GPS tags took at least 1 location on breeding grounds) and the final location as the non-breeding location.
# Mantel Test - All birds #
Breed_dist<-dist(cbind(GPSdata[c(1,11,21,31,37,46,56,66,76,86,96,106,116,126,136),1],
GPSdata[c(1,11,21,31,37,46,56,66,76,86,96,106,116,126,136),2]))
NB_dist<-dist(cbind(GPSdata[c(10,20,30,36,43,55,65,75,85,95,105,115,125,135,145),1],
GPSdata[c(10,20,30,36,43,55,65,75,85,95,105,115,125,135,145),2]))
Mantel.result1<-mantel.rtest(Breed_dist,NB_dist,nrepet=10000)
## Monte-Carlo test
## Observation: 0.8393364
## Call: mantel.rtest(m1 = Breed_dist, m2 = NB_dist, nrepet = 10000)
## Based on 10000 replicates
## Simulated p-value: 0.00029997
### Mantel Test - HBEF birds only ###
HBEFBreed_dist<-dist(cbind(GPSdata[c(1,11,21,31,37),1],
GPSdata[c(1,11,21,31,37),2]))
HBEFNB_dist<-dist(cbind(GPSdata[c(10,20,30,36,43),1],
GPSdata[c(10,20,30,36,43),2]))
HBEFMantel.result1<-mantel.rtest(HBEFBreed_dist,HBEFNB_dist,nrepet=10000)
### Mantel Test - JBWS birds only ###
JBWSBreed_dist<-dist(cbind(GPSdata[c(46,56,66,76,86,96,106,116,126,136),1],
GPSdata[c(46,56,66,76,86,96,106,116,126,136),2]))
JBWSNB_dist<-dist(cbind(GPSdata[c(55,65,75,85,95,105,115,125,135,145),1],
GPSdata[c(55,65,75,85,95,105,115,125,135,145),2]))
JBWSMantel.result1<-mantel.rtest(JBWSBreed_dist,JBWSNB_dist,nrepet=10000)
Birds breeding within Hubbard Brook Experimental Forest, NH exhibited stronger connectivity than birds breeding at Jug Bay Wetland Sanctuary, MD.
HBEFMantel.result1
## Monte-Carlo test
## Observation: 0.5893957
## Call: mantel.rtest(m1 = HBEFBreed_dist, m2 = HBEFNB_dist, nrepet = 10000)
## Based on 10000 replicates
## Simulated p-value: 0.01619838
JBWSMantel.result1
## Monte-Carlo test
## Observation: -0.1159604
## Call: mantel.rtest(m1 = JBWSBreed_dist, m2 = JBWSNB_dist, nrepet = 10000)
## Based on 10000 replicates
## Simulated p-value: 0.5906409
We define boundaries during the non-breeding season as the 95% kernel density estimate around the non-breeding locations for each breeding population.
# Use only the final non-breeding location in analysis #
WinterPoints<-GPSspdf[c(10,20,30,36,43,55,65,75,85,95,105,115,125,135,145),]
# Non breeding locations #
NBlocations<-cbind(GPSdata[c(10,20,30,36,43,55,65,75,85,95,105,115,125,135,145),2],
GPSdata[c(10,20,30,36,43,55,65,75,85,95,105,115,125,135,145),1],
GPSdata[c(10,20,30,36,43,55,65,75,85,95,105,115,125,135,145),8],
GPSdata[c(10,20,30,36,43,55,65,75,85,95,105,115,125,135,145),9])
# Estimate bandwidth using LSCV #
bwidthNH<-Hlscv(NBlocations[1:5,1:2])
bwidthMD<-Hlscv(NBlocations[6:15,1:2])
# Create a raster layer for each kernel density estimate #
NHpopBoundary<-raster(kde(NBlocations[1:5,1:2],h=bwidtNH))
MDpopBoundary<-raster(kde(NBlocations[6:15,1:2],h=bwidthMD))
# Determine the maximum value within the raster surface #
NHmax<-cellStats(NHpopBoundary,max)
MDmax<-cellStats(MDpopBoundary,max)
# Scale the surface to values between 0-1 #
NHscaleBoundary<-calc(NHpopBoundary,fun=function(x){x/NHmax})
MDscaleBoundary<-calc(MDpopBoundary,fun=function(x){x/MDmax})
# Convert cells with values < 0.05 to NA #
NHscaleBoundary[NHscaleBoundary<0.05]<-NA
MDscaleBoundary[MDscaleBoundary<0.05]<-NA
# Only keep areas of the population boundary that fall over land #
NHrelBoundary<-mask(NHscaleBoundary,worldmap)
MDrelBoundary<-mask(MDscaleBoundary,worldmap)
# Convert raster to polygon #
NHrelPoly<-rasterToPolygons(NHscaleBoundary,fun=function(x){x>0.049999})
MDrelPoly<-rasterToPolygons(MDscaleBoundary,fun=function(x){x>0.049999})
# Merge into a single polygon #
NHrelPolyD<-gUnaryUnion(NHrelPoly)
MDrelPolyD<-gUnaryUnion(MDrelPoly)
# *NOTE* This analysis takes quite a while to run - not run here
# Iterate through points to see if KDE reach an asymptote #
# Jug Bay Wetland Sanctuary - MARYLAND
# subset data to include only JBWS non-breeding locations #
JBnonbreed<-GPSdata[c(55,65,75,85,95,105,115,125,135,145),c(2,1)]
# Make empty list to store results #
bwidthMD<-MDpopBoundary<-probKDE<-KDE95<-AREA<-AreaLand<-list()
SUMarea<-rep(NA,10)
#Generate all the combinations of individuals to run through #
numBirds<-seq(from=1,to=10,by=1)
Combos<-list()
for(i in 2:10){
Combos[[i]]<-combn(numBirds,i,simplify=FALSE)
}
SIZEarray<-array(NA,c(252,10))
# Iterate through all the possible combinations and generate population size #
for(n in 3:10){
for(i in 1:length(Combos[[n]])){
bwidthMD[[i]]<-Hscv(JBnonbreed[Combos[[n]][[i]],1:2])
MDpopBoundary[[i]]<-raster(kde(JBnonbreed[Combos[[n]][[i]],],h=bwidthMD[[i]]))
probKDE[[i]]<-calc(MDpopBoundary[[i]],fun=function(x){x/cellStats(MDpopBoundary[[i]],max)})
KDE95[[i]]<-calc(probKDE[[i]],fun=function(x){x>=0.05})
AREA[[i]]<-overlay(KDE95[[i]],area(MDpopBoundary[[i]]),fun=function(x,y){return(x*y)})
AreaLand[[i]]<-mask(AREA[[i]],worldmap)
SIZEarray[i,n]<-cellStats(AreaLand[[i]],sum)
}
}
# Hubbard Brook Experimental Forest, NEW HAMPSHIRE #
# Iterate through points to see if KDE reach an asymptote #
NHnonbreed<-GPSdata[c(10,20,30,36,43),c(2,1)]
bwidthNH<-NHpopBoundary<-NHprobKDE<-NHKDE95<-NHAREA<-NHAreaLand<-list()
#Generate all the combinations of individuals to run through #
NHnumBirds<-seq(from=1,to=5,by=1)
NHCombos<-list()
for(i in 2:5){
NHCombos[[i]]<-combn(NHnumBirds,i,simplify=FALSE)
}
NHSIZEarray<-array(NA,c(252,10))
for(n in 3:5){
for(i in 1:length(NHCombos[[n]])){
bwidthNH[[i]]<-Hscv(NHnonbreed[NHCombos[[n]][[i]],1:2])
NHpopBoundary[[i]]<-raster(kde(NHnonbreed[NHCombos[[n]][[i]],],h=bwidthNH[[i]]))
NHprobKDE[[i]]<-calc(NHpopBoundary[[i]],fun=function(x){x/cellStats(NHpopBoundary[[i]],max)})
NHKDE95[[i]]<-calc(NHprobKDE[[i]],fun=function(x){x>=0.05})
NHAREA[[i]]<-overlay(NHKDE95[[i]],area(NHpopBoundary[[i]]),fun=function(x,y){return(x*y)})
NHAreaLand[[i]]<-mask(NHAREA[[i]],worldmap)
NHSIZEarray[i,n]<-cellStats(NHAreaLand[[i]],sum)
}
}
# Determine the mean size for JBWS #
meanSize<-apply(SIZEarray,2,mean,na.rm=TRUE)
sdSize<-apply(SIZEarray,2,sd,na.rm=TRUE)
seSize<-apply(SIZEarray,2,FUN=function(x){sd(x,na.rm=TRUE)/sqrt(length(x))})
# Determine the mean size of HBEF #
NHmeanSize<-apply(NHSIZEarray,2,mean,na.rm=TRUE)
NHsdSize<-apply(NHSIZEarray,2,sd,na.rm=TRUE)
NHseSize<-apply(NHSIZEarray,2,FUN=function(x){sd(x,na.rm=TRUE)/sqrt(length(x))})
1) Compare latitutde estimates on September 23rd 2013 between capture locations - Note September 23rd coincides with fall equinox, a time when light-level geolocator estimates for latitude are unreliable.
# Subset data set to only show where birds were on September 23rd 2013 #
Sept23<-subset(GPSdata,Date=="9/23/2013")
n<-length(Sept23[,6]) # number of individuals
Location<-as.numeric(as.factor(Sept23[,9]))
Latitude<-Sept23[,1]
# T-test assuming unequal variance using JAGS #
Latitude_Ttest<-function()
{
#Priors
mu1~dunif(0,20)
mu2~dunif(0,20)
tau1<-1/(sigma1*sigma1)
tau2<-1/(sigma2*sigma2)
sigma1~dunif(0,100)
sigma2~dunif(0,100)
#Likelihood
for (i in 1:n1){
y1[i]~dnorm(mu1,tau1)
}
for (i in 1:n2){
y2[i]~dnorm(mu2,tau2)
}
# Derived parameters
delta<-mu2-mu1
} # End Model
win.data<-list(y1=abs(Sept23[1:5,1]-43.94), y2=abs(Sept23[6:16,1]-38.77),n1=4,n2=11)
params<-c("mu1","mu2","delta")
LatTtest<-jags(model=Latitude_Ttest,
data=win.data,
parameters.to.save=params,
n.chains=3,
n.iter=50000,
n.burnin=10000,
n.thin=10)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 31
##
## Initializing model
## Inference for Bugs model at "C:/Users/MTHALL~1/AppData/Local/Temp/RtmpSi3C4X/model11b815f42fcb.txt", fit using jags,
## 3 chains, each with 50000 iterations (first 10000 discarded), n.thin = 10
## n.sims = 12000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## delta 5.965 3.594 -2.551 3.969 6.296 8.342 12.020 1.001 6500
## mu1 4.690 2.881 0.545 2.753 4.250 6.002 12.168 1.001 4600
## mu2 10.655 2.128 6.335 9.337 10.691 12.001 14.820 1.001 12000
## deviance 95.619 3.867 90.643 92.786 94.756 97.575 105.203 1.001 12000
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 7.5 and DIC = 103.1
## DIC is an estimate of expected predictive error (lower deviance is better).
2) Does migration distance differ between the two populations?
# Generate straight line migration distance #
NHdist<-as.matrix(dist(cbind(GPSdata[c(1,11,21,31,37,10,20,30,36,43),11],
GPSdata[c(1,11,21,31,37,10,20,30,36,43),10])))
MDdist<-as.matrix(dist(cbind(GPSdata[c(46,56,66,76,86,96,106,116,126,136,55,65,75,
85,95,105,115,125,135,145),11],
GPSdata[c(46,56,66,76,86,96,106,116,126,136,55,65,75,
85,95,105,115,125,135,145),10])))
# Determine distances *NOTE* conversion from m to km #
NHdistances<-rep(NA,5)
for(i in 1:5){
NHdistances[i]<-NHdist[i,5+i]/1000
}
MDdistances<-rep(NA,10)
for(i in 1:10){
MDdistances[i]<-MDdist[i,10+i]/1000
}
# T-test assuming unequal variance using JAGS #
Distance_Ttest<-function()
{
#Priors
mu1~dunif(0,5000)
mu2~dunif(0,5000)
tau1<-1/(sigma1*sigma1)
tau2<-1/(sigma2*sigma2)
sigma1~dunif(0,100)
sigma2~dunif(0,100)
#Likelihood
for (i in 1:n1){
y1[i]~dnorm(mu1,tau1)
}
for (i in 1:n2){
y2[i]~dnorm(mu2,tau2)
}
#Derived Parameters
delta<-mu2-mu1
} # End model
win.data<-list(y1=NHdistances, y2=MDdistances,n1=5,n2=10)
params<-c("mu1","mu2","delta")
DistTtest<-jags(model=Distance_Ttest,
data=win.data,
parameters.to.save=params,
n.chains=3,
n.iter=50000,
n.burnin=10000,
n.thin=10)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 30
##
## Initializing model
## Inference for Bugs model at "C:/Users/MTHALL~1/AppData/Local/Temp/RtmpC2JYOc/model10ac378d40db.txt", fit using jags,
## 3 chains, each with 50000 iterations (first 10000 discarded), n.thin = 10
## n.sims = 12000 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75%
## delta -1179.596 47.331 -1271.431 -1211.670 -1179.818 -1147.579
## mu1 2767.684 35.379 2697.353 2744.644 2767.909 2790.687
## mu2 1588.088 31.396 1526.890 1566.960 1587.709 1609.252
## deviance 239.664 2.800 236.051 237.612 239.052 241.080
## 97.5% Rhat n.eff
## delta -1085.992 1.001 12000
## mu1 2837.566 1.001 12000
## mu2 1649.307 1.001 12000
## deviance 246.643 1.001 12000
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 3.9 and DIC = 243.6
## DIC is an estimate of expected predictive error (lower deviance is better).
3) Determine Land-Cover Classification of non-breeding locations
# read in land-cover data set
LandCover<-raster("Spatial_Layers/global_2012_5min.asc")
# Land-cover Type #
#Value Label
#0 Water
#1 Evergreen Needleleaf forest
#2 Evergreen Broadleaf forest
#3 Deciduous Needleleaf forest
#4 Deciduous Broadleaf Forest
#5 Mixed Forest
#6 Closed shrublands
#7 Open shrublands
#8 Wood Savanas
#9 Savannas
#10 Grasslands
#11 Permanent wetlands
#12 Croplands
#13 Urban and built-up
#14 Cropland/Natural vegetation mosaic
#15 Snow and Ice
#16 Barren or sparsely vegetated
#254 Unclassified
#255 Fill Value
# Extract land-cover class to locations
WinterLandClass<-extract(LandCover,WinterPoints)
#HBEF land-cover classes
HBEFtable
## Evergreen Broadleaf forest Croplands
## 2 1
## Cropland/Natural vegetation mosaic
## 2
# JBWS land-cover classes
JBWStable
## Wood Savanas Permanent wetlands
## 3 1
## Cropland/Natural vegetation mosaic
## 6