Methods, Results, Figures

Methods

Sampling

Question: Should I use the terminology asexual and sexual or unisexual and bisexual? I realized I’ve been using asexual and sexual without thinking too hard about it. Just because there’s a unisexual population, it doesn’t mean they’re not at least facultatively parthenogenic, right?

A total of six Phasmid species were sampled from across New Zealand (Table). We selected species that contain multiple asexual and sexual populations. Need sampling details.

Environmental data

Marginality

Niche marginality is the magnitude which the center of species’ Hutchinsonian niche differs from available environmental conditions (Hirzel et al. 2002). We estimated the marginality of each species-reproductive mode combination using the adehabitatHS v0.3.14 package in R v3.6.0 (adehabitat CITE, R CITE). Statistical significance was assessed using a permutation-based approach, implemented in the R package adehabitatHS v0.3.14 (Elaborate?). We obtained 19 bioclimatic variables from the CHELSA database v1.2 to summarize current climatic conditions (Karger et al. 2017). These variables are derived from the ERA-Interim climatic analysis and capture broad-scale temperature and precipitation patterns. The variables were scaled at 1 km resolution. We sampled all environment from within a 50 km buffer placed around a minimum convex hull of each species’ sampling range to represent the available background environment.

Stability

To capture historical climate dynamics, we obtained climate data from seven additional time periods: the late-Holocene (4.2-0.3 kya), the mid-Holocene (8.326-4.2 kya), the early-Holocene (11.7-8.326 kya), the Younger Dryas Stadial (12.9-11.7 kya), the Bølling-Allerød (14.7-12.9 kya), the Heinrich Stadial 1 (17.0-14.7 kya), and the Last Glacial Maximum (~21 kya). All variables were downloaded from the Paleoclim database (Brown cite). All variables except the Last Glacial Maximum (LGM) were derived from Fordham et al. (2017). The LGM data were derived from the CHELSA v1.2 database, which uses the PMIP3 circulation model (Karger et al. 2017). All variables were scaled to 5 km resolution, and current climate was downscaled to match the resolution of the other variables. We estimated climate stability since the LGM using the average annual temperature and average annual precipitation variables. We used the climateStability v0.1.1 R package (Owens and Guralnick 2019) to calculate stability across time periods. Briefly, the method calculates the stability of a climate variable as inverse of the mean standard deviation between time slices over time elapsed (Owens and Guralnick 2019). In addition to temperature and precipitation stability, we estimated overall climate stability as the product of the two stability values. We then extracted stability values for each species-reproductive mode combination, where a stability value was recovered for each individual locality. To assess statistical differences in occupied stability between reproductive modes, we conducted Welch’s t-tests (CITE) for each species and stability metric.

Seasonality

Results

Marginality

Sexual individuals occupied more marginal habitat than asexual individuals for all species
p-values to come

Stability

Asexual individuals occupy siginficantly more stable climate space than sexual individuals in three species- Asteliaphasma jucundum, Clitarchus hookeri, and Tectarchus ovobessus. In A. jucundum and C. hookeri, this is in all three metrics. In T. ovobessus, this is for temperature.
Sexual individuals occupy significantly more stable climate space than asexual individuals in Argosarchus horridus according to all three metrics
See table for significant results

Seasonality

Figures and Tables

Sampling map

Caption: Sampling localities for the six species included in this study. Open circles signify asexual individuals, while closed circles signify sexual individuals. Sample sizes are as follows: Argosarchus horridus: 96 asexual, 213 sexual; Asteliaphasma jucundum: 21 asexual, 159 sexual; Clitarchus hookeri: 147 asexual, 351 sexual; Niveaphasma annulata: 111 asexual, 114 sexual; Tectarchus huttoni: 33 asexual, 45 sexual; Tectarchus ovobessus: 60 asexual, 27 sexual.

