Progress report for GNE20-233
Project Information
Crop damage is a persistent issue faced by farmers nationwide, resulting in staggering costs to farmers in crop loss. Ungulate species such as white-tailed deer (Odocoileus virginianus) and sika deer (Cervus nippon) are a major source of agricultural damage through extensive browsing. Often a farmer’s only solution to mitigate browsing losses is through crop damage permits to hunt deer on their farm. However, deer populations occur over larger scales than a single property, limiting farmers’ abilities to actually reduce local deer densities through hunting and leaving them ineffective in reducing crop damage by ungulates. Surrounding habitats can reduce or exacerbate impacts on crop damage – reducing damage by providing additional food that lowers the intensity of deer browsing on crops or exacerbating damage by providing refuge to deer that actively browse on agricultural fields. My research looks to examine both the effect that adjacent habitat composition and configuration has on crop damage as well as quantify the impact that dietary overlap of two co-occurring ungulates has on crop damage. High-resolution dietary analysis using metabarcoding plant DNA present in ungulate fecal pellets allows crop damage to be quantified. Combining the dietary profile of ungulates with landscape level metrics of the surrounding habitats will reveal how habitats in close proximity to agriculture influence crop damage. Through informing famers, policy makers and state agencies about the influence of adjacent habitats on crop damage, this work will develop management polices beyond deer culling to effectively mitigate crop damage.
The specific objectives of this research are to:
(1) Determine the impact that the composition and configuration of neighboring habitats (forests, marshes, fallow fields) within agricultural landscapes have on the proportion of agricultural crops in a deer’s diet.
Hypothesis 1: Overall area, number of patches, and connectivity of neighboring habitats will decrease the proportion of agricultural species present in the diet of the deer.
Hypothesis 2: The proportion of agricultural species present in the diet of the deer can be explained by the composition and configuration of neighboring habitats.
Hypothesis 3: Reliance on plant species present in neighboring habitats will vary through time depending upon the abundance and availability of crops.
(2) Evaluate the impact of dietary niche overlap of sympatric deer species has on crop damage.
Hypothesis 1: We will observe higher overall deer densities in regions where both deer species occur than regions where just white-tailed deer or sika deer occur in isolation.
Hypothesis 2: Agricultural species will make up a smaller proportion of the diet of white-tailed deer and sika deer in regions where they co-occur due to interspecific competition for crop resources.
Hypothesis 3: Despite lower per capita crop consumption, crop damage will be greater in areas where the species co-occur due to higher densities.
Hypothesis 4: Higher deer densities will be supported by a greater reliance on plant species from neighboring habitats, and we expect to see a wider dietary breadth in regions where both species co-occur.
* Objective 2 had to be modified due to restricted access to farms in regions which both species co-occur. COVID-19 limited our ability to interact with farmers in these regions as well as farms that were possible site locations being sold to new owners that did not want to participate in our study. We were only able to gain access to two farm sites in which both species co-occur resulting in not a large enough sample size to proceed with this objective. I instead reallocated this effort to pursue a related question on the spatial distribution of crop damage at the patch and landscape-level. Objective 3 below outlines this new line of enquiry.
(3) Understand landscape and patch dynamics that influence the spatial distribution of crop damage at the patch and landscape-level
Hypothesis 1: Edge habitat composition will influence the prevalence and overall distance crop damage occurs into the agricultural field (patch-level).
Hypothesis 2: Field configuration will impact the level of crop damage occurring within the agricultural field (patch-level).
Hypothesis 3: Composition and configuration of neighboring habitats will influence the proportion of crop damage the entire agricultural field experiences (landscape-level).
Hypothesis 4: Neighboring habitats not directly adjacent but interconnected to the agricultural field will impact the proportion of damage the entire agricultural field experiences (landscape-level).
Hypothesis 5: Higher proportions of agriculture in the landscape will result in lower levels of crop damage across the entire agricultural field (landscape-level).
The purpose of this project is to examine the effect that landscape configuration and composition of surrounding habitat types has on crop damage by deer. Crop damage is a persistent issue faced by farmers nationwide; surveys have shown that crop damage inflicted by deer is the most widespread form of wildlife damage to crops (Conover and Decker 2018; Wywialowsky 1994, Conover 1998). Ungulates, specifically white-tailed deer (Odocoileus virginanus) and sika deer (Cervus nippon), damage crops through browsing and are a major source of agricultural damage, both of which occur on the Eastern Shore of Maryland (Claslick and Decker 1979, Conover 2001, Conover and Decker 2018, DeVault et al. 2007, Springer et al. 2013). Often a farmer’s only solution to mitigate crop damage is through crop damage permits that allow the increased harvest of deer on their farms. However, culling deer on-site may not solve the problem, as deer populations occur at larger scales than a single property. Although crop damage by deer often correlates with deer population density, not all deer present in an agricultural landscape consume the same proportion of agricultural species (Bleier et al. 2012). Agricultural fields do not occur in isolation but are surrounded by a complex matrix of habitats. Deer frequently utilize habitats surrounding agricultural fields. How do surrounding habitats affect deer consumption of crops? Surrounding habitats have been suggested to negatively and positively impact the amount of crop damage occurring on the adjacent fields. Deer can forage on plant species present in habitats surrounding agricultural fields mitigating deer browsing on crops, or deer can utilize these habitats as refuge resulting increased deer densities and increased browsing pressure on crops. However, no research to date has tried to directly quantify how surrounding habitats mitigate or intensify crop damage.
My research will utilize a relatively new methodology of DNA metabarcoding of deer fecal pellets to construct of high-resolution diet profiles of ungulate species in agricultural landscapes. The DNA barcoding approach will allow me to identify all plant species, crop and non-crop, and their relative contribution to deer diets. Quantifying crop damage inflicted by the deer present in the agricultural landscape on an individualistic level. Combining the dietary profiles of deer on agricultural fields with landscape metrics such as patch size, patch density, and patch connectivity of adjacent habitats will allow for a better understanding of how these neighboring habitats to agriculture affect crop damage. Data on how adjacent habitats influence crop damage will provide management recommendations to mitigate crop damage for farmers and policy makers. As many agricultural landowners in my study region own adjacent habitats as well as their fields, there will be a high opportunity for implementing deer management in adjacent habitats. Strategic community partners (see letters of support) will help me to develop recommendations for best practices based on my findings and disseminate results to land owners beyond the few sites where the research will be conducted.
Research
The following objectives will be accomplished via field experiments established in Dorchester, Wicomico, Worcester, and Somerset Counties, Maryland and lab work conducted at George Washington University.
Objective 1. Evaluate the impact that the configuration and composition of adjacent habitats to agricultural fields in the landscape have on crop damage.
In order to test how varying composition and configurations of adjacent habitats influence crop damage, I will construct high-resolution diet profiles of deer. High-resolution dietary analysis provides a mechanism to quantify the level of crop damage caused by a single deer. Plant DNA can be extracted from fecal pellets. Using next generation sequencing and metabarcoding techniques, plant material in the fecal pellets can be identified to the species level (Kartzinel et al. 2015). Metabarcoding utilizes short sections of DNA from a plant gene compared against reference libraries of this gene to pinpoint the sequence to a specific species. I will barcode the internal transcribed spacer 2 (ITS2) locus utilizing new universal primers (UniPlantF and UniPlantR), which has been suggested as the 'gold standard' for barcoding plants for herbivory studies (Moorhouse-Gann et al. 2018). This locus has a high level of taxonomic resolution as well as having a rather short amplicon length (160-320 bp) in order for it to mostly stay intact after passing thorough the digestive system of an ungulate. We will then use targeted amplification of this loci during PCR, each sample will be indexed using Nextera XT indexed adaptors and sequencing will occur on the Illumina MiSeq. Metabarcoding will allow me to identify the relative frequency of occurrence for each plant species (i.e. presence/absence of a certain plant species averaged across all fecal pellets) as well as the proportion of each species present in the fecal pellet utilizing relative read abundance (RRA) (Kartzinel et al. 2015).
Fecal pellets will be collected randomly throughout the predominant habitats in the landscape (agricultural fields, forested regions, and marshes). Time-based fecal pellet accumulation curves will be utilized to ensure an exhaustive search. The number of fecal pellets discovered will be graphed by the time spent collecting these fecal pellets, making sure that the number of new fecal pellets discovered asymptotes indicating overall fecal abundance. This indicates sampling effort is adequate to collect all fecal pellets present. Fecal pellets will be collected during 3 time points throughout the year in order to understand the seasonal variation in the importance of neighboring habitat. Collections will occur in March prior to planting, June during the growing season of the crops, and lastly September-October prior to harvest.
Fecal pellets will be visually inspected for freshness based on pellet color, consistency, and sheen (Brinkman et al. 2010). Fecal pellets will be ranked for freshness ranging from high quality (sticky with mucus, slight sheen on pellets, soft texture) to low quality (firm and dry but still with a sheen of mucus) (Kalb and Bowman 2018). We will not count fecal pellets that show signs of mold, appear broken or decayed, or lack a mucosal sheen on the exterior. A total of at least 60 pellets will be analyzed for DNA per sampling effort, resulting in at least 180 fecal pellets included in dietary analyses. A sample-based species accumulation curve will ensure the number of samples accurately reflects the plant species present in the diet of the deer. Fecal pellets will be randomly sampled from those collected in the field to be included in the diet analysis. We will utilize a fecal pellet index to quantify deer density at each farm sampled (following Forsyth 2005).
A GPS unit will provide an exact location of each fecal sample collected, which will provide the spatial data necessary to conduct landscape analysis. Spatial buffers will be created around the point locations of the 60 fecal pellets included in the dietary analysis. Buffer size will be the average distance traveled by deer in this region (Eyler 2001) multiplied by the food passage rate (Mautz and Petrides 1971), which will approximate the region in which foraging occurred for that fecal pellet. Within buffers, I will measure the landscape variables of patch size (overall area of each habitat), patch density (number of habitat patches), and patch connectivity (links or overlap between similar habitat patches).
In order to evaluate hypothesis 1, I will conduct a multiple regression to estimate the quantity of agricultural species present in the diet of the deer as a function of the season of sampling, deer density, and the calculated landscape metrics as predictor variables. This will allow us to directly access how well the composition and configuration of neighboring habitats predicts crop damage. Hypothesis 2 will be evaluated by Akaike information criterion (AIC) comparing models only using deer density to explain the percentage of agricultural species in their diet with models using landscape metrics. Conducting a multivariate multiple regression to predict the proportion of each species present in the diet of the deer as a function of deer density, season of sampling, and landscape metrics will assess hypothesis 3.
Objective 2. Evaluate the impact of dietary niche overlap of sympatric deer species on crop damage.
To investigate how dietary niche overlap of co-occurring deer species (white-tailed and sika) impacts crop damage, we will take advantage of the disjunct and overlapping spatial distribution of these two species in the Lower Eastern Shore of Maryland. Using many of the same methods as in the previous objective, we will identify the relative frequency of occurrence for each plant species as well as the proportion of each species present in the fecal pellet when the two deer species occur alone or together. As above, metabarcoding of the rbcL and ITS plant genes in fecal pellets will enable us to determine high-resolution diets of both species of deer. However, in addition to the methods described above, I will also sequence DNA of a “minibarcode” of the mitochondrial gene Cytochrome Oxidase 1 (CO1) used to differentiate animal species. I will use this minibarcode sequencing to discern the deer species from which the fecal pellets originated.
Sika deer and white-tailed deer on the Eastern Shore of Maryland occur in isolation and sympatrically. Fecal pellets will be collected throughout regions in which both species occur (Buck Range Farm; see letter of support from Donald Webster), and regions in which they occur in isolation (e.g. Blackwater National Wildlife Refuge for sika deer; see letter of support from Matt Whitbeck) (Somerset County for white-tailed deer; see letter of support from Sarah Hirsch). As in Objective 1, time-based fecal pellet accumulation curves, seasonal sampling, and pellet freshness quality assurance will be used in a rigorous and informative sampling protocol. A total of 30 fecal pellets will be analyzed for plant DNA for each of the 3 regions (only sika, both species, and only white-tailed) per time point in four time points, resulting in 360 fecal pellet samples for dietary analysis. Again, we will utilize a fecal pellet index to quantify deer density at each farm sampled.
Statistically, I will compare whether deer density or deer diet preferences are a stronger predictor of deer crop damage. In order to test Hypothesis 1, a one-way ANOVA will be conducted using total deer density as a response variable and deer species present in the region (whitetail deer, sika deer, or both) as a predictor. To evaluate Hypothesis 2, I will utilize a two-way ANOVA with the quantity of agricultural species present in the diet of the deer as a response variable and predictor variables of season of fecal pellets collection (pre-planting, growing season, pre-harvest, cover crop) and deer species present. This will allow me to test if the amount of agricultural species present in the diet of each deer varies depending upon if they occur in isolation or together. In order to evaluate hypothesis 3, I will calculate total crop damage inflicted by taking deer density at each site and multiplying it by the average percentage of agricultural species present in their diet. We will then conduct a one-way ANOVA using this total crop damage as the response variable and whether the region contains whitetail deer, sika deer, or both species co-occur as a predictor variable. Lastly, I will use non-metric dimensional scaling (NMDS), an ordination technique, to understand how the plant community present in each species’ diet varies and test hypothesis 4. I will use Analysis of Similarity (ANOSIM) to test for differences between each species’ diet and vector fitting to test how deer density affects diet composition.
* Objective 2 had to be modified due to restricted access to farms in regions which both species co-occur. I instead reallocated this effort to pursue a related question on the spatial distribution of crop damage at the patch and landscape-level. The methods for Objective 3 are outlined below.
Objective 3. Understand landscape dynamics that influence the spatial distribution of crop damage at the patch and landscape-level.
To investigate how landscape dynamics such as edge composition, field configuration, adjacent habitat connectivity and more influence the spatial distribution of crop damage, we will survey agricultural fields for crop damage. We will select 10 agricultural fields that have been planted with soybeans (Glycine max) as this is the predominant crop planted on the Eastern Shore of Maryland that experiences high levels of ungulate crop damage. Also, all farmers participating in this study plant soybeans whereas only two farmers plant corn (Zea mays). We will code a random number generator that will give us row and distance into the field for each sampling point. We intend to sample 10 points from each corner and edge of the field, for the interior we will sample an equal number of samples as collected for the edges and corners. This is due to the fact some fields will have more than 4 corners and edges due to the shape of the field and we intent to have equal sample sizes across all 3 categories. A corner will be defined as a 40 x 40-meter square from where two edges meet, an edge will be defined as an area 40 meters into the field that wasn’t already defined as a corner and the interior will be everything that was not defined as a corner or edge (Figure 1). At each sampling point we shall record a GPS location, the height of the 10 nearest plants as well as if there are any signs of ungulate herbivory on these plants. The same individual will sample all 10 fields and had a high level of experience deciphering ungulate herbivory from other types of herbivory. We intend to return to 30 of these sampling points just prior to harvest to survey for yield data. Yield data will be assessed by sampling 5 soybean plants at each point for the number of pods present on the plant, on one of the five plants the number of seeds present will be counted to calculate a seeds per pod average.
We will map the proportion of herbivory and average plant height utilizing the GPS data collected. Utilizing ArcGIS we plant to use Empirical Bayesian Kriging to map a smooth surface of proportion of crop damage and average plant heights across the entire field. Statistically, we plant to utilize spatial analyses such as Local Moran’s I and Getis-Ord Gi* to highlight any clustering or hot and cold spots in the data. We will then construct a hierarchical explanatory model using a generalized linear mixed effects model that incorporates patch-level, class-level and landscape-level metrics to explain spatial variability in crop damage present both within the field and between fields. Patch-level metrics that will be included in the model are distances to both the nearest and second nearest edge as well as the landcover class of these edges. Class-level and landscape-level metrics will be calculated from landcover classifications of buffer zones around each field. We will utilize the USDA's 2021 National Cropland Data Layer for landcover classification, landcover types utilized were corn, sorghum, soybean, CRP, forest, shrub, developed, wetland and open water. Utilizing ArcGIS we will clip this layer to a constructed 220 hectare buffer zone (average home-range of female white-tailed deer in this region (Eyler 2001)). Each farm's clipped landcover classification file will then be imported into R studio, utilizing the 'landscapemetrics' package both class-level and landscape-level metrics can be calculated for each field (Figure 2). Class-level metrics calculated for each field include core area index, patch density, connectivity, patch cohesion index and shape index for each landcover class present in the landscape (CRP, forest, etc.). Landscape-level metrics calculated for each field include core area index, relative mutual information criteria, Shannon’s diversity index, patch density, percentage of like adjacencies and contiguity index. These metrics were included to incorporate a wide range of different types of metrics (shape, diversity, complexity, aggregation, area and shape) in order to highlight what characteristics of the surround landscape influence crop damage. Allowing us to better understand how factors within the field as well as factors in the landscape surrounding the field influence how crop damage occurs across the field.
Initial trial sampling efforts were conducted in late August to early September 2020 at 2 farm sites (Barneck Farm and Anderson Farm) a total of 17 fecal pellets were collected. In March 2021 to May 2021 with the help of extension agents in Somerset County, we were able to connect with 5 other farmers experiencing heavy crop damage due to deer that were willing to allow us to use their farms as sites in our study. Starting in Late May/ early June 2021 we began collecting fecal pellets across 7 farmers on the Lower Eastern Shore of Maryland. We conducted another sampling effort in late August/ early September prior to the start of hunting season. A total of 62 fecal pellet samples were collected in June and 61 fecal pellet samples in late August/ early September across all 7 farms. A final sampling effort was conducted in late February/ early March of 2022 prior to any vegetation being present on the agricultural fields and 62 fecal pellet samples were collected across all 7 farms sampled. At the end of 2022 we have collected a total lc9d9of 185 fecal pellet samples. All 185 samples have been processed and DNA extractions have occurred. These samples have been sent off to the George Washington Genomics Core to get prepared for amplification and sequencing as well as quality control checks. Sequencing will be completed by March 2023 and computational analysis of the resulting data will begin in which relative read frequency for both agricultural species and non-agricultural species will be determined. This data will be incorporated with landscape metrics already calculated from Objective 3 in order to assess how adjacent habitats influence the proportion of agricultural species present in the diet of the deer
Due to the impasses relating to Objective 2, we decided to not only collect fecal pellets on the farm sites but also survey a select number of fields on each farm for crop damage in 2021. A total of 10 field were surveyed across the same 7 farms that we collect fecal pellets from. Across each field we sampled 120-180 randomly distributed points (depending upon the number of edges and corners the field had). We then utilized this data and Empirical Bayesian Kriging techniques is ArcGIS to be able to produce smoothed maps of crop damage and plant height of the entire field. We found that there is a high level of spatial variability present across the field, with certain regions experiencing higher levels of crop damage then others (Figure 1). In particular, the edges and corners experienced highest levels of herbivory (45-80% of plants experiencing herbivory). Crop damage was not limited to only the edges and corners but in some fields extended well into the interior (0-55% of plants experiencing herbivory) (Figure 2). In order to understand variability in crop damage within the field we constructed generalized linear mixed effects models that incorporated patch-level metrics which included distances to both the nearest and second nearest edges as well as the edges landcover composition. We found that the best model based off Akaike Information Criterion (AIC) included the distance to both edges as well as the landcover composition (Figure 3). The output of this model showed that herbivory decreased dramatically as distance to the first edge increased whereas it decreased more much less rapidly as distance to the second edge increased. It also showed the edge composition of the first edge impacted herbivory with forest and CRP edges had a more pronounced impact on herbivory than developed or agricultural edges (Figure 4).
However, there was a high level of spatial variability in crop damage between fields, highlighting the possibility that the surrounding landscape in which a field occurs could influence the level of crop damage on that field. In order to address how the surrounding landscape could influence crop damage we utilized a generalized linear mixed effect model framework to create hierarchical models that incorporated the patch-level model above (which we refer to as our field model) as well as landscape-level metrics for the entire field. We included core area index (landscape area metric), percent of like adjacencies (landscape connectivity), contiguity index (landscape shape) and relative mutual information (landscape complexity) as out metrics. We found the best model based of AIC included landscape area and landscape complexity (Figure 5). The output of this model showed that as core area present in the landscape and the relative mutual information increased so did the level of herbivory on the field (Figure 6). This highlights that as the landscape gets more fragmented and less interconnected with a greater amount of different land cover classes (more complex) the level of herbivory or crop damage increases. We also saw a shift in the effect of distance to second edge from just including the patch-level or "field model" with herbivory increasing as distance to second edge increased, the rest of the results remained consistent.
Lastly, we wanted to look at how including class-level metrics impacts the herbivory experienced by the field, we again utilized a generalized linear mixed effect model framework to create a hierarchical model that incorporates the patch-level model as well as class-level metrics. We primarily examined the 4 most prevalent land classes in our landscape (Forest, CRP, Soybean and Corn). We included an area metric and a connectivity metric for each land class in our model. We found the best model based of AIC included CRP mean area and CRP connectivity (Figure 7). This model showed that as CRP mean area increased so did herbivory and as CRP cohesion index or connectivity decreased so did herbivory experienced by the field (Figure 8). This highlighted as the amount of CRP on the landscape increases so does crop damage as well as when CRP present on the landscape is not well connected but fragmented crop damage increases as well. The patch-level metrics included in the model remained consistent with the output from the initial patch-level model.
Overall, when comparing all of the models created four models preformed equally as well based of having a ΔAIC of less than 2 (Figure 9). These include models that incorporate both landscape level and class-level metrics. Our results highlight some key aspects in regards to understanding crop damage, at the field or patch-level distance to edge is important in understanding the distribution of crop damage throughout the field. The landcover classification of the edge is also important, we found that forest edge and CRP edge influence crop damage experienced by the field greater than developed and agricultural edges. At the landscape-level, the more complex the landscape and the more core area present in the landscape both increase crop damage experienced by the field. Finally, at the land cover class level, when looking at the influence of certain land cover classes on crop damage including metrics beyond overall area of that land cover class is important. Often studies looking at utilizing landscape metrics or spatial pattern indices to explain crop damage neglect connectivity primarily focusing on overall area of each habitat in the landscape. We found including CRP mean area in combination with CRP connectivity better explained the crop damage experienced by the field, this also applied to looking at how forest land cover influenced crop damage.
Education & Outreach Activities and Participation Summary
Participation Summary:
This project will leverage the outreach program of a larger project funded by the USDA’s National Institute for Food and Agriculture (NIFA) that is researching the impacts of saltwater intrusion on agricultural fields in the Lower Eastern Shore of Maryland. The NIFA project has annual stakeholder meetings each spring in which farmers, extension agents, and other agricultural professionals participate, and I will have an opportunity to share my research results with this group. Unfortunately, due to COVID guidelines this stakeholders meeting did not occur in 2021, however in 2022 this meeting was held virtually with farmers, agricultural extension agents, and other researchers. I presented preliminary results during this meeting and talked with farmers about the issues they face in regards to crop damage. I plan on presenting completed results from my crop damage surveys in person to this stakeholder meeting in early 2023. In 2022, I also participated in an informal gathering of local farmers that one of the farmers that I collaborate with held. In which we discussed the issues of crop damage, what my research is trying to address, and issues they are facing in regards to the problem of crop damage. I also presented results from my crop damage field surveys to the 29th Annual Wildlife Society Conference in Spokane, Washington. At which I talked with researchers and land managers from both federal and state wildlife and agricultural agencies about the impacts of crop damage as well as utilizing the methodology I developed to survey crop damage more accurately.
Future outreach includes utilizing the University of Maryland agricultural extension agents to help disseminate my final results of this study. Sarah Hirsch in Somerset will assist me in this aim (see letter of support from S. Hirsch). As well as we are also currently plan on publishing the results of crop damage field surveys and the resulting models (Objective 3) in an academic journal in the coming year. My research will be shared at a state level in collaboration with the Maryland Department of Agriculture, which is responsible for developing Best Management Practices for agriculture. I plan on presenting my findings with Maryland Farm Bureau at their annual meeting, a leading voice of Maryland agriculture and predominant force in agriculture policy development in the state. Lastly, I will communicate my findings directly to local farmers during farm field days and extension meetings on the Eastern Shore of Maryland after the completion of my project in 2023. Enabling farmers experiencing extensive crop damage due to deer herbivory to make informed decisions on how to mitigate the intensity of crop damage occurring on their farm is the motivation for this research.