Introduction

The spotted lanternfly (Lycorma delicatula) is an invasive plant pest that has spread into several areas of North America. Temperature is an important factor that can influence its survival and potential distribution. Previous research showed that winter temperature can strongly affect the survival of overwintering eggs (Lee et al., 2014), while climatic suitability models have identified large areas where the species may potentially establish (Wakie et al., 2020).

A recent CBC News article discussed new research suggesting that spotted lanternfly eggs may be able to survive Canadian winters. This research was based on a study by Turbelin et al. (2025), who examined the cold tolerance of spotted lanternfly egg masses. Their experiments showed that some eggs were able to hatch after exposure to temperatures below -20°C, suggesting that the species may tolerate colder conditions than previously expected.

In this project, I used Dataset S5 from Turbelin et al. (2025), which contains hatch data from egg masses exposed to different temperatures and exposure durations. I used quantitative methods covered in BIOL 9915 to investigate how treatment temperature, exposure duration, and egg-mass class were related to hatch proportion. The analyses included simple linear regression, one-way ANOVA, two-way ANOVA, and ANCOVA.

Research questions

This project addressed four questions:

  1. Does treatment temperature affect hatch proportion after one hour of exposure?
  2. Does exposure duration affect hatch proportion at -15°C?
  3. Do exposure duration and egg-mass class interact to affect hatch proportion at -15°C?
  4. Does the relationship between treatment temperature and hatch proportion differ between single and multiple egg masses after one hour of exposure?

Data import and inspection

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(ggfortify)
slf <- read.csv("SLF_S5.csv")
glimpse(slf)
## Rows: 318
## Columns: 14
## $ Source          <chr> "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", …
## $ YR              <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ TreatmentDate   <chr> "18/01/2023", "18/01/2023", "18/01/2023", "18/01/2023"…
## $ TreatmentTemp   <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, …
## $ TreatmentTime   <chr> "control", "control", "control", "control", "control",…
## $ EggMassClass    <chr> "single", "multiple", "single", "single", "single", "m…
## $ EggMass         <chr> "281", "112", "267", "285", "264", "105", "282", "276"…
## $ EggMass2        <chr> "281", "112a", "267", "285", "264", "105a", "282", "27…
## $ HatchDate       <chr> "4/18/2023", "4/6/2023", "4/18/2023", "4/22/2023", "4/…
## $ TotalEggs       <chr> "0", "4", "15", "21", "24", "27", "29", "30", "31", "3…
## $ TotalHatch      <int> 0, 3, 1, 8, 0, 26, 4, 0, 19, 30, 27, 2, 0, 6, 26, 34, …
## $ PercentageHatch <dbl> NA, 0.75, 0.07, 0.38, 0.00, 0.96, 0.14, 0.00, 0.61, 0.…
## $ HatchStartDate  <chr> NA, "4/27/2023", NA, NA, NA, "4/20/2023", NA, NA, "4/2…
## $ HatchEndDate    <chr> NA, "4/28/2023", NA, NA, NA, "4/21/2023", NA, NA, "4/2…
head(slf)
##   Source YR TreatmentDate TreatmentTemp TreatmentTime EggMassClass EggMass
## 1     PA  2    18/01/2023             6       control       single     281
## 2     PA  2    18/01/2023             6       control     multiple     112
## 3     PA  2    18/01/2023             6       control       single     267
## 4     PA  2    18/01/2023             6       control       single     285
## 5     PA  2    18/01/2023             6       control       single     264
## 6     PA  2    18/01/2023             6       control     multiple     105
##   EggMass2 HatchDate TotalEggs TotalHatch PercentageHatch HatchStartDate
## 1      281 4/18/2023         0          0              NA           <NA>
## 2     112a  4/6/2023         4          3            0.75      4/27/2023
## 3      267 4/18/2023        15          1            0.07           <NA>
## 4      285 4/22/2023        21          8            0.38           <NA>
## 5      264 4/18/2023        24          0            0.00           <NA>
## 6     105a 3/31/2023        27         26            0.96      4/20/2023
##   HatchEndDate
## 1         <NA>
## 2    4/28/2023
## 3         <NA>
## 4         <NA>
## 5         <NA>
## 6    4/21/2023
table(slf$TreatmentTemp)
## 
##   -25 -22.5   -20   -15   -10    -5     0     6 
##    30    20    50    53    58    54    29    24
table(slf$TreatmentTime)
## 
##       1      12     240     360 control 
##     103      70      56      65      24
table(slf$EggMassClass)
## 
## multiple   single 
##       79      239
table(slf$YR)
## 
##   1   2 
##  35 283

