3 DIfference in diet between H. curtus and H. schsitosus

3.1 No. of samples collected by sea snake species

source("./Functions/setup.R")

# importing gut content data

gc = read.csv("./Data/Sea_snakes_gut_content_2018-19.csv") %>%
    filter(Source == "Gut content")  # removing specimens collected from fisheries landings

snakes%>%
  filter(Species == "Hydrophis curtus" | Species == "Hydrophis schistosus")%>%
  group_by(Species)%>%
  summarise(N = n(), # total number of individuals sampled
            n = sum((Gut.Content) == "Yes"), # number of individuals with gut content
            n_male = sum((Gut.Content) == "Yes" & Sex == "Male"), # males with gut content
            n_female = sum((Gut.Content) == "Yes" & Sex == "Female"), # females with content
            prop.gc = n/N) # propotion of individuals with gut content

Species	N	n	n_male	n_female	prop.gc
Hydrophis curtus	179	36	18	18	0.2011173
Hydrophis schistosus	605	93	55	36	0.1537190

A very low proportion of individuals had gut content present at the time of sampling. Sampling across species and sexes may be adequate for comparison.

3.2 Total number of prey specimens collected

gc %>%
    group_by(Snake.Species) %>%
    summarise(N = length(unique(GC.Code)))

Snake.Species	N
Hydrophis curtus	35
Hydrophis schistosus	89

3.3 Proportion of unidentified specimens

# calculating richness

gc%>%
  summarise(N_sp = length(unique(Prey.Species)), # number of prey species
            N_fam = length(unique(Prey.Family)), # number of prey families
            # proportion of unidentitified specimens
            unid_sp = sum(Prey.Species == "Unidentified" | Prey.Species == "", na.rm = T)*100/n(), ## species
            unid_fam = sum(Prey.Family == "Unidentified" | Prey.Family == "", na.rm = T)*100/n())%>% ## family
  # creating clean table
  gather()%>%
  mutate(Metric = ifelse(substr(key, 1, 1) == "N", "Richness", "% Unidentified"),
         Unit = ifelse(grepl(key, pattern = "sp", fixed = T), "Prey Species", "Prey Family"))%>%
  dplyr::select(Unit, Metric, value)%>%
  spread(Metric, value)

Unit	% Unidentified	Richness
Prey Family	29.45736	17
Prey Species	40.31008	26

A large portion of gut content specimens were unidentifiable. This is an unavoidable consquence of VGCA.

3.4 No. of snakes with mode than one prey specimen

gc %>%
    filter(Field.Code != "") %>%
    count(Field.Code, name = "Prey") %>%
    count(Prey > 1)

Prey > 1	n
FALSE	96
TRUE	8

3.5 Prey preference

We modified and used the the Index of relative importance according to Pinkas et al. 1971 to determine the difference in prey preference between H. curtus and H. schistosus

# calulating index of relative importance (Pinkas et al. 1971)

IRI <- gc%>%
  filter(Prey.Family != "")%>%
  filter(Snake.Species == "Hydrophis schistosus" | Snake.Species == "Hydrophis curtus")%>%
  group_by(Snake.Species)%>%
  mutate(Fr = length(unique(Field.Code)), # total number of snakes sampled
         W = sum(Weight..g., na.rm = T), # total weight of prey sampled
         N = n())%>% # total number of prey sampled
  group_by(Snake.Species, Prey.Family)%>%
  summarise(f = length(unique(Field.Code)), # number snakes in which prey family occured
            Fr = last(Fr),
            w = sum(Weight..g., na.rm = T), # weight of prey family
            W = last(W),
            n = n(), # number of individuals of prey family
            N = last(N))%>%
  group_by(Snake.Species, Prey.Family)%>%
  # caluclating percentages
  summarise(per.F = f*100/Fr,
            per.W = w*100/W,
            per.N = n*100/N,
            # caluclating IRI
            IRI = (per.N + per.W)*per.F)%>%
  ungroup()%>%
  complete(Prey.Family, nesting(Snake.Species), fill = list(0))%>%
  replace_na(replace = list(IRI = 0, per.F = 0, per.W = 0, per.N = 0))%>%
  arrange(Snake.Species, desc(IRI))%>%
  mutate(rank = 1:n())

# clean table

IRI%>%
  filter(IRI > 0)