locs <-
  readxl::read_excel(here("data/all species New_6-14-19.xlsx")) %>%
  rename(reproductive_mode = `Reproductive mode`) %>%
  filter(
    genus != "acanthoxyla",
    genus != "micrarchus",
    genus != "spinotectarchus",
    genus != "Tepakiphasma",
    species != "nov. sp.",
    species != "rakauwhakanekeneke",
    species != "tepaki",
    species != "salebrosus"
  ) %>%
  mutate(genus = str_to_sentence(genus)) %>%
  mutate(genus_sp = str_c(genus, species, sep = " ") %>% str_to_sentence()) %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)

nz_basemap <-
  st_read(here("data", "NewZealand_Boundary.kml"), quiet = TRUE) %>% st_crop(
    xmin = 165,
    xmax = 179,
    ymin = -48,
    ymax = -34
  ) %>%
  rmapshaper::ms_simplify()

full_map <- ggplot() +
  geom_sf(data = nz_basemap) +
  facet_wrap(~ genus_sp) +
  geom_sf(data = locs, aes(shape = reproductive_mode), show.legend = "point") +
  expand_limits(x = 180) +
  scale_shape_manual(values = c(1, 19), name = "Reproductive mode") +
  scale_color_colorblind() +
  ggsn::scalebar(
    data = locs,
    dist = 250,
    dist_unit = "km",
    height = 0.03,
    st.bottom = FALSE,
    st.dist = 0.05,
    st.size = 4,
    transform = TRUE,
    location = "bottomright" ,
    model = "WGS84",
    facet.var = c("genus_sp"),
    facet.lev = c("Tectarchus ovobessus")
  ) +
  theme_map() +
  theme(legend.position = "bottom", 
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        strip.text = element_text(size = 12, face = "bold.italic"))

pdf(file = here("output", "manuscript_figs", "map.pdf"),
    width = 10,
    height= 7)
ggsn::north2(ggp = full_map, x = 0.77, y = 0.2, symbol = 4)
dev.off()

## quartz_off_screen 
##                 2

ggsn::north2(ggp = full_map, x = 0.77, y = 0.2, symbol = 4)

Marginality scores

Caption: Marginality estimates for each species-reproductive mode combination. In all cases, sexual individuals occupy more marginal environmental space than asexual individuals. Ecological niche factor analysis was performed to determine if a reproductive mode occupies more marginal environmental space than the background environment of the entire species. Statistical significance is indicated by asterisks, where * indicates p < 0.05 and ** indicates p < 0.01.

Note: going to include asterisks to indicate statistical significance or mention in the caption and reference the table if many values are significant. Will insert significance asterisks outside of R.

# get all files to read in
dir_list <-
  dir_ls(spread_path, glob = "*_marginality_score.csv",
         recurse = TRUE) 

dir_list <- dir_list[!str_detect(dir_list, "subset")]

# get reproductive mode. A bit hacky, but I don't want to spend too much time with this. _sexual = sexual, while sexual = asexual
m <-
  map(dir_list, function(x)
    str_extract(x, "_*sexual")) %>% unlist()

# get genus and species numbers for all of the values. Have to repeat each twice since we're reading in two spreadsheets per species
genus_sp_names <-
  dir_ls(spread_path) %>% 
  str_extract(pattern = "[^/]+$") %>% 
  rep(each = 2)

  # remove elevation csv from list
genus_sp_names <- genus_sp_names[genus_sp_names != "elevation.csv"]


# read in all csv files and bind the rows. make a new column with reproductive mode and genus_sp name.
all_marg <- map(dir_list, read_csv) %>%
  bind_rows() %>%
  mutate(
    repro_mode = case_when(m == "sexual" ~ "asexual",
                           m == "_sexual" ~ "sexual") %>%
      str_to_sentence(),
    genus = str_split_fixed(genus_sp_names, pattern = "_", n = 2)[,1] %>% 
      str_extract("^[a-z]") %>% 
      str_to_sentence() %>% 
      str_c("."),
    species = str_split_fixed(genus_sp_names, pattern = "_", n = 2)[,2],
    genus_sp = str_c(genus, species, sep = " ")
  ) %>%
  tidyr::pivot_wider(names_from = repro_mode, values_from = marginality) %>%
  rowwise() %>%
  mutate(mean_marg = mean(c(Asexual, Sexual))) %>%
  ungroup() %>%
  mutate(genus_sp = fct_reorder(genus_sp, mean_marg))