Data cleaning

filter(slf, TotalEggs == "?")
##   Source YR TreatmentDate TreatmentTemp TreatmentTime EggMassClass EggMass
## 1     PA  2    18/01/2023             6       control       single     265
##   EggMass2 HatchDate TotalEggs TotalHatch PercentageHatch HatchStartDate
## 1      265 4/11/2023         ?          0               0       5/5/2023
##   HatchEndDate
## 1         <NA>
filter(slf, is.na(PercentageHatch))
##   Source YR TreatmentDate TreatmentTemp TreatmentTime EggMassClass EggMass
## 1     PA  2    18/01/2023             6       control       single     281
## 2     PA  2    24/01/2023            -5             1       single     110
## 3     PA  2    25/01/2023           -10            12       single      92
##   EggMass2 HatchDate TotalEggs TotalHatch PercentageHatch HatchStartDate
## 1      281 4/18/2023         0          0              NA           <NA>
## 2      110 4/22/2023         0          0              NA           <NA>
## 3       92 4/19/2023         0          0              NA           <NA>
##   HatchEndDate
## 1         <NA>
## 2         <NA>
## 3         <NA>

One observation had an unknown total egg count, recorded as “?”. Three observations had zero total eggs and therefore had missing hatch proportions. These three observations were excluded from analyses using hatch proportion.

slf_clean <- filter(slf, !is.na(PercentageHatch))

glimpse(slf_clean)
## Rows: 315
## Columns: 14
## $ Source          <chr> "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", …
## $ YR              <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ TreatmentDate   <chr> "18/01/2023", "18/01/2023", "18/01/2023", "18/01/2023"…
## $ TreatmentTemp   <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, …
## $ TreatmentTime   <chr> "control", "control", "control", "control", "control",…
## $ EggMassClass    <chr> "multiple", "single", "single", "single", "multiple", …
## $ EggMass         <chr> "112", "267", "285", "264", "105", "282", "276", "112"…
## $ EggMass2        <chr> "112a", "267", "285", "264", "105a", "282", "276", "11…
## $ HatchDate       <chr> "4/6/2023", "4/18/2023", "4/22/2023", "4/18/2023", "3/…
## $ TotalEggs       <chr> "4", "15", "21", "24", "27", "29", "30", "31", "33", "…
## $ TotalHatch      <int> 3, 1, 8, 0, 26, 4, 0, 19, 30, 27, 2, 0, 6, 26, 34, 45,…
## $ PercentageHatch <dbl> 0.75, 0.07, 0.38, 0.00, 0.96, 0.14, 0.00, 0.61, 0.91, …
## $ HatchStartDate  <chr> "4/27/2023", NA, NA, NA, "4/20/2023", NA, NA, "4/27/20…
## $ HatchEndDate    <chr> "4/28/2023", NA, NA, NA, "4/21/2023", NA, NA, "4/28/20…
slf_summary_clean <- slf_clean %>%
  group_by(TreatmentTemp, TreatmentTime) %>%
  summarise(mean_hatch = mean(PercentageHatch))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by TreatmentTemp and TreatmentTime.
## ℹ Output is grouped by TreatmentTemp.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(TreatmentTemp, TreatmentTime))` for per-operation
##   grouping (`?dplyr::dplyr_by`) instead.
slf_summary_clean
## # A tibble: 23 × 3
## # Groups:   TreatmentTemp [8]
##    TreatmentTemp TreatmentTime mean_hatch
##            <dbl> <chr>              <dbl>
##  1         -25   1                0.0558 
##  2         -25   12               0      
##  3         -22.5 1                0      
##  4         -22.5 12               0      
##  5         -20   1                0.135  
##  6         -20   12               0.00273
##  7         -20   240              0.073  
##  8         -20   360              0      
##  9         -15   1                0.352  
## 10         -15   12               0.345  
## # ℹ 13 more rows
ggplot(slf_summary_clean,
       aes(x = TreatmentTemp,
           y = mean_hatch,
           colour = TreatmentTime,
           group = TreatmentTime)) +
  geom_point() +
  geom_line() +
  xlab("Treatment temperature (°C)") +
  ylab("Mean hatch proportion") +
  theme_bw()

