Chapter 3 Results
Rhys Jenson
2025-10-09
| Sample ID | δ13C (‰) | δ15N (‰) |
|---|---|---|
| HW-20240911-GN-F1-FEC-SIA | -18.35 | 10.84 |
| HW-20240822-GN-F1-FEC-SIA | -17.78 | 10.54 |
| HW-20240822-GN-F4-FEC-SIA | -17.87 | 11.53 |
| HW-20240822-GN-F5-FEC-SIA | -17.93 | 11.06 |
| HW-20240823-GN-F1-FEC-SIA | -17.57 | 11.56 |
| HW-20240823-CC-F2-FEC-SIA | -16.43 | 10.64 |
| HW-20240912-GN-F1-FEC-SIA | -18.39 | 10.42 |
| HW-20240712-GN-F2-FEC-SIA | -19.11 | 11.51 |
| HW-20240821-GN-F1-FEC-SIA | -16.47 | 11.53 |
| HW-20240821-GN-F2-FEC-SIA | -16.86 | 11.40 |
| HW-20240821-GN-F3-FEC-SIA | -17.74 | 14.29 |
| HW-20240717-GN-F1-FEC-SIA | -19.22 | 12.18 |
| HW-20240815-GN-F1-FEC-SIA | -17.37 | 13.05 |
Part 1) Stable Isotope Modelling
- Check the effects of lipid extraction on δ13C and δ15N isotopes within fecal samples.
Show t-test details (lipid extraction and δ13C)
Paired t-test
data: wide$LE and wide$NLE
t = 8.4323, df = 12, p-value = 2.185e-06
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
1.070203 1.815951
sample estimates:
mean difference
1.443077
Show t-test details (Lipid extraction and δ15N)
Paired t-test
data: wide$LE and wide$NLE
t = 8.4323, df = 12, p-value = 2.185e-06
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
1.070203 1.815951
sample estimates:
mean difference
1.443077
Lipid extraction helps correct artificially depleted fecal ¹³C values. However, δ¹⁵N appears affected by lipid extraction as well, perhaps some mechanism is unintentionally extracting nitrogen-containing compounds. Therefore, we will use LE values for δ13C and NLE values for δ15N.
- Lets check if our humpback whale fecal values (δ13C and δ15N) vary significantly between months. The data set violated at least one ANOVA assumption of normality and equal variance with between N and C values, so Kruskal-Wallis test was used.
| Analyte | Kruskal–Wallis χ² | df | p-value |
|---|---|---|---|
| δ13C (‰) | 10.04 | 2 | 0.007 |
| δ15N (‰) | 4.05 | 2 | 0.132 |
Month appears to affect our lipid-extracted δ13C values. However, we have extremely uneven group sizes (N=10 for august, N=2 for July, N=2 for September) and low statistical power. Same is the case with the NLE δ15N values. Therefore, cannot rely on statistical inference in this case. There is not enough power to say isotopic signature of feces changed significantly between months.
- We dont have many isotope values for prey from our study zone, so we will apply a baseline shift correction to a larger 2023 Strait of Georgia (SOG) dataset. Pacific Krill act as an isotopic baseline in both the SOG and our study zone, Juan de Fuca (JdF). Using the shift in 13C and 15N between krill in both sites, we can apply an adjustion to prey species collected from the SOG in 2023, to represent what isotopic prey signatures would be seen in JdF 2024.
Our time and location adjust values for krill, juvenile herring, and adult herring are:
Show Krill-Adjusted Data
| d13C_mean_adj | d13C_sd_adj | d15N_mean_adj | d15N_sd_adj |
|---|---|---|---|
| -17.75 | 1.37 | 8.95 | 1.35 |
Show Herring-Adjusted Data
| d13C_mean_adj | d13C_sd_adj | d15N_mean_adj | d15N_sd_adj |
|---|---|---|---|
| -14.75 | 1.08 | 13.11 | 0.46 |
Show Juvenile Herring-Adjusted Data
| d13C_mean_adj | d13C_sd_adj | d15N_mean_adj | d15N_sd_adj |
|---|---|---|---|
| -16.4 | 1.15 | 12.18 | 0.39 |
- Now that we have three sources (with k independent tracers, you can uniquely solve for k + 1 sources), lets generate a Bayesian stable isotope mixing model (SIMM) using lipid-extracted δ13C values and non-lipid extracted δ15N values.
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 26
## Unobserved stochastic nodes: 31
## Total graph size: 152
##
## Initializing model
##
## Summary for 1## mean sd
## deviance 82.725 4.506
## Krill 0.754 0.081
## Adult Herring 0.095 0.055
## Juvenile Herring 0.151 0.085
## sd[d13C] 0.953 0.443
## sd[d15N] 1.064 0.485##
## Summary for 1## 2.5% 25% 50% 75% 97.5%
## deviance 76.497 79.458 81.963 85.097 93.588
## Krill 0.586 0.701 0.756 0.809 0.907
## Adult Herring 0.016 0.051 0.087 0.129 0.219
## Juvenile Herring 0.024 0.086 0.139 0.205 0.341
## sd[d13C] 0.281 0.638 0.892 1.208 2.004
## sd[d15N] 0.321 0.722 0.987 1.334 2.180
## Most popular orderings are as follows:## Probability
## Krill > Juvenile Herring > Adult Herring 0.6643
## Krill > Adult Herring > Juvenile Herring 0.3350
## Juvenile Herring > Krill > Adult Herring 0.0007
- Sensitivity analysis of our SIMM. Check how sensitive the model outputs are to the fecal correction factors. We are not going to test priors since we do not have diet data for this region. Rorqual diet is expected to vary drastically based on region around Vancouver Island so we believe using a generalist prior is best suited to let the actual data speak for themselves.
By increasing the 15N fecal-correction factor by 50%, the mean predicted proportion of krill consumed increases 11.4513981 %. By decreasing the 15N fecal-correction factor by 50%, the mean predicted proportion of krill consumed decreases 22.3701731 %. Therefore, the SIMM is fairly sensitive to the 15N fecal-correction factor. This value is predicted from previous studies. Accuracy of model would improve via study of wild humbpack. A deceased individual would be needed to compare the SI of full prey items from stomach and fecal samples from distal end of colon - this study was performed on a fin whale but they failed to included NLE 15N values.
This is the extent of SIA fecal analyses available to us. Trophic position cannot be determined since feces is not an assimilated tissue. Nitrogen isnt affect by lipid content, so this works well.
Part 2) Lets look at the results from our shotgun metagenomics.
Figure 7: Stacked bar plot showing the proportional composition (%) of prey DNA reads identified by shotgun metagenomics in each fecal sample. Each bar represents one whale fecal sample, with colors denoting different prey taxa or functional prey groups (e.g., krill, forage fish)
| SampleID | Top1 | Top2 | Top3 |
|---|---|---|---|
| mn24_194F1 | krill (92.26%) | Oncorhynchus keta (5.20%) | Clupea pallasii (2.23%) |
| mn24_194F2 | krill (69.18%) | Clupea pallasii (30.82%) | Ammodytes personatus (0.00%) |
| mn24_199F1 | krill (84.02%) | Clupea pallasii (15.98%) | Ammodytes personatus (0.00%) |
| mn24_228F1 | Clupea pallasii (55.97%) | krill (21.23%) | Engraulis mordax (4.20%) |
| mn24_234F1 | krill (66.59%) | Sardinops sagax (6.18%) | Clupea pallasii (5.08%) |
| mn24_234F2 | krill (100.00%) | Ammodytes personatus (0.00%) | Clupea pallasii (0.00%) |
| mn24_234F3 | Clupea pallasii (78.25%) | krill (19.74%) | Sardinops sagax (0.41%) |
| mn24_235F1 | krill (87.43%) | Clupea pallasii (1.87%) | Sardinops sagax (1.42%) |
| mn24_235F2 | krill (78.92%) | Ammodytes personatus (2.64%) | Sardinops sagax (2.49%) |
| mn24_235F3 | krill (98.88%) | Clupea pallasii (0.96%) | Sardinops sagax (0.16%) |
| mn24_235F4 | krill (52.61%) | Clupea pallasii (47.39%) | Ammodytes personatus (0.00%) |
| mn24_235F5 | krill (73.80%) | Clupea pallasii (26.20%) | Ammodytes personatus (0.00%) |
| mn24_236F1 | krill (71.19%) | Clupea pallasii (28.46%) | Ammodytes personatus (0.07%) |
| mn24_236F2 | krill (85.37%) | Engraulis mordax (2.63%) | Ammodytes personatus (2.15%) |
| mn24_255F1 | krill (77.36%) | Clupea pallasii (22.64%) | Ammodytes personatus (0.00%) |
| mn24_256F1 | krill (100.00%) | Ammodytes personatus (0.00%) | Clupea pallasii (0.00%) |
| mn24_258F1 | Clupea pallasii (88.67%) | krill (11.32%) | Ammodytes personatus (0.00%) |
| mn24_261F1 | Clupea pallasii (83.88%) | krill (11.08%) | Sardinops sagax (1.31%) |
| mn24_261F2 | Clupea pallasii (100.00%) | Ammodytes personatus (0.00%) | Engraulis mordax (0.00%) |
| mn24_263F1 | krill (56.03%) | Clupea pallasii (13.45%) | Ammodytes personatus (5.30%) |
Figure 8: Mean relative abundance (%) of prey items identifed with WSM across 20 humpback whale fecal samples
Across 20 fecal samples, krill made up 62.9% of the dietary composition revealed by WSM. Herring made up 30.2%.
If we want to assess variance of samples, we must acknowledge the compositional (not absolute) nature of shotgun metagenomic data
“In summary, the analysis of compositional data by traditional methods can appear to give satisfactory results. However, these results can be misleading and unpredictable. Compositionally appropriate tools exist as drop-in replacements at each stage of the analysis as shown in figure below.”
(Gloor et al., 2017, PLOS Computational Biology)
More info
In an ecological study it is possible for many different species to co-exist, and their absolute abundance may be important. For example, in an area containing only tigers, it is important to know if the population size is sufficient to maintain needed genetic diversity for long-term survival (Shaffer, 1981). However, the abundance of one species may not influence the abundance of another; the area may contain both tigers and ladybugs, and the migration of several ladybugs into the area would not be expected to affect the number of tigers. The assumption of true independence can not hold in high-throughput sequencing (HTS) experiments because the sequencing instruments can deliver reads only up to the capacity of the instrument. Thus, it is proper to think of these instruments as containing a fixed number of slots which must be filled. Returning to our tiger and ladybug analogy, the migration of ladybugs into an area containing a fixed number of slots that are already filled must displace either tigers or ladybugs from the occupied slots. This analogy extents, without restriction, to any number of taxa, and to any fixed capacity instrument (Aitchison, 1986; Lovell et al., 2011;Friedman and Alm, 2012; Fernandes et al., 2013, 2014; Lovell et al., 2015; Mandal et al., 2015; Gloor et al., 2016a,b; Gloor and Reid, 2016; Tsilimigras and Fodor, 2016). Thus, the total read count observed in a HTS run is a fixed-size, random sample of the relative abundance of the molecules in the underlying ecosystem. Moreover, the count can not be related to the absolute number of molecules in the input sample as shown in Figure 1. This is implicitly acknowledged when microbiome datasets are converted to relative abundance values, or normalized counts,or are rarefied (McMurdie and Holmes, 2014; Weiss et al., 2017) prior to analysis. Thus the number of reads obtained is irrelevant, and contains only information on the precision of the estimate (Fernandes et al., 2013). Data that are naturally described as proportions or probabilities, or with a constant or irrelevant sum, are referred to as compositional data.
Figure 9: PCoA plot of Aitchison distance for beta diversity comparisons of prey OTUs at the species level detected using whole metagenomic shotgun sequencing (WMS). Statistical analysis for beta diversity was performed using PERMANOVA to determine significance differences (P-value) and percentage of the variance explained (R²) between the groups. The right plot is grouped by month and the right plot is grouped by relative krill proportion
Extra plot
Figure 9b: PCoA plot of Aitchison distance for beta diversity comparisons of prey OTUs at the species level detected using whole metagenomic shotgun sequencing (WMS) grouped by sample volume. Statistical analysis for beta diversity was performed using PERMANOVA to determine significance differences (P-value) and percentage of the variance explained (R²) between the groups
Principal component analysis of CLR-transformed prey compositions showed that the first two components explained 74.7% of total variance, indicating that most compositional variability among samples was captured by a single dominant gradient. However, samples did not form discrete clusters by month, suggesting that prey composition varied continuously among whales rather than being structured by sampling period. Upon grouping by krill proportion, we see that Axis 1 can be explained well be the proportion of krill in diet. Therefore, variance in diet composition is mainly attributed to relative amount of krill in diet. If whales were more generalist feeders, the proportion of a single prey item would not cause much composition change in diet, however, it is clear here that whales are mainly eating one or two prey species. The amount of krill consumed is not explained by sampling month. The variance in krill consumed can is not explained by any measured variables in this study, but is likley controlled by local oceanographic factors such as upwellings and additional interspecific differences in whale feeding preferences. Lastly, points were grouped into sample volume to assess the affect of diet composition on collected fecal volume. There was no apparent relationship.
Average distance to the global centroid in Aitchison space was moderate (2.6630674 ± 0.7115093), indicating that some fecal samples were compositionally distinct but overall diet composition among individuals was somewhat homogeneous.
Figure 10. Aitchison distance to the global centroid for each sample. Bars show per-sample distinctness; dashed line indicates the mean across samples; shaded band is mean ± 1 SD. Higher values indicate samples whose prey composition (relative read ratios) deviates more from the dataset’s average diet (global centroid) in Aitchison space.
The mean global distance to centroid (weighted by sample size (N), not group size, is 2.6706746
Lets test whether group means differ (are monthly diets different on average?)Show PERMANOVA (Month) table
| Term | Df | Sum_of_Squares | R2 | F_value | Pr(>F) |
|---|---|---|---|---|---|
| Model | 2 | 27.734 | 0.184 | 1.922 | 0.084 |
| Residual | 17 | 122.680 | 0.816 | NA | NA |
| Total | 19 | 150.414 | 1.000 | NA | NA |
An assumption of PERMANOVA is within-group variance (dispersion) is roughly equal across groups. Below, this is confirmed. Our group sizes are unequal however, so this could affect robsutness. Might be worth using equal group weight PERMANOVA. <- I did this in another script (include eventually) and it still appears month is insignificant (even with equal month weighting). Therefore, I chose to leave global centroid as is - rather than weighting it equally across months.
Test for within-group dispersion (Are monthly diets equally variable?)Show PERMDISP (permutest) results
| Term | Df | Sum_Sq | Mean_Sq | F | N.Perm | Pr(>F) |
|---|---|---|---|---|---|---|
| Groups | 2 | 1.380 | 0.690 | 0.864 | 999 | 0.411 |
| Residuals | 17 | 13.577 | 0.799 | NA | NA | NA |
Show Tukey HSD pairwise comparisons
| Group_Comparison | diff | lwr | upr | p_adj |
|---|---|---|---|---|
| July-August | -0.508 | -2.001 | 0.985 | 0.664 |
| September-August | 0.320 | -0.844 | 1.483 | 0.764 |
| September-July | 0.828 | -0.793 | 2.449 | 0.409 |
| Category | Sample1 | Sample2 | Aitchison_Distance |
|---|---|---|---|
| Most Similar | mn24_255F1 | mn24_194F2 | 0.77 |
| Most Similar | mn24_194F2 | mn24_255F1 | 0.77 |
| Most Similar | mn24_235F5 | mn24_235F4 | 0.82 |
| Most Different | mn24_261F2 | mn24_234F2 | 7.17 |
| Most Different | mn24_234F2 | mn24_261F2 | 7.17 |
| Most Different | mn24_261F2 | mn24_194F1 | 6.69 |
Figure 11. Heat map showing pairwise Aitchison distances among humpback whale fecal samples, illustrating compositional dissimilarity in prey DNA profiles.
| SampleID | LAT | LON | Location | Date | Month | WhaleID | Richness | DietEvenness | DistToCentroid |
|---|---|---|---|---|---|---|---|---|---|
| mn24_194F1 | 48.516 | -124.909 | Bamfield | 2024-07-12 | July | No | 7 | 0.612 | 3.049 |
| mn24_194F2 | 48.514 | -124.918 | Bamfield | 2024-07-12 | July | No | 3 | 0.612 | 2.317 |
| mn24_199F1 | 48.523 | -124.892 | Bamfield | 2024-07-17 | July | No | 4 | 0.793 | 1.870 |
| mn24_228F1 | 48.513 | -124.853 | Bamfield | 2024-08-15 | August | No | 18 | 0.607 | 2.638 |
| mn24_234F1 | 48.541 | -124.817 | Bamfield | 2024-08-21 | August | No | 18 | 0.761 | 2.445 |
| mn24_234F2 | 48.540 | -124.812 | Bamfield | 2024-08-21 | August | No | 2 | 0.977 | 3.507 |
| mn24_234F3 | 48.532 | -124.820 | Bamfield | 2024-08-21 | August | No | 17 | 0.292 | 2.102 |
| mn24_235F1 | 48.510 | -124.870 | Bamfield | 2024-08-22 | August | No | 18 | 0.432 | 2.086 |
| mn24_235F2 | 48.511 | -124.872 | Bamfield | 2024-08-22 | August | No | 18 | 0.637 | 2.604 |
| mn24_235F3 | 48.511 | -124.871 | Bamfield | 2024-08-22 | August | No | 4 | 0.308 | 2.567 |
| mn24_235F4 | 48.511 | -124.873 | Bamfield | 2024-08-22 | August | No | 3 | 0.850 | 2.143 |
| mn24_235F5 | 48.510 | -124.865 | Bamfield | 2024-08-22 | August | No | 3 | 0.872 | 2.087 |
| mn24_236F1 | 48.502 | -124.858 | Bamfield | 2024-08-23 | August | No | 12 | 0.413 | 1.585 |
| mn24_236F2 | 48.529 | -124.877 | Bamfield | 2024-08-23 | August | No | 17 | 0.506 | 3.252 |
| mn24_255F1 | 48.522 | -124.852 | Bamfield | 2024-09-11 | September | No | 2 | 0.772 | 2.672 |
| mn24_256F1 | 48.533 | -124.852 | Bamfield | 2024-09-12 | September | No | 1 | NaN | 3.405 |
| mn24_258F1 | 48.516 | -124.832 | Bamfield | 2024-09-14 | September | No | 3 | 0.379 | 3.019 |
| mn24_261F1 | 48.528 | -124.810 | Bamfield | 2024-09-17 | September | No | 17 | 0.271 | 2.864 |
| mn24_261F2 | 48.522 | -124.795 | Bamfield | 2024-09-17 | September | No | 1 | NaN | 4.256 |
| mn24_263F1 | 48.505 | -124.853 | Bamfield | 2024-09-19 | September | No | 18 | 0.733 | 2.944 |
Individual fecal samples exhibited low mean dietary evenness (mean = 0.6014725 ± 0.2168613), indicating that humpbacks fed predominantly on one or two prey taxa per feeding event.
Could be worth using a weighted centroid due to uneven group (month) sizes!!!