marg_lollipop <-
  ggplot(data = all_marg, aes(x = repro_mode, y = marginality, group = genus_sp)) +
  geom_segment(aes(
    x = genus_sp,
    xend = genus_sp,
    y = Asexual,
    yend = Sexual
  )) +
  geom_point(
    aes(x = genus_sp, y = Asexual),
    shape = 21,
    size = 4,
    fill = "white"
  ) +
  geom_point(aes(x = genus_sp, y = Sexual), shape = 19, size = 4) +
  labs(x = "Species", y = "Marginality") +
  coord_flip() +
  theme_bw() + 
  theme(axis.title = element_text(size = 14),
        axis.text = element_text(size = 10),
        axis.text.y = element_text(face = "italic"))

pdf(file = here("output", "manuscript_figs", "marginality.pdf"),
    width = 7,
    height = 4.32623792125) # divide 7 by the golden ratio
marg_lollipop
dev.off()

## quartz_off_screen 
##                 2

marg_lollipop

Stability

Note: below the statistics are stability maps to put these results into context

Caption: Distribution of environmental stability since the last glacial maximum (LGM) for each species-reproductive mode combination. Stability was estimated for average annual temperature, average annual precipitation, and overall stability, which is temperature multiplied by precipitation. Time periods for the Statistical significance is indicated by asterisks, where * indicates p < 0.05 and ** indicates p < 0.01.

Note: I’m putting asterisks in outside of R. The top row are all significantly different across metrics (temp, precip, overall). T. ovobessus is significantly different for temperature.

# get files to read in for each metric (temperature, precipitation, overall)
temp_stab_list <-
  dir_ls(spread_path, glob = "*temp_stability_df.csv",
         recurse = TRUE)

precip_stab_list <-
  dir_ls(spread_path, glob = "*precip_stability_df.csv",
         recurse = TRUE)

overall_stab_list <-
  dir_ls(spread_path, glob = "*overall_stability_df.csv",
         recurse = TRUE)


read_csv_sp_name <- function(file) {
  s <- read_csv(file)
  
  genus_sp_file <- str_split(file, pattern = "/") %>% unlist()
  
  genus_sp <- genus_sp_file[length(genus_sp_file) - 1]
  
  mutate(s, genus_sp = genus_sp)
}

temp_stab <- map(temp_stab_list, read_csv_sp_name) %>% bind_rows()
precip_stab <- map(precip_stab_list, read_csv_sp_name) %>% bind_rows()
overall_stab <- map(overall_stab_list, read_csv_sp_name) %>% bind_rows()


# this assumes all individuals are in the same order. A glance at each data set confirms this, but make sure not to rearrange observation
all_stab <- bind_cols(temp_stab, 
                      precip_stab %>% dplyr::select(contains("stability")),
                      overall_stab %>% dplyr::select(contains("stability"))) %>% 
  pivot_longer(cols = contains("stability"), names_to = "metric", values_to = "stability") %>% 
  mutate(metric = str_remove(metric, "_stability_scaled"),
         genus_sp = str_replace(genus_sp, "_", " ") %>% str_to_sentence(),
         reproductive_mode = str_to_sentence(reproductive_mode)) %>% 
  mutate(metric = case_when(
    metric == "overall" ~ "Overall",
    metric == "temp" ~ "Temperature",
    metric == "precip" ~ "Precipitation"
  ))

stab_box <- ggplot(data = all_stab, aes(x = metric, y = stability, fill = reproductive_mode)) + 
  geom_boxplot() +
  scale_fill_manual(values = c("transparent", "darkgrey")) +
  labs(y = "Scaled Stability", fill = "Reproductive mode") +
  facet_wrap(~ genus_sp) + 
  theme_few() +
  theme(legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        strip.background = element_rect(fill = "lightgrey", color = "black"),
        strip.text = element_text(size = 12, face = "bold.italic"),
        panel.grid.major.y = element_line( size=.1, color="gray" ),
        legend.position = "top",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 16))