Mean hatch proportion was generally lower at colder temperatures, although the pattern differed among exposure durations. # Simple Linear Regression ## Effect of Temperature on Hatch Success After One Hour of Exposure

slf_1h <- filter(slf_clean, TreatmentTime == "1")

glimpse(slf_1h)
## Rows: 102
## Columns: 14
## $ Source          <chr> "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", …
## $ YR              <int> 2, 2, 1, 2, 2, 1, 2, 2, 2, 1, 1, 1, 1, 2, 1, 2, 2, 2, …
## $ TreatmentDate   <chr> "24/01/2023", "24/01/2023", "1-Mar-22", "24/01/2023", …
## $ TreatmentTemp   <dbl> -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5, -5…
## $ TreatmentTime   <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1",…
## $ EggMassClass    <chr> "single", "multiple", "single", "single", "multiple", …
## $ EggMass         <chr> "204", "153", "102_yr1", "62", "153", "79_yr1", "91", …
## $ EggMass2        <chr> "204", "153a", "102_yr1", "62", "153b", "79_yr1", "91"…
## $ HatchDate       <chr> "4/22/2023", "4/6/2023", "3/28/2022", "4/18/2023", "4/…
## $ TotalEggs       <chr> "16", "18", "25", "26", "31", "32", "33", "35", "35", …
## $ TotalHatch      <int> 0, 0, 0, 14, 29, 18, 17, 0, 5, 28, 0, 37, 33, 25, 33, …
## $ PercentageHatch <dbl> 0.00, 0.00, 0.00, 0.54, 0.94, 0.56, 0.52, 0.00, 0.14, …
## $ HatchStartDate  <chr> NA, "4/30/2023", NA, "10-05-2023", "4/30/2023", "18-Ap…
## $ HatchEndDate    <chr> NA, "5/1/2023", NA, "12-05-2023", "5/1/2023", "19-Apr-…
ggplot(slf_1h,
       aes(x = TreatmentTemp,
           y = PercentageHatch)) +
  geom_point() +
  xlab("Treatment temperature (°C)") +
  ylab("Hatch proportion") +
  theme_bw()

### Linear regression model

model_1h <- lm(PercentageHatch ~ TreatmentTemp,
               data = slf_1h)

Checking model assumptions

autoplot(model_1h, smooth.colour = NA)
## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggfortify package.
##   Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 102 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 102 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 102 rows containing missing values or values outside the scale range
## (`geom_line()`).

### Model results

anova(model_1h)
## Analysis of Variance Table
## 
## Response: PercentageHatch
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## TreatmentTemp   1 2.4875 2.48749  31.286 1.944e-07 ***
## Residuals     100 7.9507 0.07951                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_1h)
## 
## Call:
## lm(formula = PercentageHatch ~ TreatmentTemp, data = slf_1h)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.48681 -0.15307 -0.04182  0.17413  0.93818 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.598058   0.068715   8.703 6.81e-14 ***
## TreatmentTemp 0.022250   0.003978   5.593 1.94e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.282 on 100 degrees of freedom
## Multiple R-squared:  0.2383, Adjusted R-squared:  0.2307 
## F-statistic: 31.29 on 1 and 100 DF,  p-value: 1.944e-07

Interpretation

A simple linear regression was used to examine the relationship between treatment temperature and hatch proportion after one hour of exposure. Hatch proportion increased with increasing temperature (slope = 0.022), and the relationship was statistically significant (p < 0.001). However, inspection of the diagnostic plots indicated that some linear-model assumptions, particularly normality of the residuals, were not fully satisfied. Therefore, the regression results should be interpreted cautiously. ### Fitted regression

ggplot(slf_1h,
       aes(x = TreatmentTemp,
           y = PercentageHatch)) +
  geom_point() +
  geom_smooth(method = "lm") +
  xlab("Treatment temperature (°C)") +
  ylab("Hatch proportion") +
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

One-way ANOVA

slf_m15 <- filter(slf_clean, TreatmentTemp == -15)

table(slf_m15$TreatmentTime)
## 
##   1  12 240 360 
##  19  11  12  11
ggplot(slf_m15,
       aes(x = TreatmentTime,
           y = PercentageHatch)) +
  geom_boxplot() +
  xlab("Exposure duration (hours)") +
  ylab("Hatch proportion") +
  theme_bw()

### Model

model_m15 <- lm(PercentageHatch ~ TreatmentTime,
                data = slf_m15)