Prey.Family	Snake.Species	per.F	per.W	per.N	IRI	rank
Unidentified	Hydrophis curtus	44.827586	15.2823920	40.625000	2506.193436	1
Clupeidae	Hydrophis curtus	6.896552	29.9003322	9.375000	270.864360	2
Engraulidae	Hydrophis curtus	10.344828	11.9601329	9.375000	220.708271	3
Carangidae	Hydrophis curtus	13.793103	3.3222591	12.500000	218.238057	4
Scombridae	Hydrophis curtus	6.896552	6.6445183	6.250000	88.927712	5
Cynoglossidae	Hydrophis curtus	6.896552	3.9867110	6.250000	70.598007	6
Serranidae	Hydrophis curtus	3.448276	16.6112957	3.125000	68.056192	7
Leiognathidae	Hydrophis curtus	6.896552	2.9900332	6.250000	63.724367	8
Tetraodontidae	Hydrophis curtus	3.448276	7.6411960	3.125000	37.124814	9
Nemipteridae	Hydrophis curtus	3.448276	1.6611296	3.125000	16.503895	10
Tetraodontidae	Hydrophis schistosus	34.246575	38.7307267	35.869565	2554.804519	17
Unidentified	Hydrophis schistosus	23.287671	12.9453244	21.739130	807.720182	18
Ariidae	Hydrophis schistosus	9.589041	8.8407094	8.695652	168.156891	19
Clupeidae	Hydrophis schistosus	6.849315	10.9456402	6.521739	119.639584	20
Plotosidae	Hydrophis schistosus	8.219178	7.7882440	6.521739	117.616299	21
Teraponidae	Hydrophis schistosus	6.849315	6.5673841	5.434783	82.206621	22
Serranidae	Hydrophis schistosus	5.479452	7.6303742	5.434783	71.589900	23
Leiognathidae	Hydrophis schistosus	4.109589	1.0524654	3.260870	17.726034	24
Synodontidae	Hydrophis schistosus	2.739726	1.1840236	2.173913	9.199826	25
Sillaginidae	Hydrophis schistosus	1.369863	2.4206704	1.086957	4.804968	26
Scombridae	Hydrophis schistosus	1.369863	1.2103352	1.086957	3.146975	27
Pleuronectidae	Hydrophis schistosus	1.369863	0.3683629	1.086957	1.993588	28
Carangidae	Hydrophis schistosus	1.369863	0.3157396	1.086957	1.921502	29

# plotting prey preference

IRI %>%
    ggplot(aes(reorder(Prey.Family, IRI), IRI, fill = Snake.Species)) + geom_col(col = "black", position = position_dodge(preserve = "single")) + scale_y_sqrt(name = "IRI (sq.rt.)") + 
    labs(x = "Prey Family") + scale_fill_brewer(palette = "Greys", name = "Snake Species") + theme(axis.text.x = element_text(hjust = 1, angle = 60), legend.text = element_text(face = "italic"))

# saving high res plot

ggsave(last_plot(), filename = "./Figures and Tables/figure2.tiff", width = 8, height = 6, dpi = 300)

Figure 2: Index of relative importance (IRI, square root transformed) (Pinkas et al., 1970) of prey families in sea snake diet.

H. schistosus showed preference for the Tetraodontid fish, while H. curtus did not exhibit any prey preference.

3.6 Range of prey found in H. curtus and H. shistosus gut

We compared the dietary breadth between snakes species using the richness and diversity of prey families.

# Creating data matrix for prey community analysis analysis

fam_spade = gc%>%
  filter(Prey.Family != "Unidentified", Prey.Family !="", Field.Code != "")%>% # removing unidentified specimens
  group_by(Snake.Species, Prey.Family)%>%
  summarise(n = n())%>% # total number of specimens per family in each snake species
  spread(Snake.Species, n, fill = 0)%>%
  column_to_rownames("Prey.Family")

3.6.1 Prey richness

gc%>%
  filter(Snake.Species == "Hydrophis schistosus" | Snake.Species == "Hydrophis curtus",
         Prey.Family != "Unidentified", Prey.Family != "")%>%
  group_by(Snake.Species)%>%
  summarise(Prey.Species = length(unique(Prey.Species)) - 1, # removing unidentified
            Prey.Families = length(unique(Prey.Family)))

Snake.Species	Prey.Species	Prey.Families
Hydrophis curtus	11	9
Hydrophis schistosus	19	12

3.6.2 Prey family diversity

# loading libraries

library(SpadeR)  ## for community indices

# calculating shanon diversity

HC_div = as.data.frame(Diversity(fam_spade[, 1])$Shannon_diversity) %>%
    rownames_to_column("Estimator") %>%
    mutate(Species = "Hydrophis curtus")

HS_div = as.data.frame(Diversity(fam_spade[, 2])$Shannon_diversity) %>%
    rownames_to_column("Estimator") %>%
    mutate(Species = "Hydrophis schistosus")

HC_div %>%
    full_join(HS_div) %>%
    dplyr::select(Species, everything()) %>%
    arrange(Estimator)

Species	Estimator	Estimate	s.e.	95%Lower	95%Upper
Hydrophis curtus	Chao & Shen	11.671	3.155	5.487	17.855
Hydrophis schistosus	Chao & Shen	7.427	1.158	5.158	9.696
Hydrophis curtus	Chao et al. (2013)	12.029	3.207	5.744	18.314
Hydrophis schistosus	Chao et al. (2013)	7.509	1.179	5.198	9.820
Hydrophis curtus	Jackknife	11.776	2.583	6.713	16.839
Hydrophis schistosus	Jackknife	7.301	1.126	5.093	9.508
Hydrophis curtus	MLE	8.073	1.181	5.758	10.388
Hydrophis schistosus	MLE	6.428	0.911	4.642	8.214

While H. schistosus fed on greater number of prey families, H. curtus showed greater evenness in prey preference.