pdf(file = here("output", "manuscript_figs", "stability.pdf"),
    width = 9,
    height= 7)
stab_box
dev.off()

## quartz_off_screen 
##                 2

stab_box

Statistics

Calculating some statistics from the data (sample size, assumptions, pairwise comparisons, etc.)

# counts per species per reproductive mode
paste0("Sample size per species per reproductive mode.")

## [1] "Sample size per species per reproductive mode."

all_stab %>% 
  group_by(genus_sp, reproductive_mode) %>% 
  tally()

## # A tibble: 12 x 3
## # Groups:   genus_sp [6]
##    genus_sp               reproductive_mode     n
##    <chr>                  <chr>             <int>
##  1 Argosarchus horridus   Asexual              96
##  2 Argosarchus horridus   Sexual              213
##  3 Asteliaphasma jucundum Asexual              21
##  4 Asteliaphasma jucundum Sexual              159
##  5 Clitarchus hookeri     Asexual             147
##  6 Clitarchus hookeri     Sexual              351
##  7 Niveaphasma annulata   Asexual             111
##  8 Niveaphasma annulata   Sexual              114
##  9 Tectarchus huttoni     Asexual              33
## 10 Tectarchus huttoni     Sexual               45
## 11 Tectarchus ovobessus   Asexual              60
## 12 Tectarchus ovobessus   Sexual               27

# df is all_stab or all_seasonality , genus_species is formatted "Genus species", metric is "Temperature", "Precipitation", or "Overall", scale is "stability" or "seasonality
check_assumptions <- function(df, genus_species, metric, scale = "stability") {
  stab <- df %>% 
    filter(genus_sp == genus_species, metric == metric)
  
  asexual <- stab %>%
    filter(reproductive_mode == "Asexual") %>%
    dplyr::select(scale) %>%
    as_vector()

  sexual <- stab %>%
    filter(reproductive_mode == "Sexual") %>%
    dplyr::select(scale) %>%
    as_vector()

  # # shapiro wilk test for normality for both reproductive modes
  shapiro_asexual <-
    shapiro.test(asexual) %>%
    broom::tidy() %>%
    dplyr::select(p.value, method) %>%
    mutate(genus_sp = genus_species,
           reproductive_mode = "Asexual",
           metric = metric)

  shapiro_sexual <-
    shapiro.test(sexual) %>%
    broom::tidy() %>%
    dplyr::select(p.value, method) %>%
    mutate(genus_sp = genus_species,
           reproductive_mode = "Sexual",
           metric = metric)

  var_df <- var.test(asexual, sexual) %>%
    broom::tidy() %>%
    dplyr::select(p.value,
           method) %>% 
    mutate(genus_sp = genus_species,
           reproductive_mode = "Both",
           metric = metric)

  all_tests <- bind_rows(shapiro_asexual,
                         shapiro_sexual,
                         var_df)

  all_tests

}


temp_assumptions_stab <-
  map_df(
    unique(all_stab$genus_sp),
    ~ check_assumptions(
      df = all_stab,
      genus_species = .x,
      metric = "Temperature",
      scale = "stability"
    )
  )


precip_assumptions_stab <-
  map_df(
    unique(all_stab$genus_sp),
    ~ check_assumptions(
      df = all_stab,
      genus_species = .x,
      metric = "Precipitation",
      scale = "stability"
    )
  )

overall_assumptions_stab <- 
  map_df(
    unique(all_stab$genus_sp),
    ~ check_assumptions(
      df = all_stab,
      genus_species = .x,
      metric = "Overall",
      scale = "stability"
    )
  )


all_assumptions_stab <-
  bind_rows(temp_assumptions_stab, 
            precip_assumptions_stab, 
            overall_assumptions_stab)

# see which assumptions are violated
paste0("Assumptions test. Too many are violated- going to go with a Welch's t-test. Sampling is uneven anyways.")

## [1] "Assumptions test. Too many are violated- going to go with a Welch's t-test. Sampling is uneven anyways."

all_assumptions_stab %>% 
  filter(p.value < 0.05)