Checking model assumptions

autoplot(model_m15, smooth.colour = NA)
## Warning: Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).

### Model results

anova(model_m15)
## Analysis of Variance Table
## 
## Response: PercentageHatch
##               Df Sum Sq Mean Sq F value   Pr(>F)   
## TreatmentTime  3 1.1667 0.38888  5.7197 0.001943 **
## Residuals     49 3.3316 0.06799                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model_m15)
## 
## Call:
## lm(formula = PercentageHatch ~ TreatmentTime, data = slf_m15)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35211 -0.09083 -0.01455  0.09789  0.65545 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.35211    0.05982   5.886  3.5e-07 ***
## TreatmentTime12  -0.00756    0.09879  -0.077  0.93931    
## TreatmentTime240 -0.26127    0.09615  -2.717  0.00907 ** 
## TreatmentTime360 -0.33756    0.09879  -3.417  0.00128 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2608 on 49 degrees of freedom
## Multiple R-squared:  0.2594, Adjusted R-squared:  0.214 
## F-statistic:  5.72 on 3 and 49 DF,  p-value: 0.001943

Group means

slf_m15_means <- slf_m15 %>%
  group_by(TreatmentTime) %>%
  summarise(mean_hatch = mean(PercentageHatch))

slf_m15_means
## # A tibble: 4 × 2
##   TreatmentTime mean_hatch
##   <chr>              <dbl>
## 1 1                 0.352 
## 2 12                0.345 
## 3 240               0.0908
## 4 360               0.0145

Interpretation

At -15°C, hatch proportion differed significantly among exposure durations (F = 5.72, p = 0.0019). Using 1 h as the reference group, hatch proportion after 12 h was not significantly different from 1 h (p = 0.939), while hatch proportion was significantly lower after 240 h (p = 0.009) and 360 h (p = 0.001). Diagnostic plots showed some departures from normality and constant variance, so the results were interpreted cautiously. ### Final plot

ggplot(slf_m15_means,
       aes(x = TreatmentTime,
           y = mean_hatch)) +
  geom_point() +
  xlab("Exposure duration (hours)") +
  ylab("Mean hatch proportion") +
  theme_bw()

# Two-way ANOVA

table(slf_m15$TreatmentTime,
      slf_m15$EggMassClass)
##      
##       multiple single
##   1          4     15
##   12         2      9
##   240        4      8
##   360        2      9
slf_m15_two_summary <- slf_m15 %>%
  group_by(TreatmentTime, EggMassClass) %>%
  summarise(mean_hatch = mean(PercentageHatch))
## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by TreatmentTime and EggMassClass.
## ℹ Output is grouped by TreatmentTime.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(TreatmentTime, EggMassClass))` for per-operation
##   grouping (`?dplyr::dplyr_by`) instead.
ggplot(slf_m15_two_summary,
       aes(x = TreatmentTime,
           y = mean_hatch,
           colour = EggMassClass,
           group = EggMassClass)) +
  geom_point() +
  geom_line() +
  xlab("Exposure duration (hours)") +
  ylab("Mean hatch proportion") +
  theme_bw()

### Model

model_m15_two <- lm(PercentageHatch ~ TreatmentTime * EggMassClass, data = slf_m15)

Checking model assumptions

autoplot(model_m15_two, smooth.colour = NA)
## Warning: Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).

drop1(model_m15_two, ~ ., test = "F")
## Single term deletions
## 
## Model:
## PercentageHatch ~ TreatmentTime * EggMassClass
##                            Df Sum of Sq    RSS     AIC F value Pr(>F)
## <none>                                  3.0397 -135.50               
## TreatmentTime               3  0.297892 3.3376 -136.55  1.4700 0.2354
## EggMassClass                1  0.036267 3.0760 -136.87  0.5369 0.4675
## TreatmentTime:EggMassClass  3  0.240986 3.2807 -137.46  1.1892 0.3246

Model without interaction

model_m15_two_add <- lm(PercentageHatch ~ TreatmentTime + EggMassClass,
                        data = slf_m15)

drop1(model_m15_two_add, ~ ., test = "F")
## Single term deletions
## 
## Model:
## PercentageHatch ~ TreatmentTime + EggMassClass
##               Df Sum of Sq    RSS     AIC F value   Pr(>F)   
## <none>                     3.2807 -137.46                    
## TreatmentTime  3   1.18902 4.4697 -127.07  5.7988 0.001822 **
## EggMassClass   1   0.05084 3.3316 -138.64  0.7438 0.392733   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model results