Statistical comparison is difficult as each individual snake fed on only one type of specimen at time.

3.7 Diet overlap between H. curtus and H. shistosus

Number of overlapping prey families:

gc %>%
    filter(Snake.Species == "Hydrophis schistosus" | Snake.Species == "Hydrophis curtus", Prey.Family != "Unidentified", Prey.Family != "") %>%
    group_by(Prey.Family) %>%
    summarise(N_pred = length(unique(Snake.Species))) %>%
    ungroup() %>%
    summarise(Overlap = sum(N_pred > 1)) %>%
    as.numeric()

## [1] 6

Shoener’s (1971) resource overlap index (D):

gc%>%
  filter(Prey.Species != "Unidentified")%>%
  group_by(Snake.Species, Prey.Family)%>%
  summarise(n = n())%>% # number of speciemens per prey family in each snake species
  group_by(Snake.Species)%>%
  mutate(N = sum(n))%>% # total number of specimens collected from each snake species
  group_by(Snake.Species, Prey.Family)%>%
  summarise(p = n/N)%>% # relative proportion of each prey family
  spread(Snake.Species, p, fill = 0)%>%
  rename("HC" = 'Hydrophis curtus', "HS" = 'Hydrophis schistosus')%>%
  mutate(P = abs(HC - HS))%>% # difference in relactive proportion prey items 
  summarise(D = 1 - 0.5*sum(P))%>% # Schoener's resource overlap index
  as.numeric()

## [1] 0.2380952

Morista - Horn overlap:

prey_sim <- SimilarityPair(X = fam_spade, datatype = "abundance")

prey_sim$Empirical_relative

##                            Estimate      s.e.   95%.LCL   95%.UCL
## C12=U12(q=1,Horn)         0.4419859 0.1103973 0.2256071 0.6583647
## C22(q=2,Morisita)         0.2777801 0.1131096 0.0560853 0.4994749
## U22(q=2,Regional overlap) 0.4347854 0.1461632 0.1483056 0.7212653
## ChaoJaccard               0.4946000 0.1537000 0.1933480 0.7958520
## ChaoSorensen              0.6619000 0.1736000 0.3216440 1.0021560

There was limited overlap in prey between H. schistosus and H. curtus

3.7.1 Testing segregation in prey between H. curtus and H. schistosus

# loading libraries

library(vegan)

# formating data for vegan

fam_simboo <- gc %>%
    filter(Prey.Family != "Unidentified", Prey.Family != "", Field.Code != "") %>%
    group_by(Snake.Species, Field.Code, Prey.Family) %>%
    summarise(n = n()) %>%
    spread(Prey.Family, n, fill = 0) %>%
    ungroup()

# PERMANOVA to compare composition

set.seed(2)

permanova_sp <- adonis2(fam_simboo[, 3:length(fam_simboo)] ~ fam_simboo$Snake.Species, data = fam_simboo[, 1])

## clean table

broom::tidy(permanova_sp)

term	df	SumOfSqs	R2	statistic	p.value
fam_simboo$Snake.Species	1	1.556349	0.0511228	3.771403	0.002
Residual	70	28.886984	0.9488772
Total	71	30.443333	1.0000000

Prey composition was significantly different between the two snake species.

3.8 Difference in size selectivity of prey between H. curtus and H. schistosus

# summary

gc%>%
  filter(Snake.Species == "Hydrophis schistosus" | Snake.Species == "Hydrophis curtus",
         Condition < 3)%>% # removing very digested specimens 
  group_by(Snake.Species)%>%
  skimr::skim(Maximum.Body.Girth..cm.)%>%
  skimr::yank("numeric")%>%
  dplyr::select(Snake.Species:p100)

Variable type: numeric

Snake.Species	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100
Hydrophis curtus	2	0.86	2.77	0.76	2.0	2.08	2.65	3.2	4.5
Hydrophis schistosus	0	1.00	3.46	1.17	0.6	3.00	3.60	4.3	5.5

# testing diff max prey width by species

gc%>%
  filter(Snake.Species == "Hydrophis schistosus" | Snake.Species == "Hydrophis curtus",
         Condition < 3, # removing digested specimens 
         !is.na(Maximum.Body.Girth..cm.))%>% # removing unrecorded data
  dplyr::select(Snake.Species, Maximum.Body.Girth..cm.)%>%
  droplevels()%>%
  nest()%>%
  mutate(test = map(data, ~t.test(Maximum.Body.Girth..cm. ~ Snake.Species, data = .)), # t test
         sumry = map(test, broom::tidy),
         d = map(data, ~lsr::cohensD(Maximum.Body.Girth..cm. ~ Snake.Species, data = .)))%>% # effect size
  dplyr::select(sumry, d)%>%
  unnest()%>%
  dplyr::select(estimate1:parameter, d)%>%
  rename(`Hydrophis curtus` = estimate1,
         `Hydrophis schistosus` = estimate2)

Hydrophis curtus	Hydrophis schistosus	statistic	p.value	parameter	d
2.76775	3.463946	-2.386272	0.0237889	28.90876	0.6420566

H. schistosus fed on larger prey than H. curtus