## # A tibble: 36 x 5
##     p.value method                  genus_sp          reproductive_mo… metric   
##       <dbl> <chr>                   <chr>             <chr>            <chr>    
##  1 5.76e- 5 Shapiro-Wilk normality… Argosarchus horr… Asexual          Temperat…
##  2 9.86e- 8 Shapiro-Wilk normality… Argosarchus horr… Sexual           Temperat…
##  3 1.20e- 2 Shapiro-Wilk normality… Asteliaphasma ju… Asexual          Temperat…
##  4 3.77e-10 Shapiro-Wilk normality… Asteliaphasma ju… Sexual           Temperat…
##  5 6.55e- 4 F test to compare two … Asteliaphasma ju… Both             Temperat…
##  6 3.58e- 9 Shapiro-Wilk normality… Clitarchus hooke… Asexual          Temperat…
##  7 2.69e-14 Shapiro-Wilk normality… Clitarchus hooke… Sexual           Temperat…
##  8 7.43e- 5 F test to compare two … Clitarchus hooke… Both             Temperat…
##  9 9.95e- 8 Shapiro-Wilk normality… Niveaphasma annu… Asexual          Temperat…
## 10 2.43e-10 Shapiro-Wilk normality… Niveaphasma annu… Sexual           Temperat…
## # … with 26 more rows

# Since most species have unequal sample sizes and have at least one assumption violated, I'm going to use Welch's t-tests

run_ttests_stability <- function(df, genus_species, met) {
  scale_df <- df %>%
    filter(genus_sp == genus_species,
           metric == met)
  test <- t.test(stability ~ reproductive_mode,
         data = scale_df,
         var.equal = FALSE) %>% 
    broom::tidy() %>% 
    mutate(genus_sp = genus_species,
           metric = met)
  test
}

temp_ttests_stab <- map_df(
    unique(all_stab$genus_sp),
    ~ run_ttests_stability(
      df = all_stab,
      genus_species = .x,
      met = "Temperature"
    )
  )

precip_ttests_stab <- map_df(
    unique(all_stab$genus_sp),
    ~ run_ttests_stability(
      df = all_stab,
      genus_species = .x,
      met = "Precipitation"
    )
  )

overall_ttests_stab <- map_df(
    unique(all_stab$genus_sp),
    ~ run_ttests_stability(
      df = all_stab,
      genus_species = .x,
      met = "Overall"
    )
  )

all_ttests_stab <- bind_rows(temp_ttests_stab, 
                        precip_ttests_stab,
                        overall_ttests_stab)

Table

All statistically significant (P < 0.05) comparisons.

all_ttests_stab %>% 
  filter(p.value < 0.05) %>% 
  dplyr::select(Species = genus_sp, Metric = metric, Estimate = estimate, Statistic = statistic, Pvalue = p.value) %>% 
  group_by(Species) %>% 
  gt() %>% 
  tab_header(title = "Welch's two sample t-tests comparing occupied stability between reproductive modes") %>% 
  tab_footnote(footnote = "Direction of the estimate and statistic indicate the difference in stability from sexual to asexual reproductive modes.", locations = cells_column_labels(
      columns = vars(Estimate, Statistic))) %>% 
  fmt_number(columns = vars(Estimate, Statistic, Pvalue),
             decimals = 3) %>% 
  tab_style(style = list(cell_text(style = "italic"),
                         cell_text(weight = "bold")),
            locations = cells_row_groups())

Welch's two sample t-tests comparing occupied stability between reproductive modes

Metric	Estimate¹	Statistic¹	Pvalue
Argosarchus horridus
Temperature	−0.075	−2.862	0.006
Precipitation	−0.050	−2.277	0.025
Overall	−0.036	−3.209	0.002
Asteliaphasma jucundum
Temperature	0.214	3.475	0.012
Precipitation	0.091	2.561	0.041
Overall	0.063	4.324	0.004
Clitarchus hookeri
Temperature	0.173	7.904	0.000
Precipitation	0.052	3.113	0.002
Overall	0.046	5.660	0.000
Tectarchus ovobessus
Temperature	0.132	2.803	0.013
¹ Direction of the estimate and statistic indicate the difference in stability from sexual to asexual reproductive modes.