summary(model_m15_two_add)
## 
## Call:
## lm(formula = PercentageHatch ~ TreatmentTime + EggMassClass, 
##     data = slf_m15)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41112 -0.10637 -0.05637  0.10888  0.66905 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.411122   0.090994   4.518 4.07e-05 ***
## TreatmentTime12    -0.005414   0.099080  -0.055  0.95665    
## TreatmentTime240   -0.270452   0.096986  -2.789  0.00757 ** 
## TreatmentTime360   -0.335414   0.099080  -3.385  0.00143 ** 
## EggMassClasssingle -0.074755   0.086678  -0.862  0.39273    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2614 on 48 degrees of freedom
## Multiple R-squared:  0.2707, Adjusted R-squared:  0.2099 
## F-statistic: 4.453 on 4 and 48 DF,  p-value: 0.003847

Checking final model assumptions

autoplot(model_m15_two_add, smooth.colour = NA)
## Warning: Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 53 rows containing missing values or values outside the scale range
## (`geom_line()`).

### Interpretation At -15°C, there was no significant interaction between exposure duration and egg-mass class (F = 1.19, p = 0.325). After removing the interaction, exposure duration had a significant effect on hatch proportion (F = 5.80, p = 0.0018), while egg-mass class did not have a significant effect (F = 0.74, p = 0.393). Using 1 h as the reference group, hatch proportion after 12 h was not significantly different from 1 h (p = 0.957), while hatch proportion was significantly lower after 240 h (p = 0.0076) and 360 h (p = 0.0014). Diagnostic plots showed some departures from normality and constant variance, so the results were interpreted cautiously.

ANCOVA

table(slf_1h$TreatmentTemp, slf_1h$EggMassClass)
##        
##         multiple single
##   -25          4     15
##   -22.5        0     10
##   -20          0     17
##   -15          4     15
##   -10          6     14
##   -5           2     15
ggplot(slf_1h,
       aes(x = TreatmentTemp,
           y = PercentageHatch,
           colour = EggMassClass)) +
  geom_point() +
  xlab("Treatment temperature (°C)") +
  ylab("Hatch proportion") +
  theme_bw()

model_ancova <- lm(PercentageHatch ~ TreatmentTemp * EggMassClass,
                   data = slf_1h)

Checking model assumptions

autoplot(model_ancova, smooth.colour = NA)
## Warning: Removed 102 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 102 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 102 rows containing missing values or values outside the scale range
## (`geom_line()`).

drop1(model_ancova, ~ ., test = "F")
## Single term deletions
## 
## Model:
## PercentageHatch ~ TreatmentTemp * EggMassClass
##                            Df Sum of Sq    RSS     AIC F value Pr(>F)
## <none>                                  7.7718 -254.60               
## TreatmentTemp               1  0.054965 7.8268 -255.88  0.6931 0.4071
## EggMassClass                1  0.090222 7.8620 -255.42  1.1377 0.2888
## TreatmentTemp:EggMassClass  1  0.158817 7.9306 -254.53  2.0026 0.1602

Model without interaction

model_ancova_add <- lm(PercentageHatch ~ TreatmentTemp + EggMassClass,
                       data = slf_1h)

drop1(model_ancova_add, ~ ., test = "F")
## Single term deletions
## 
## Model:
## PercentageHatch ~ TreatmentTemp + EggMassClass
##               Df Sum of Sq     RSS     AIC F value    Pr(>F)    
## <none>                      7.9306 -254.53                      
## TreatmentTemp  1   2.43041 10.3610 -229.27 30.3395 2.871e-07 ***
## EggMassClass   1   0.02008  7.9507 -256.27  0.2506    0.6177    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Checking final model assumptions

autoplot(model_ancova_add, smooth.colour = NA)
## Warning: Removed 102 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 102 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 102 rows containing missing values or values outside the scale range
## (`geom_line()`).

Interpretation

After one hour of exposure, there was no significant interaction between treatment temperature and egg-mass class (F = 2.00, p = 0.160). After removing the interaction, treatment temperature had a significant effect on hatch proportion (F = 30.34, p < 0.001), while egg-mass class did not have a significant effect (F = 0.25, p = 0.618). This indicates that hatch proportion changed with temperature, but there was no evidence that the relationship differed between single and multiple egg masses. Diagnostic plots showed some departures from normality and constant variance, so the results were interpreted cautiously.

