Climate-assisted persistence of tropical fish vagrants in temperate marine ecosystems
Population genomicsDNA was sourced from fin clips or gill tissue sampled from 223 individuals of Siganus fuscescens from 2013 to 2017. From the northwest to the southwest of Australia, 40 individuals were sampled from the Kimberley, 36 from the Pilbara, nine from Exmouth Gulf, seven from Coral Bay, 40 from Shark Bay, 51 from Cockburn Sound, and 40 from Wanneroo Reef (Supplementary Data 3). However, following quality filtering of these DNA sequences, three rabbitfish individuals were excluded (see below), resulting in 220 rabbitfish individuals used in all remaining analyses (Supplementary Table S4). These tissue samples were extracted using the DNeasy Blood & Tissue Kit (Qiagen, Germany) based on a modified protocol, which included an in-house binding buffer, 1.4× volume of both wash buffers, and a partial automation of the extractions on a QIAcube (Qiagen) platform to minimize human handling and cross-contamination.SNP genotyping was conducted using the DArTseq protocol at the Diversity Arrays Technology (University of Canberra, Australia), which is a reduced representation genomic library preparation method that uses two restriction enzymes46,47. Genomic DNA was digested with the enzymes PstI–SphI and PstI–NspI and small fragments (0.75) or rare (allele frequency 1% and those 620. OTUs not assigned to bacterial or eukaryotic kingdoms were removed from the dataset and the accuracy of taxonomic assignment was assessed through the use of Australian databases for marine flora and diatoms25,26. This resulted in a table containing 86 OTUs, but we only retained OTUs with at least 10 read sequences given that these are less likely to be erroneous sequences that can arise from index-tag jumping. These 78 OTUs—used in downstream statistical analyses—corresponded to cyanobacteria (Cyanophyceae), unknown Eukaryota, dinoflagellates (Dinophyceae), diatoms (Coscinodiscophyceae and Fragilariophyceae), microalgae (microscopic algae of cell size ≤20 µm including Cryptophyceae, Haptophyceae, Mediophyceae, and Chlorarachniophyceae), green macroalgae (Chlorophyta with cell size >20 µm), red macroalgae (Rhodophyta with cell size >20 µm), and brown macroalgae (Ochrophyta with cell size >20 µm) and were represented by silhouettes from PhyloPic (http://phylopic.org/about/) on Figs. 4 and 5, and Supplementary Fig. S2. We then calculated the relative abundance of the 78 OTUs (based on the total number of sequence reads from each individual stomach content, which was visualized in the figure) using a circular plot that was generated with the R-package Circlize57. We also represented the 30% most abundant OTUs across all stomach content samples with a heatmap using a Bray–Curtis distance matrix, which was computed with the R-package phyloseq73 (Supplementary Fig. S2).To investigate differences in stomach contents between tropical residents and vagrants to temperate environments, we performed a non-metric multidimensional scaling ordination (nMDS) in two dimensions based on the Bray–Curtis dissimilarity of individuals. The nMDS plot, whose stress value was 0.12, was plotted using the R-package ggplot274. To further test the dissimilarity in diet composition among tropical residents and temperate vagrants, a permutational analysis of variance (PERMANOVA) was conducted on the same distance matrix with 100,000 permutations. We also tested the homogeneity of group dispersions using the PERMDISP2 procedure with 100,000 permutations as well. The nMDS plot, PERMANOVA, and PERMDISP2 were done with the R-package Vegan60. Finally, to highlight food sources that were unique or significantly associated to a single region or a combination of regions, we used the indicator species (IndVal) analysis in the R-package Indicspecies75 with 100,000 permutations and a significance level corrected with the Benjamini and Hochberg (BH) method76 (Supplementary Data 1 and 2). Significant results were illustrated using colored Venn diagrams on Fig. 5.The 23S rRNA sequence of the kelp species, Ecklonia radiata, from the Western Australian region was not available in the NCBI database, and so three samples were collected in November 2018 at Dunsborough (southwest Australia) and their DNA was extracted with the Miniplant Kit (Qiagen) according to manufacturer’s instructions. Prior to extraction, kelp tissues were rinsed with a continuous flow of tap water for 30 min, then soaked in a solution of 70% ethanol, and finally thoroughly rinsed with Milli-Q water. Tissues were also bead-bashed twice with the Tissue Lyzer II (Qiagen) for 30 s on each cycle. The optimal yield of template DNA was estimated with qPCR following the same method as described above. Each kelp sample was prepared for single‐step fusion‐tag library build using unique index tags following the methods of DiBattista et al.77 and pooled to form an equimolar library. Size selection was also conducted with a Pippin Prep instrument using the same size range as above, and cleaning was done with QIAQuick PCR purification kit (Qiagen). Final libraries were quantified using a Qubit 4.0 Fluorometer (Invitrogen) and sequenced on the Illumina Miseq platform using 500 cycles and V2 chemistry (for paired-end sequencing).Paired-end reads were stitched together using the Illumina Miseq analysis software (MiSeq Reporter V. 2.5) under the default settings. Sequences were assigned to samples using MID tag combinations in Geneious v.10.2.6 and reads strictly matching the MID tags, sequencing adapters, and template-specific primers were retained. Each of the three samples was dereplicated into unique sequences. The unique sequence with the highest number of reads (86,000–120,000) was identical in the three samples, and it did not match any 23S rRNA gene sequences available in the NCBI database based on BLASTn. This sequence was thus designated the 23S rRNA voucher sequence of Ecklonia radiata from southwestern Australia, blasted against all OTUs found in the stomach of rabbitfish individuals in this study, and deposited on GenBank (accession number MW752516).Past and current observations, and climate modelsHistorical sea surface temperature (SST) data were acquired from two sources, each with different temporal coverage and spatial resolution. The present-day (2008–2017) and 1900–1909 SST climatologies were calculated from HadISST78, which is resolved monthly and at 1° spatially. Additionally, the National Oceanic and Atmospheric Administration (NOAA) Coral Reef Watch “CoralTemp v1.0” (daily and 5-km resolution)79 was used to assess SST anomalies during the 2011 marine heatwave.Historical and projected SST data were extracted from outputs of a suite of Coupled Model Intercomparison Project Phase 5 (CMIP5) models. We used the monthly-resolution SST model outputs that included historical greenhouse gas (Historical GHG), and representative concentration pathways of 4.5 and 8.5 W m−2 forcings (“RCP4.5” and “RCP8.5”) runs of the r1i1p1 (designation of initial conditions) ensemble member80. These models included ACCESS, CanESM, CMCC, CNRM, CSIRO, GFDL, GISS-E2-H, INMCM, MIROC, MRI, and NorESM80. The model SST data for each run (historical GHG, RCP4.5, and RCP8.5) were converted to anomalies relative to a 2008–2017 base period, and these anomalies were added to the HadISST 2008–2017 climatology. This analysis was conducted separately for both mean annual and minimum monthly mean (MiMM). Finally, we calculated ensemble means by averaging the SST anomalies from the 11 models. Ensemble means are plotted in Fig. 1 as decadal averages (thick lines) and decadal ranges (shading) of the mean annual 20 °C contour and the MiMM 17 °C contour. The historical GHG run is used to compare the observed and GHG-forced rates of warming between 1900–1909 and 2018–2017, while the two RCP runs are used to project future (2090–2099) SST scenarios. The observed 1900–1909 contours (from HadISST) fall within the ranges of those from the CMIP5 historical GHG ensembles, indicating that anthropogenic emissions are responsible for warming in this region over the past century.Surface ocean currents during the 2011 heatwave were assessed using Simple Ocean Data Assimilation (SODA) v.3.3.181, a state-of-the-art ocean model constrained by observations when and where they are available. We calculated the near-surface (0–25 m) current anomalies (relative to 1980–2015 mean) for the austral summer (January, February, March, or “JFM”) of 2011, which was the peak of the 2010–2011 Western Australia marine heatwave7. These current anomalies are plotted on top of SST anomalies in Fig. 1b. All climate analyses were performed in MATLAB2012b.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