Maps

Overall stability:

Temperature stability:

Precipitation stability:

Seasonality

I z-normalized the temperature and precipitation seasonality values for comparison and calculated “overall seasonality” as precip_seasonality * temp_seasonality.

# get files to read in for each metric (temperature, precipitation, overall)
temp_seasonality_list <-
  dir_ls(spread_path, glob = "*temp_seasonality_df.csv",
         recurse = TRUE)

precip_seasonality_list <-
  dir_ls(spread_path, glob = "*precip_seasonality_df.csv",
         recurse = TRUE)

overall_seasonality_list <-
  dir_ls(spread_path, glob = "*overall_seasonality_df.csv",
         recurse = TRUE)


read_csv_sp_name <- function(file) {
  s <- read_csv(file)
  
  genus_sp_file <- str_split(file, pattern = "/") %>% unlist()
  
  genus_sp <- genus_sp_file[length(genus_sp_file) - 1]
  
  mutate(s, genus_sp = genus_sp)
}

# function to z-normalize data
z_norm <- function(x) {
  (x - min(x)) / (max(x) - min(x))
}

# the NA I'm removing is the same observation for both data sets.
# normalize seasonality data for comparison and add in overall seasonality column.
temp_seasonality <- map(temp_seasonality_list, read_csv_sp_name) %>% 
  bind_rows() %>% 
  na.omit() %>% 
  mutate(temp_seasonality_norm = z_norm(temp_seasonality))


precip_seasonality <- map(precip_seasonality_list, read_csv_sp_name) %>%
  bind_rows() %>%
  na.omit() %>% 
  mutate(precip_seasonality_norm = z_norm(precip_seasonality))

overall_seasonality <- tibble(
  overall_seasonality_norm = precip_seasonality$precip_seasonality_norm * temp_seasonality$temp_seasonality_norm)

# this assumes all individuals are in the same order. A glance at each data set confirms this, but make sure not to rearrange observation
all_seasonality <- bind_cols(temp_seasonality %>% 
                               dplyr::select(-temp_seasonality), 
                             precip_seasonality %>%
                               dplyr::select(precip_seasonality_norm),
                             overall_seasonality %>% 
                               dplyr::select(overall_seasonality_norm)) %>% 
  pivot_longer(cols = contains("seasonality"),
               names_to = "metric",
               values_to = "seasonality") %>%
  mutate(metric = str_remove(metric, "_seasonality_norm"),
         genus_sp = str_replace(genus_sp, "_", " ") %>%
           str_to_sentence(),
         reproductive_mode = str_to_sentence(reproductive_mode)) %>%
  mutate(metric = case_when(
    metric == "temp" ~ "Temperature",
    metric == "precip" ~ "Precipitation",
    metric == "overall" ~ "Overall"
  ))

seasonality_box <- ggplot(data = all_seasonality, aes(x = metric, y = seasonality, fill = reproductive_mode)) + 
  geom_boxplot() +
  scale_fill_manual(values = c("transparent", "darkgrey")) +
  labs(y = "Seasonality", fill = "Reproductive mode") +
  facet_wrap(~ genus_sp) + 
  theme_few() +
  theme(legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        strip.background = element_rect(fill = "lightgrey", color = "black"),
        strip.text = element_text(size = 12, face = "bold.italic"),
        panel.grid.major.y = element_line( size=.1, color="gray" ),
        legend.position = "top",
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 16))

pdf(file = here("output", "manuscript_figs", "seasonality.pdf"),
    width = 9,
    height= 7)
seasonality_box
dev.off()

## quartz_off_screen 
##                 2

seasonality_box

Statistics

Calculating some statistics from the data (sample size, assumptions, pairwise comparisons, etc.)

# counts per species per reproductive mode
paste0("Sample size per species per reproductive mode.")

## [1] "Sample size per species per reproductive mode."

all_seasonality %>% 
  group_by(genus_sp, reproductive_mode) %>% 
  tally()