Discussion

The analyses showed that temperature and exposure duration were important factors affecting spotted lanternfly egg hatch. After one hour of exposure, hatch proportion increased as treatment temperature increased. This agrees with earlier research showing that winter temperature can strongly affect the survival of overwintering spotted lanternfly eggs (Lee et al., 2014). It is also consistent with the cold-tolerance results reported by Turbelin et al. (2025).

Exposure duration was also important. At -15°C, hatch proportion was much lower after 240 and 360 hours of exposure than after one hour, while 12 hours was not significantly different from one hour. This suggests that the effect of cold depends not only on how low the temperature is, but also on how long the eggs are exposed to it. This pattern is consistent with the experimental work of Turbelin et al. (2025), which examined both temperature and duration when evaluating cold tolerance.

Egg-mass class had a different pattern. In both the two-way ANOVA and ANCOVA, there was no evidence that single and multiple egg masses responded differently under the conditions examined here. Therefore, temperature and exposure duration appeared to be more important predictors of hatch proportion in this dataset than egg-mass class.

Temperature is also important beyond the egg stage. Kreitman et al. (2021) showed that temperature affects development and survival of spotted lanternfly nymphs. At a larger geographic scale, Wakie et al. (2020) identified temperature-related variables as important predictors of potential establishment. More recently, Zhao et al. (2024) identified mean temperature of the driest quarter as an important variable affecting predicted global suitability and reported a possible northward expansion of suitable areas under future climate scenarios.

Together, these studies show why cold tolerance is important when considering the potential spread of spotted lanternfly into colder regions. The results of the present analysis support the idea that cold temperatures reduce hatch success, particularly when exposure is prolonged. However, the results should not be interpreted as direct evidence that spotted lanternfly will establish in Canada because establishment also depends on other climatic and biological factors.

There were also limitations in the statistical analyses. Some diagnostic plots showed departures from normality and constant variance, and the experimental design was unbalanced for some combinations of temperature, exposure duration, and egg-mass class. Therefore, the results from the linear models should be interpreted cautiously.

Conclusion

This project used Dataset S5 from Turbelin et al. (2025) to examine how cold exposure affected spotted lanternfly egg hatch. Treatment temperature and exposure duration were the most important factors in the analyses, while egg-mass class showed no significant effect under the conditions examined.

Overall, colder temperatures and longer exposure reduced hatch proportion. These results support previous research showing that temperature is an important factor affecting spotted lanternfly survival and potential distribution. However, because some statistical assumptions were not fully met and the experimental design was unbalanced, the results should be interpreted cautiously.

The findings show how quantitative methods can be used to connect experimental data with broader biological questions about cold tolerance and the potential spread of an invasive species.

References

Greco, F. (2026, June 29). New research suggests invasive spotted lanternfly may be able to survive Canadian winters. CBC News.

Kreitman, D., Keena, M. A., Nielsen, A. L., & Hamilton, G. (2021). Effects of temperature on development and survival of nymphal Lycorma delicatula (Hemiptera: Fulgoridae). Environmental Entomology, 50(1), 183–191. doi:10.1093/ee/nvaa155

Lee, Y. S., Jang, M. J., Kim, J. Y., & Kim, J. R. (2014). The effect of winter temperature on the survival of lantern fly, Lycorma delicatula (Hemiptera: Fulgoridae) eggs. Korean Journal of Applied Entomology, 53(3), 311–315. doi:10.5656/KSAE.2014.07.0.034

Turbelin, A. J., Sinclair, B. J., Rost, J., & Roe, A. D. (2025). Cold tolerance strategy and lower temperature thresholds of Lycorma delicatula (Hemiptera: Fulgoridae) egg masses. Environmental Entomology, 54(4), 890–898. doi:10.1093/ee/nvaf007

Wakie, T. T., Neven, L. G., Yee, W. L., & Lu, Z. (2020). The establishment risk of Lycorma delicatula (Hemiptera: Fulgoridae) in the United States and globally. Journal of Economic Entomology, 113(1), 306–314. doi:10.1093/jee/toz259

Zhao, Z., Yang, L., & Chen, X. (2024). Globally suitable areas for Lycorma delicatula based on an optimized Maxent model. Ecology and Evolution, 14(9), e70252. doi:10.1002/ece3.70252