Species 1: Spodoptera frugiperda (Fall Armyworm Moth)

Species 2: Helicoverpa armigera (Cotton Bollworm Moth)

Song, Y. et al. (2023) Interspecific Competition between Invasive ISpodoptera frugiperda/I and Indigenous IHelicoverpa armigera/I in Maize Fields of China. Agronomy (Basel). [Online] 13 (3).

(Note loading species data from gbif may take a few minutes… or a while)

1. Run linear models to predict the present-day distribution of species 1 using climate variables and use them to present a map of its current distribution. Which set of climatic variables best explain the current distribution of the species?

## Best model AIC: 3266.345

To generate this plot, a function runs a linear model to predict present-day distributions of S. frugiperda. Using imputed climate data, the function fits the model for each climate variable (19 of them; provided by WorldClim) and selects the best fitting model using its Akaike Information Criterion value. It then plots the probability of occurrence using the best climate model and a prediction of the S. frugiperda distribution by generating presence/absence data using randomly generated points as a null model. This interaction of the model shows ‘bio12’ as the best predictor of S. frugiperda occurrence. Bio12 is annual precipitation and shows that precipitation near the equator is a good predictor of S. frugiperda presence.

(due to the randomly generated variable in the model, the ‘best model’ may be different when ran multiple times. If used in study this should be repeated, and a statically significant result should be found.)

2. Run linear models to predict the present-day distribution of species 2 using climate variables and use them to present a map of its current distribution. Which set of climatic variables best explain the current distribution of the species?

## Best model AIC: 3363.926

To generate this plot, the same function for Q.1 is run however this time with data for H.armigera imputed. This results in bio7 being selected as the best fitting model. This is the annual temperature range, given by the maximum temperature in the warmest month (bio5) - the minimum temperature in the coldest month (bio6). When used to work out presence/absence the model shows H.armigera is distributed across Europe, South-East Asia, coastal Australia, coastal Antarctica (unlikely), South America, and across Sub-Saharan Africa. 

(due to the randomly generated variable in the model, the ‘best model’ may be different when ran multiple times. If used in study this should be repeated, and a statically significant result should be found.)

3. Plot the overlap in distribution of the two species. Devise and calculate a metric for the degree of overlap between their ranges, explaining how you calculated it.

To generate this plot, I use a function that takes coordinate data from both species and plots. Red shows the distribution of S.frugiperda (n=7986); blue shows the distribution of H.armigera (n=29,577); green (n=634) shows overlaps. To work out overlaps, I used mean flight distance data for (female) H.armigera of 41.25 km reported by Wu et al. (2006). As a threshold radius. The radius around each point (chosen due to smaller dataset) is worked out and each point of the other species that is found within it is recorded as an overlap. This automatically ignores overlap in ranges across oceans, as ocean data is removed by a clean function run on all data at the start of this code. This is used to construct an overlap coordinate data frame that is fed into the same plotting function as the individual species data.

Wu, H. et al. (2006) Flight Potential of Pink Bollworm, Pectinophora gossypiella Saunders (Lepidoptera: Gelechiidae). Environmental entomology. [Online] 35 (4), 887–893.

4. Use a linear model to test whether the distribution of species 1 at the present time depends on the distribution of species 2, while also taking account of the effects of climatic variables.

## Best model AIC: 596.9147

## 
## Call:
## glm(formula = modelFormula2, data = envtrain2)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.703e-01  1.338e-02 -50.109  < 2e-16 ***
## bio4         3.326e-04  1.191e-05  27.936  < 2e-16 ***
## bio15       -9.169e-04  1.390e-04  -6.594 4.38e-11 ***
## bio1         6.550e-02  6.351e-04 103.128  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1403945)
## 
##     Null deviance: 4797.3  on 19893  degrees of freedom
## Residual deviance: 2792.4  on 19890  degrees of freedom
##   (1 observation deleted due to missingness)
## AIC: 17405
## 
## Number of Fisher Scoring iterations: 2

The same model function as Q1. and Q2. is run however this time the effects of S.frugiperda is compared to H.armigera with the climate variables best fitting H.armigera (bio7) and S.frugiperda (bio12). Additionally the climate variables that best fit the overlaps are considered. This is reported as bio11, which is the mean temperature in the coldest quarter. This suggests that the possible interaction between S.frugiperda and H.armigera is dependent on a minimum temperature level — this possibly represents where temperature is hot enough to maximise precipitation. I think this is backed up by ‘quarters’ being the measurement of time (as opposed to year for bio12 and month for bio7). This may be a better scale to measure a rainy season that happens during temperatures above a minimum value.

5. Predict the future distribution of each species using CMIP6 data for future climate and predict how the degree of overlap in ranges change will change in that time. Do you expect the two species to co-occur more often or less often than at the present?

Future prediction: S. frugiperda

## Best model AIC: 3196.151

## 
## Call:  glm(formula = modelFormula, family = binomial(link = "logit"), 
##     data = envtrain)
## 
## Coefficients:
## (Intercept)        bio12  
##     0.30446      0.00238  
## 
## Degrees of Freedom: 8099 Total (i.e. Null);  8098 Residual
##   (475 observations deleted due to missingness)
## Null Deviance:       3677 
## Residual Deviance: 3192  AIC: 3196

Future prediction: H. armigera

## Best model AIC: 3058.553

## 
## Call:  glm(formula = modelFormula, family = binomial(link = "logit"), 
##     data = envtrain)
## 
## Coefficients:
## (Intercept)        bio15  
##     4.96518     -0.04182  
## 
## Degrees of Freedom: 9946 Total (i.e. Null);  9945 Residual
##   (2373 observations deleted due to missingness)
## Null Deviance:       3894 
## Residual Deviance: 3055  AIC: 3059

Future Predictions: Overlap

## Best model AIC: 548.707

## 
## Call:  glm(formula = modelFormula, family = binomial(link = "logit"), 
##     data = envtrain)
## 
## Coefficients:
## (Intercept)         bio4  
##    1.039741    -0.004983  
## 
## Degrees of Freedom: 664 Total (i.e. Null);  663 Residual
##   (25 observations deleted due to missingness)
## Null Deviance:       790.1 
## Residual Deviance: 544.7     AIC: 548.7

These plots are generated using the same climate prediction model function however this time the input data for climate is changed to future projections between 2061-2080 provided by the CMIP6 database. This is given for  S.frugiperda that retains bio12 as its climate variable, H.armigera is driven to be more precipitation dependent given by its new climate value of bio17 (Precipitation of Driest Quarter). This however, does translate to a blank presence/absence plot possibly suggesting extinction at current rates of climate change. Bio 6 is retained in the overlapping data, however a slight shift towards the equator from both sides is observed. This may be due to dryer spells above and below the interaction range pushing both species more equatorial — suggesting the interaction between S.frugiperda and H.armigera is precipitation dependent.