## # A tibble: 12 x 3
## # Groups:   genus_sp [6]
##    genus_sp               reproductive_mode     n
##    <chr>                  <chr>             <int>
##  1 Argosarchus horridus   Asexual              96
##  2 Argosarchus horridus   Sexual              213
##  3 Asteliaphasma jucundum Asexual              21
##  4 Asteliaphasma jucundum Sexual              150
##  5 Clitarchus hookeri     Asexual             147
##  6 Clitarchus hookeri     Sexual              315
##  7 Niveaphasma annulata   Asexual             111
##  8 Niveaphasma annulata   Sexual              114
##  9 Tectarchus huttoni     Asexual              33
## 10 Tectarchus huttoni     Sexual               45
## 11 Tectarchus ovobessus   Asexual              60
## 12 Tectarchus ovobessus   Sexual               27

temp_assumptions_seasonality <-
  map_df(
    unique(all_seasonality$genus_sp),
    ~ check_assumptions(
      df = all_seasonality,
      genus_species = .x,
      metric = "Temperature",
      scale = "seasonality"
    )
  )


precip_assumptions_seasonality <-
  map_df(
    unique(all_seasonality$genus_sp),
    ~ check_assumptions(
      df = all_seasonality,
      genus_species = .x,
      metric = "Precipitation",
      scale = "seasonality"
    )
  )

overall_assumptions_seasonality <- 
  map_df(
    unique(all_seasonality$genus_sp),
    ~ check_assumptions(
      df = all_seasonality,
      genus_species = .x,
      metric = "Overall",
      scale = "seasonality"
    )
  )


all_assumptions_seasonality <-
  bind_rows(temp_assumptions_seasonality, 
            precip_assumptions_seasonality, 
            overall_assumptions_seasonality)

# see which assumptions are violated
paste0("Assumptions test. Too many are violated- going to go with a Welch's t-test. Sampling is uneven anyways.")

## [1] "Assumptions test. Too many are violated- going to go with a Welch's t-test. Sampling is uneven anyways."

all_assumptions_seasonality %>% 
  filter(p.value < 0.05)

## # A tibble: 42 x 5
##        p.value method                genus_sp          reproductive_mo… metric  
##          <dbl> <chr>                 <chr>             <chr>            <chr>   
##  1     5.84e-6 Shapiro-Wilk normali… Argosarchus horr… Asexual          Tempera…
##  2     2.31e-7 Shapiro-Wilk normali… Argosarchus horr… Sexual           Tempera…
##  3     1.24e-2 F test to compare tw… Argosarchus horr… Both             Tempera…
##  4     5.70e-3 Shapiro-Wilk normali… Asteliaphasma ju… Asexual          Tempera…
##  5     3.80e-4 Shapiro-Wilk normali… Asteliaphasma ju… Sexual           Tempera…
##  6     2.36e-2 F test to compare tw… Asteliaphasma ju… Both             Tempera…
##  7     1.85e-3 Shapiro-Wilk normali… Clitarchus hooke… Asexual          Tempera…
##  8     7.20e-8 Shapiro-Wilk normali… Clitarchus hooke… Sexual           Tempera…
##  9     3.90e-5 F test to compare tw… Clitarchus hooke… Both             Tempera…
## 10     1.66e-9 Shapiro-Wilk normali… Niveaphasma annu… Asexual          Tempera…
## # … with 32 more rows

# Since most species have unequal sample sizes and have at least one assumption violated, I'm going to use Welch's t-tests

run_ttests_seasonality <- function(df, genus_species, met) {
  seasonality <- df %>%
    filter(genus_sp == genus_species,
           metric == met)
  test <- t.test(seasonality ~ reproductive_mode,
         data = seasonality,
         var.equal = FALSE) %>% 
    broom::tidy() %>% 
    mutate(genus_sp = genus_species,
           metric = met)
  test
}

temp_ttests_seasonality <- map_df(
    unique(all_seasonality$genus_sp),
    ~ run_ttests_seasonality(
      df = all_seasonality,
      genus_species = .x,
      met = "Temperature"
    )
  )

precip_ttests_seasonality <- map_df(
    unique(all_seasonality$genus_sp),
    ~ run_ttests_seasonality(
      df = all_seasonality,
      genus_species = .x,
      met = "Precipitation"
    )
  )

overall_ttests_seasonality <- map_df(
    unique(all_seasonality$genus_sp),
    ~ run_ttests_seasonality(
      df = all_seasonality,
      genus_species = .x,
      met = "Overall"
    )
  )

all_ttests_seasonality <- bind_rows(temp_ttests_seasonality, 
                        precip_ttests_seasonality,
                        overall_ttests_seasonality)

Table

All statistically significant (P < 0.05) comparisons.

all_ttests_seasonality %>% 
  filter(p.value < 0.05) %>% 
  dplyr::select(Species = genus_sp, Metric = metric, Estimate = estimate, Statistic = statistic, Pvalue = p.value) %>% 
  group_by(Species) %>% 
  gt() %>% 
  tab_header(title = "Welch's two sample t-tests comparing occupied seasonality between reproductive modes") %>% 
  tab_footnote(footnote = "Direction of the estimate and statistic indicate the difference in seasonality from sexual to asexual reproductive modes.", locations = cells_column_labels(
      columns = vars(Estimate, Statistic))) %>% 
  fmt_number(columns = vars(Estimate, Statistic, Pvalue),
             decimals = 3) %>% 
  tab_style(style = list(cell_text(style = "italic"),
                         cell_text(weight = "bold")),
            locations = cells_row_groups())

Welch's two sample t-tests comparing occupied seasonality between reproductive modes

Metric	Estimate¹	Statistic¹	Pvalue
Clitarchus hookeri
Temperature	0.050	2.348	0.022
Niveaphasma annulata
Temperature	0.328	5.166	0.000
Overall	0.071	2.652	0.010
Argosarchus horridus
Precipitation	−0.095	−2.326	0.024
Tectarchus huttoni
Precipitation	−0.144	−2.324	0.030
Overall	−0.061	−2.608	0.016
¹ Direction of the estimate and statistic indicate the difference in seasonality from sexual to asexual reproductive modes.

Elevation

This plot illustrates elevational range differences among species and between reproductive modes within species to place niche and stability results in context.

# extract elevations (in meters) per locality from SRTM DEM at 80 m resolution. 864.8 m resolution at -45 degrees (closest to NZ). This takes a while, so I'm only doing it if the file isn't there already
if (!file.exists(here("output", "spreadsheets", "elevation.csv"))) {
  nz_dem <- raster(here("data", "NZ_80m_dem", "nz-80m-digital-elevation-model.tif"))
  locs_sp <- as_Spatial(locs)
  elev <- raster::extract(nz_dem, locs_sp) %>% 
    as_tibble() %>% 
    rename(elevation = value) %>% 
    bind_cols(locs)
  write_csv(elev, here("output", "spreadsheets", "elevation.csv"))
} else elev <- read_csv(here("output", "spreadsheets", "elevation.csv"))


ggplot() +
  geom_density_ridges(data = elev, 
                      aes(x = elevation, 
                          y = genus_sp, 
                          fill = reproductive_mode),
                      color = "black",
                      alpha = 0.7) +
  scale_fill_manual(values = c("white", "black")) +
  labs(x = "Elevation (m)", fill = "Reproductive mode") +
  theme_minimal() +
  theme(axis.title.y = element_blank())

Environmental space

Should I make a plot of environmental PCA space to illustrate marginality further? Is it superfluous?

Idea 1: Just use the hex plot with species env. space scatterplot on top that I’ve already made. Might choose a different way to represent density of background environmental space.

Idea 2: contour plot in PC space (the dimension reduction ENFA outputs or PCA): unique contours for each reproductive mode and background e-space.

Publication Figures

Connor French

Methods, Results, Figures

Methods

Sampling

Environmental data

Marginality

Stability

Seasonality

Results

Marginality

Stability

Seasonality

Figures and Tables

Sampling map

Marginality scores

Stability

Statistics

Table

Maps

Seasonality

Statistics

Table

Elevation

Environmental space