Genomics of rafting crustaceans reveals adaptation to climate change in tropical oceans

Domain: neuroscience
Published: 2026-02-06
Editorial outcome: minor_revision

Abstract
Rafting dispersal has been proposed as a way for coastal species to track climate-driven niche shifts. However, little information exists on how rafting species disperse and adapt to shifting environmental conditions, particularly ocean currents and salinity. Here, we integrate dispersal simulations, ecological genomics, and salinity stress experiments to investigate rafting dynamics and adaptive shifts in widely distributed crustaceans across the Indo-Australian Archipelago. We develop a quantified model to examine asymmetric gene flow between populations driven by recent seasonal oceanographic shifts. Our climatic and dispersal models suggest that rafting populations must cope with increasing salinity fluctuation caused by rapidly-shifting oceanic connectivity patterns. Our genomic data provide evidence for recent selective sweeps at osmoregulatory loci, and key duplications at glycoside hydrolase gene families. Our experimental data reveal plastic expression of osmoregulatory genes required for survival during long-distance rafting voyages. These synergies between rafting dispersal and genomic change highlight the potential for rafting species to adapt to rapidly shifting oceanographic conditions.

Introduction
Rafting dispersal is thought to play a fundamental role in determining distributions of diverse marine and coastal organisms, and can facilitate direct responses to ecological disturbance and environmental change1–3. Such long-distance voyages involve fluctuating climatic and oceanographic conditions and can potentially promote both biogeographic and adaptive shifts in rafting taxa4,5. However, opportunities to understand rafting dispersal and its role in climate-driven biological shifts have been scarce, because such oceanographic events and associated genomic shifts remain challenging to directly observe and quantify6.

For species hitchhiking on rafting substrates, ocean currents can act as vehicles for directional dispersal, potentially leading to asymmetric gene flow and source-sink biodiversity patterns7–10. Quantifying such dispersal dynamics is crucial for understanding how oceanographic processes can influence the directionality and magnitude of gene flow. In many cases, long-distance rafting species must cope with thermal and salinity shifts occurring during voyages, and obtain energy from rafting substrates for months4,5. In addition, under global change, shifting oceanographic conditions (accelerating currents; fluctuating salinity) present additional physical and metabolic challenges11–14. Thus, there is an urgent need to address how rafting species disperse, adapt, and respond to shifting climates. Here, we address these questions by integrating oceanographic, physiological, and ecological genomic analyses.

Intertidal talitrid amphipods, which are broadly distributed in the tropical Indo-Australian Archipelago (IAA), are sensitive to oceanographic change because they lack autonomous larval dispersal capacity and instead disperse primarily as adults with driftwood15. These taxa have been inferred to disperse passively via ocean currents for several thousands of kilometers1,16, and during such journeys must cope with drastic salinity shifts (15 to 34 ppt). In addition, ocean currents in the IAA connecting the South China Sea and Indonesian waters show seasonal variation, shifting northward in summer and southward in winter17,18, which have potential to impact gene flow between populations19,20. Therefore, wide-distributed intertidal talitrids in this region present a strong system for assessing adaptive shifts in response to rapid environmental change.

In this study, we examine 203 talitrid samples across the IAA and sequence the genome of the widely distributed talitrid Talorchestia martensii. We also resequence the whole genome of 104 individuals to infer their population structure and demographic histories. Then, we integrate gene flow patterns and genetic simulation models to infer their potential dispersal patterns between different geographical regions. We develop and apply oceanographic models to further test the role of ocean currents in rafting dispersal. In addition, we incorporate genomic variants into environmental modeling to assess the relative maladaptation of different populations and predict dispersal dynamics under future oceanographic forecasts. To assess genetic responses to rafting dispersal, we experimentally quantify genomic and transcriptomic shifts linked to altered salinity and identify the duplication events of genes involved in wood-feeding. Overall, our results provide a model system for elucidating rafting dispersal driven by ocean currents and shed light on the genetic adaptation mechanisms that can underpin rafting dispersal.

Results
Population dynamics of intertidal crustaceans in Indo-Australian Archipelago
A total of 203 talitrid samples from 80 localities were examined, encompassing a substantial portion of the IAA coastline (Fig. 1a, Supplementary Fig. 1a and Supplementary Data 1). Based on the monophyletic relationship and COI distance threshold delimitation (16%)21, we identified T. martensii as the most abundant and widely distributed species (154 samples; Supplementary Fig. 1b, c). We first assembled a 647.1-Mb contig-level T. martensii genome and resequenced whole genomes from 104 samples to a mean genome-wide depth of coverage of 15× (range: 10.2–47.8×) per sample (Supplementary Table 1). Single-nucleotide polymorphism (SNP) calling and stringent quality filtering revealed 6,690,030 whole-genome SNPs. Based on retained SNPs after stringent quality-filtering, phylogenetic analysis, ADMIXTURE analyses (K=5), and principal component analysis (PCA) revealed five geographical populations: Vietnam-Thailand (VT), North Borneo (NB), East Borneo (EB), South Borneo (SB), and South Sulawesi (SS) (Fig. 1b, c, and Supplementary Fig. 2).

Divergence times and long-term effective population sizes for T. martensii were estimated by Pairwise Sequentially Markovian Coalescent (PSMC) method22 and G-PhoCS23 using SNP data (Supplementary Figs. 3 and 4 and Supplementary Table 2). Our results suggest T. martensii originated ca. 2.5 million years ago (Ma) with a large effective population size (Ne), and then split into two major clades with smaller Ne (Fig. 1d and Supplementary Fig. 4a). One of these clades is distributed around South Sulawesi and East Borneo, and the other in Vietnam-Thailand, South Borneo, and North Borneo.

Genetic simulations indicate asymmetric northward gene flow
D-statistics revealed substantial gene flow between South Sulawesi and Vietnam-Thailand populations (Supplementary Table 3). The G-PhoCS23 and TreeMix24 indicated that northward gene flow (SS to VT) was stronger than southward gene flow (VT to SS) (Mann-Whitney U test, P=0.024; Fig. 1e, f). To further infer detailed population migration and gene flow patterns, we performed genetic simulations twice using Fastsimcoal2.725 with 25 simulation models in five populations (K=5), or 20 simulation models in three populations (K=3) (Supplementary Figs. 5 and 6 and Supplementary Table 4). Both simulations suggest asymmetric gene flow, with a northward dispersal model favored over random and southward dispersal models (Fig. 1g–i). The best northward model suggested that talitrids dispersed from the southern region (East Borneo and South Sulawesi) to the northern region (Vietnam-Thailand) through the central region (South Borneo and North Borneo) in a stepping-stone manner. In addition, a Mantel test performed on the pairwise geographical and genetic distances rejected isolation by geographical distance (Spearman’s P-value=0.5, Fig. 1g).

Quantified model suggests association of summer ocean currents with northward rafting
To test relationships between gene flow and ocean currents, we developed a novel model to quantify ocean currents between two sites based on the least-cost dispersal path (Fig. 2a). The direction and velocity of ocean currents play important roles in explaining gene flow between populations26. We defined direction deviation (ΔD) as the angle between the direction of ocean currents and the least-cost path line. The ΔD was calculated as the average of each grid cell along the least-cost path line, and the larger value represents the mismatch between directions of ocean currents and potential dispersal direction. The dispersal velocity (Vds) was the average value of the grid’s effective dispersal velocity, which was calculated as the projection of the ocean current velocity (Voc) onto the direction of the path. The transit time between sites was estimated by dividing the least-cost distance by dispersal velocity (Fig. 2a).

We tested the relationship between these quantified parameters and gene flow by correlation tests based on the following assumptions: (1) concordant dispersal direction with ocean currents (smaller ΔD) will contribute to higher gene flow; (2) faster dispersal velocity will contribute to higher gene flow, showing a positive relationship between velocity and gene flow; and (3) reduced transit time between two sites will contribute to higher gene flow (Fig. 2a).

As seasonal monsoons in IAA shift ocean current direction (northward flow from April to September, southward flow from October to March), correlation tests were conducted between gene flow and ocean current dynamics for each month from 2014–2016. The results suggest an opposite correlation between gene flow and ocean currents in winter and summer (Fig. 2b), in detail, gene flow was significant negative related to ocean current directions (P=0.02) and positive to ocean current velocity (P=0.03) in June (Fig. 2b, c). Due to the difference between the least-cost path and actual dispersal path, we simultaneously calculated transit times under two scenarios using equations in Fig. 2a and Lagrangian particle trajectory model27 (Fig. 2d, e), and tested the correlation between transit time and gene flow. Both analyses yielded negative correlations (Fig. 2c and Supplementary Fig. 7).

To visualize the role of ocean currents in species dispersal and distribution, we simulated the raft trajectories of T. martensii by an individual-based Lagrangian particle trajectory model, with two drift simulations over eight months, one starting in April (summer drift) and the other in October (winter drift) (Supplementary Figs. 8 and 9). According to the results from models that integrated ocean currents and environmental temperature and salinity, intraspecific connectivity occurred with both winter and summer simulations. Nonetheless, the summer drift simulations (Fig. 2d) corresponded more closely to the contemporary distribution of T. martensii than did the winter simulations (Fig. 2e). For example, 14.6% of particles from SS reached VT in summer simulations, whereas only 0.7% reached SS from VT in winter simulations (Supplementary Table 5).

To investigate the effects of driftwood mass and shape, we varied the value of “wind_drift_factor” in OpenDrift28. The results indicate that very large pieces of driftwood are likely to become stranded near their source locations, whereas smaller wind drift factors (e.g., 1%) yielded summer dispersal patterns comparable to the contemporary distribution of T. martensii (Supplementary Figs. 10 and 11 and Supplementary Tables 6 and 7).

When we simulated the rafting trajectories of each population, we found that the VT population could reach SS within eight months regardless of release timing, whereas SS could only disperse northward to VT during summer (Supplementary Figs. 8 and 9a–c). From the dispersal trajectories and quantification of dispersal percentage between populations, we found the EB population first connected with South Sulawesi (30.1% reached) and then reached SB (20.2% reached) (Supplementary Table 5), consistent with the close phylogenetic relationship between EB and SS. The quantified ocean current dynamics also showed that transoceanic dispersal from EB to SS owns consistent dispersal direction and higher dispersal velocity compared with land-connected populations (EB and SB) (Supplementary Fig. 9g).

Oceanographic change increases genetic offset
To test how rafting crustaceans will respond to future oceanographic change (in the absence of plasticity), we calculated local, forward and reverse genetic offsets29. Specifically, local offset represents the mismatch between current and predicted future genomic variation in the absence of migration. By contrast, forward genetic offset represents the minimum genetic change a population would experience if migrating to a relatively favorable future location, with high values indicating strong mismatches between contemporary and future genetic compositions. Thirdly, reverse genetic offset seeks to assess whether the genotypes required to adapt to future environments already exist within contemporary populations. To this end, we extracted a total of 117,113 SNPs with minor allele frequency > 10%, to predict genetic offset using the gradientForest (GF) approach, a regression tree-based machine-learning algorithm30.

The ranked importance of environmental variables based on GF analysis indicated that genomic variation in T. martensii was strongly associated with sea surface salinity, which shows increasing fluctuation with climate warming (Supplementary Fig. 12). In our analysis, we found that populations in Souh Sulawesi, south Borneo, the coast along the Gulf of Thailand, and many of the areas along the dispersal route exhibited higher local, lower forward, and higher reverse genetic offsets, except the northern Vietnam range (Fig. 3a–c, Supplementary Fig. 13 and Supplementary Table 8, Supplementary Data 2–4). By analyzing the dispersal potential and forward genetic offsets, we found that increased dispersal distance could mitigate large-scale maladaptation, as most populations will experience reduced habitat suitability if dispersal distances decrease to 100km (Supplementary Fig. 13e–g and Supplementary Data 3). According to the prediction with unlimited dispersal (Fig. 3b), the future suitable niche for 98.55% of populations is located north of their current distributions (a total of 22,640 grids; Supplementary Data 4).

Northward rafting supports habitat tracking but dispersal deviations increase
Based on the GF results, we integrated forward genetic offset with ocean current modeling to calculate directional deviation (ΔD) between ocean currents and least-cost dispersal paths, to assess rafting potential of T. martensii to suitable habitats. In the absence of land barriers, the least-cost path was calculated as the Euclidean distance between two sites. By contrast, when this path was interrupted by land barriers, the least-cost path was inferred from oceanographic connectivity. Based on annual mean ocean current data under 2100 SSP5-8.5 scenarios, we found that 86.95% of populations will face dispersal-direction deviations (ΔD) exceeding 90°, preventing them from reaching suitable future habitats (Fig. 3d, Supplementary Fig. 14 and Supplementary Data 5). Given that currents in the IAA are predominantly southward, but with strong seasonal variation (northward flow in summer and southward flow in winter), we further tested the species’ dispersal potential using seasonal northward ocean current data under present, 2050, and 2100 SSP5-8.5 scenarios. The results suggest that contemporary conditions facilitate northward dispersal in 81% of populations (ΔD<90°) (Fig. 3e, Supplementary Fig. 15 and Supplementary Data 6). However, compared to present ocean current conditions, 84% of populations under the 2050 SSP5-8.5 scenario show increased ΔD (63% with slower dispersal velocity and longer transit time), while 81% of populations under the 2100 SSP5-8.5 scenario show similar trends (44% with slower dispersal velocity and with longer transit time) (Fig. 3f and Supplementary Fig. 16).

Genetic changes associated with salinity shifts
Salinity shifts pose osmotic challenges to rafting species31,32. Accordingly, we found sea surface salinity to be strongly associated with genomic variation in T. martensii, with a total of 575 significant salinity-associated SNPs, located within 459 genes, consistently identified across multivariate latent-factor mixed model (LFMM)33, univariate LFMM, and redundancy analysis (RDA)34 (Fig. 4a and Supplementary Fig. 17a). Enriched genes included those involved in plasma membrane (77 genes, GO:0005886), ion transport (five genes, GO:0006811) and SLC-mediated transmembrane transport (19 genes, R-DME-425407) (Supplementary Table 9).

We then performed selective sweep analysis to detect selective signals between the VT population, which experienced the lowest salinity, and the SS population, which faced the highest salinity (Supplementary Fig. 17b). Based on the top 5% FST and θπ ratio, we identified 126 genes with selective signals in the VT group, and 53 genes in the SS group (Fig. 4b). Many of these genes are similarly implicated in osmoregulatory pathways, e.g., plasma membrane (10 genes, GO:0005886); transmembrane transport (four genes, GO:0055085); potassium ion transport (two genes, GO:0006813) (Supplementary Table 10). For example, GRIK2, a glutamate-gated ion channel gene that is known to maintain ionic homeostasis in euryhaline fish35, showed striking associations with salinity (P=3.164e−5) (Fig. 4c and Supplementary Fig. 18a). We used BLASTp to detect 1:1 orthologs in the genomes of 11 additional amphipod taxa, encompassing a variety of freshwater, marine and terrestrial habitats. Phylogenetic analysis revealed that these GRIK2 sequences clustered strongly by habitat type (freshwater; marine; terrestrial) rather than reflecting the underlying species phylogeny (Supplementary Fig. 18b), providing further support for this gene’s role in amphipod salinity adaptation.

Plastic expression of osmoregulatory genes
To further detect the expression patterns related to salinity shifts, we designed salinity stress experiments for T. martensii, and sequenced transcriptomes from the whole body, gills and legs across four salinity gradients (5, 10, 20, and untreated 35 ppt), and five duration gradients (24, 48, 72h, 15 days, and 30 days) (Fig. 4d). Measurements of hemolymph osmolality reveal significant osmolality differences between salinity gradients (Fig. 4e).

The expression cluster analysis suggests a clear divergence between gills and legs (Supplementary Fig. 19a). According to the differentially expressed genes (DEGs) and gene enrichment analysis between 5 ppt and 35 ppt treatments, the top 10 GO pathways in gills are primarily related to osmoregulation (Fig. 4f and Supplementary Tables 11, 12). Notably, pathways such as amino acid transmembrane transport (five genes, GO:0003333), potassium ion transport (four genes, GO:0006813), and cell volume homeostasis (three genes, GO:0006884) were differentially expressed exclusively in gills rather than legs. Thus, we propose that amphipod gills may be crucial osmoregulatory organs, with a comparable differential expression pattern between tissues also detected in the terrestrial talitrid, M. aosen (Supplementary Fig. 19b).

Comparisons between salinity treatments revealed that the number of DEGs increased as treatment duration increased from 24 h to 15 d and 30 d, and also as salinity decreased to 10 and 5 ppt (Fig. 4g). These dynamic expression changes highlight transcriptional plasticity36 (Fig. 4g–l and Supplementary Fig. 19c–f), with the stronger shifts linked to the larger salinity shifts (Fig. 4k) and longer treatments (Fig. 4l).

To further test for salinity plasticity relevant to long-lived natural rafting voyages, we undertook extended stress experiments lasting for 15 and 30 days. Relatively strong expression shifts were detected at low salinities (Fig. 4h) and over the longest treatments (Fig. 4i). Moreover, the number of genes enriched in ion channel transport pathways increased as treatment time increased (Fig. 4j, Supplementary Table 13). Specifically, the V-type H+ ATP genes related to proton transport showed decreased expression with increasing treatment duration (VATA, VATB, VATE, VATG, VATO). Notably, genes involved in the transport of small molecules (Tret1, ABCA3, HXT1) and inorganic ions (Sodium/hydrogen exchanger, Sodium/Potassium/Chloride Transporter and KCAB2) showed increased expression in these treatments (Fig. 4m, Supplementary Fig. 20 and Supplementary Table 14).

Genomic adaptations to driftwood feeding
To test for additional genomic adaptations enabling long-distance rafting events on driftwood, we focused on the lignocellulose digestion process (from cellulose to glucose37) (Fig. 5a) that is hypothesized to be essential for wood-feeding taxa. We first constructed a phylogenomic tree for six rafting and three non-rafting amphipod species (Fig. 5b), and identified genes involved in lignocellulose degradation mediated by Carbohydrate-Active Enzymes (CAZymes)38 including glycosyl hydrolases (GHs), carbohydrate esterases, auxiliary activitie, along with associated carbohydrate-binding modules and glycosyl transferases. The CAZyme family composition was largely conserved between rafting and non-rafting species; however, the copy number of CAZyme genes varied notably (Supplementary Fig. 21). In particular, rafting species exhibited significantly higher proportions of GH genes compared to non-rafting species (Wilcoxon signed-rank test, P=0.044) (Fig. 5c and Supplementary Table 15).

Gene identification revealed 139 gene copies in the T. martensii genome that can be classified into 39 GH families. Three of GHs participate in major steps of wood digestion: (1) Glycoside hydrolase family 9 (GH9) encodes β-1,4-endoglucanases (EC 3.2.1.4) that break cellulose (crystal) bonds to produce individual cellulose polymers; (2) GH7 encodes β-1,4-exoglucanases (EC 3.2.1.91) that cleave cellulose polymers from either end to produce 2-4-mer subunits; and (3) the GH30_1 family encodes β-glucosidases (EC 3.2.1.21) that break cellobiose (glucose 2-mer) into glucose units37 (Fig. 5a).

The gene family expansion and extraction analysis suggest that the GH7 gene family show significant expansion in the rafting species compared with non-rafting species (Fig. 5b and Supplementary Data 7). Phylogenetic analysis shows that crustacean GH9 and GH30_1 gene families are closely related to insect genes, while the monophyletic crustacean clade (including nine talitrids and non-talitrid crustaceans) of GH7 is a sister clade to fungi, suggesting it likely originated from fungi via horizontal gene transfer (Fig. 5d, e).

By comparing the species tree of the nine talitrids with the gene trees of GH proteins, we identified two gene duplication models in the evolution of these genes (Fig. 5d): one within rafting talitrids (Model I), and another at the talitrid ancestor, including non-rafting talitrids (Model II). The GH30_1 gene family consists of three talitrid clades and non-talitrid relatives (Fig. 5e), with six duplication events of Model II and three of Model I. The GH9 gene clade includes three groups: two containing both talitrids and non-talitrids, and one exclusive to talitrids. Among these, two duplications belong to Model I and seven to Model II. For the GH7 family, five duplications are in Model I and four in Model II (Fig. 5e).

To detect the potential function difference between these copies, we compared gene expression, protein sequences and structure generated by gene duplication. Using Clade II of the GH7 gene family as a case study, we found some amino acid variations in different duplication events (Fig. 5e). Protein structure analyses suggest these copies share broadly similar structures (Supplementary Table 16). Moreover, transcriptome analysis revealed that GH7, GH9, and GH30_1 gene copies are highly expressed in the gut, with comparable expression levels in rafting species despite their greater gene copy numbers relative to non-rafting species (Fig. 5f and Supplementary Fig. 22).

Methods
Sample collection and genome sequencing
A total of 203 talitrid samples from 80 localities were collected from the IAA coastline from August 2017 to September 2018. We first performed species delimitation followed Liu et al. (2023)71 and inferred phylogenetic relationships based on the concatenated dataset comprising the mitochondrial (COI, COII, 16S) and nuclear (18S rRNA, 28S rRNA, H3, NaK) gene regions. All the primers were followed by Liu et al. (2023)71. Then, we collected T. martensii specimens (voucher number: T2811) from Gerigat, Sarawak, Malaysia (E111.20°, N2.04°) for genome sequencing. We isolated genomic DNA from adult talitrids using a Qiagen Blood & Cell Culture DNA Kit (QIAGEN, Hilden, Germany). We sequenced short-insert paired-end (PE, 150 bp) libraries on the Novaseq 6000 platform (Illumina). A total of 239.88 Gb clean reads (~370×) was obtained after removing low-quality reads with poor sequencing accuracy and adapters with Trimmomatic v.0.3972. Genome assembly and annotation followed the method of Liu et al. (2023)71. We further scaffolded the contig-level assembly to a pseudochromosome scale with RagTag v2.0.073 using Morinoia aosen (GCA_030386875.1)71 as a reference.

Talorchestia martensii population resequencing and SNP calling
We used T. martensii samples from 44 sampling sites for genome resequencing. The obtained clean reads (average 15×) were mapped to the reference genome using BWA v0.7.1274. The PCR duplicates of each of the samples were removed with Picard MarkDuplicates (http://broadinstitute.github.io/picard/). We used the HaplotypeCaller engine in GATK v3.7.075 to call SNPs from the BAM alignment file for each sample. The option –ERC GVCF was used to convert the variants into gVCF format, and the gVCF of each sample was merged to allow a multi-sample joint genotype procedure using the GATK GenotypeGVCFs engine. For SNPs, we further used VCFtools v0.1.1376 and BCFtools v1.977 to remove indels and keep only biallelic SNPs with the following filtering criteria: (1) minQ > 30, (2) read depth > 5, (3) max-missing counts = 5, (4) SNPs at least 5 bp away from indel regions. We further excluded highly correlated SNPs by performing linkage disequilibrium based SNP pruning in PLINK v.1.9.078 with a default pairwise r2 threshold 0.2. This threshold reduces LD-induced bias while preserving enough independent genetic markers to accurately infer population structure, an approach previously used in marine genomic studies79.

Population structure analysis within T. martensii
To characterize population structure, we first used SNP data from all samples and constructed phylogenies using RAxML v.8.2.980. Then, we used a PCA implemented in the R package snprelate81 to identify clusters of individuals. We also conducted genetic structure analysis in ADMIXTURE82, which also allowed us to infer individual ancestry proportions. We used default parameters and defined the assumed genetic clusters K as ranging from 1 to 13.

Demographic history and gene flow analysis
We inferred the T. martensii population dynamics and demographic changes over large timescales with high resolution of different samples using the Pairwise Sequentially Markovian Coalescent (PSMC) method22. We also used the population genetic modeling framework in the program G-PhoCS v1.2.323 to estimate the lineage divergence for the T. martensii clade based on SNP datasets. Details about these methods can be found in Supplementary Note 3.

We performed another two analyses to calculate the migration rates and direction between different groups. First, the migration events were inferred using TreeMix24 by assuming zero to five migration edges. Then, to quantify gene flow between different populations, we calculated D-statistics and f4-ratio statistics with D-suite v0.483. We considered gene flow to be significant with a P value < 0.05 and |Z score| > 3.

Dispersal hypothesis and gene flow model testing
We used Fastsimcoal2.725 to simulate T. martensii population migration patterns based on the following dispersal hypotheses: (1) random dispersal from one location to surrounding places without consistent directions, and (2) directional dispersal from source to sink, including (2a) northward dispersal from the southern to northern group, and (2b) southward dispersal from the northern to southern group (Supplementary Fig. 5). The SNP data were converted to site frequency spectrum (SFS) data using easySFS.py (https://github.com/isaacovercast/easySFS). We first used five populations (K=5) to infer the potential dispersal hypothesis, and listed 25 models for simulation. The model fitness was evaluated by simulating SFS for each demographic model and comparing the simulated and observed likelihoods. The best demographic model was determined based on the maximum likelihood value. During this process, we found the dispersal hypothesis and model prior could be very complex.

To simplify the model prior, we defined groups based on population structure (K=3), with the EB and SS populations designated as the south group, the SB and NB populations as the central group, and the VT population as the north group. We also listed 20 alternative models for genetic simulation. Among these models, four tested for southward dispersal, four for northward dispersal, and 12 for random dispersal between separate regions. For each model, we carried out 50 runs with 40 expectation-maximization cycles and one million coalescent simulations to simulate migration rates, Ne, and divergence time.

Quantification of ocean current dynamics
We developed a quantified model by quantifying the direction deviation, species dispersal velocity, and transit times between two sites based on the least-cost path and gridded ocean current data (1/12° horizontal resolution). The least-cost path between each pair of individuals was measured using the R package marmap84 with bathymetric data from the ETOPO1 database85, which is hosted by NOAA. We slightly adjusted some sample positions to avoid samples on land because there was limited resolution of bathymetric data in coastal waters. We identified each grid cell that intersected the least-cost path.

We calculated the direction deviation value (ΔD) between the dispersal route direction (Ddispersal) and ocean current direction (Docean) in the least-cost path from source to destination.

ΔD = |Ddispersal − Docean| (1)

ΔD is a score between 0 and 180°, where greater than 90° indicates that the ocean current direction of the grid cell is inverse to the dispersal path connecting the populations. We then used the averaged values across all grid cells of the path as direction deviations between sites.

Then, we calculated the dispersal velocities (Vds) for each grid cell according to the gridded velocities of ocean currents (Voc) and direction deviation value (ΔD) of the least-cost path:

Vds = Voc × cos(ΔD) (2)

The dispersal velocities between sites were calculated as the averaged value of all grid cells on the least-cost path. Finally, the transit times were calculated as distance/velocity between sites.

Transit times = Least cost distance / mean Vds (3)

Ocean current dispersal trajectory simulations
We used an individual-based Lagrangian particle trajectory model26 to predict the possible trajectories and stranding events of passively drifting individuals. The passive drift trajectories of individuals were computed using the OceanDrift model implemented in the Python-based Lagrangian trajectory simulation framework OPENDRIFT.

To investigate the influence of environmental variables on the rafting pathways, we used the Pelagic Egg Drift model and added the environmental temperature and salinity into the simulations. We recorded the stranded locations and calculated the percentage of particles reached each region to evaluate population connectivity quantitatively. We used modeled surface current fields and seawater temperature and salinity from Global Ocean Ensemble Physics Reanalysis (GOEPR), which was developed by the Copernicus Marine Environment Monitoring Service (CMEMS; https://marine.copernicus.eu/). As we collected most of the samples during 2017, we used GLORYS12V1 ocean current data (1/12° horizontal resolution) over three years (January 2014–December 2016) as input to account for seasonal variation in oceanographic dynamics.

While driftwood typically follows currents, the shape and mass of these particles influence their buoyancy and sensitivity to wind slippage. Therefore, to investigate the effects of driftwood mass and shape, we changed the value of “wind_drift_factor” (1%, 1.5%, 2%, 2.5%, 3%) in our passive ocean surface drift model. The ocean–wind variables were obtained from the global atmospheric model from the National Centers for Environmental Prediction (NCEP; https://www.ncep.noaa.gov/).

We ran general drift simulations using 500 particles randomly distributed within a 3 km radius around the central coordinates of each collection site. In total, 22,000 particles were released from 44 sites, representing the full extent of the current T. martensii distribution. As seasonal monsoons in IAA shift ocean current direction and talitrid have an eight-month lifecycle, in this study, we set two releasing dates as 1 April and 1 October for passive drift simulation in summer and winter, respectively. We simulated each drifting process for a conservative eight-month period and recorded the trajectories hourly for the dispersal trajectory plot.

Identification of genomic variants associated with salinity adaptation
To identify the environment-associated variants across the whole genome, we retained SNPs with a minor allele frequency > 10%, which refers to alleles that are common enough to be reliable for the detection of traits associations50. We obtained the present oceanographic data of the four annual mean environmental variables (sea surface temperature, sea surface salinity, ocean current direction, and ocean current velocity) from Bio-ORACLE (https://www.bio-oracle.org/). To constrain the species distribution near the coastal area, we added the variable “distance from shore” from the MARSPEC database (http://marspec.weebly.com) for environment analysis86. To obtain a more accurate comparison of climate influence, we used Z-Score standardization methods X’=(X-μ)/δ to scale the climate variables, while X is original value, μ is mean of the original data, δ is standard deviation of the original data, and X’ is the scaled standardized value.

Based on the SNP data and oceanographic variables, we first calculated the weighted R2 value for each climatic variable to assess its importance and selected the top variables by ranking their importance as the major environmental variables for further analysis. Then, we used the turnover function in GF models to convert complex climate data into patterns of genetic variation (allele frequencies), subsequently to select and weight these variables to best summarize genomic variation.

To identify environment-associated variants across the whole genome, we first used the multivariate LFMM2, univariate LFMM implemented in the R package LEA v3.12.287 separately to search for associations between environmental variables and allele frequencies. Based on the five ancestry clusters inferred with ADMIXTURE, we ran univariate LFMM with five latent factors to account for population structure in the genotyped data. For each environmental variable, we ran five independent MCMC runs using 5,000 iterations as burn-in followed by 10,000 iterations. P values from all five runs were then averaged for each variant and adjusted for multiple tests using a false discovery rate correction of 5%. Second, we performed RDA using the R package vegan to identify genetic variants showing an especially strong relationship with multivariate environmental axes88. Environment-associated variants were identified as significant if their loading value was greater or less than three standard deviations above or below the average value.

Genetic offset modeling under future climate scenarios
Based on the identified environment-associated SNPs, we used the gradient forest models to calculate the genetic offset and predict how allele frequencies would shift under future oceanographic change. The future climate data [2040–2050 (2050) and 2090–2100 (2100)] of the environmental variables were obtained from Bio-ORACLE. Two shared socioeconomic pathway (SSP) scenarios were applied: SSP1-2.6 and SSP5-8.5. To estimate the spatial regions where the genotype–environment relationship will be most disrupted under future oceanographic conditions, we transformed the oceanographic variables from each of the four climate scenarios (2050 SSP1-2.6, 2100 SSP1-2.6, 2050 SSP5-8.5, and 2100 SSP5-8.5) into genetic importance and calculated the Euclidean distance between the present and future genomic importance values for each grid, which served as a metric of genetic offset.

Following the methods in Gougherty et al., 202129, we integrated migration to predict potential maladaptation to future oceanographic change and calculated three different forms of genetic offset by the gradientForest approach: the local offsets for each climate scenario, and forward and reverse genetic offsets for 2100 SSP5-8.5. To assess the sensitivity of dispersal constraints, we tested how forward offset varied under migration limits of 100, 500, 1000km, and unlimited dispersal within the IAA marine ecoregion.

As the GF algorithm and forward genetic offset methods can estimate the target location of species dispersal, we integrated the quantified ocean dynamics with forward offset prediction to test whether the ocean currents can support such dispersal. Based on the future suitable location estimated by forward genetic offset with unlimited dispersal under 2100 SSP5-8.5 scenarios, we calculated ΔD, dispersal velocity, and transit times needed to that location caused by annual mean ocean currents and the northward ocean current dynamics. Finally, we compared ΔD, dispersal velocity, and transit time differences between present and future climate scenarios by subtracting present values from future values.

Selective sweep analysis between populations
We defined the salinity of seawater surrounding each collection site as their salinity exposure and chose two populations with the most dramatic salinity exposure differences to test the selection signals related to salinity adaptation. We used a sliding window method (50 kb windows with 25 kb increments) to calculate the genome-wide distribution of the θπ ratios (θπ_high salinity/θπ_low salinity) and FST values between population pairs. The windows with significantly low and high θπ ratios (lowest 5% regions and top 5% regions) and significantly high FST values (top 5% regions) were selected as regions with strong selection signals along the genome. We defined the candidate genes under selection as those within or overlapping with regions showing selective sweeps. The Gene Ontology and Reactome pathway enrichment analysis was implemented in KOBAS 3.0 using Fisher’s exact test (P<0.05) with Drosophila melanogaster as a reference89.

Salinity stress experiments and expression plasticity analysis
To test talitrids’ genetic responses to salinity shifts during rafting events, we cultured adults of the coastal talitrid Talorchestia cf. martensii and terrestrial talitrid Morinoia aosen in different salinity seawater and evaluated the expression changes of selected genes. For coastal T. cf. martensii, we set four salinity gradients (35 ppt as control, and 20, 10, 5 ppt) with three duration gradients (24, 48, and 72 hours). To investigate the long-term expression changes of T. cf. martensii, we also treated the talitrids in 20 ppt for 15 and 30 days. For terrestrial M. aosen, we set four salinity gradients (5, 10, 20, 35 ppt) with three duration gradients (24, 48, and 72h), using terrestrial treatment as control. To detect the hemolymph osmolality difference between different treatments, we measured the hemolymph osmolality with a Fiske 210 Micro-Sample Osmometer (Advanced Instruments). The hemolymph was collected from 10–15 individuals using a micro-capillary tube to extract fluid from between the second and third abdominal. All data were obtained from three replicates.

Whole body RNA from a single individual and tissue RNA from a mixture of tissues from 10–15 individuals was extracted using the RNAsimple Total RNA Kit (TIANGEN, Beijing, China). RNA-seq was performed on the Novaseq 6000 platform, and we obtained at least 6-Gb clean data from each transcriptome after removing reads containing adapter, ploy-N, and low-quality reads. All data were obtained from three replicates in each assay. The expression level (fragments per kilobase million) of each gene was calculated with RSEM v1.3.3 and log-transformed (N+1) following quantile normalization90.

Differentially expressed genes (DEGs) between different salinity and duration gradients were identified using “DESeqDataSetFromMatrix” and “DESeq” functions in the R package DESeq2 v1.24.091. The threshold value of DEGs was set to P-value < 0.05 and absolute of fold change ≥ 2. Log2(FC) ≥ 1 indicates up-regulation, and log2(FC) ≤ −1 indicates down-regulation.

The DAPC approach was used to compare the expression of DEGs using the R package adegenet92. The gene expression plasticity in response to salinity and duration gradients was quantified as the magnitude of the shift away from the reference points.

Genetic changes related to wood-feeding in T. martensii and talitrids
We searched for GH family proteins in the T. martensii genome and eight other talitrids using the Carbohydrate Active enZymes (CAZy) database (dbCAN), which uses hidden Markov models to define the signature domains for each CAZy family93, including Hyalella azteca (GCA_000764305.4)94, Parhyale hawaiensis (GCA_001587735.2)37, Hirondellea gigas (CNGBdb Project ID CNP0005374)95, Trinorchestia longiramus (GCA_006783055.1)96, Morinoia aosen (GCA_030386875.1)71, Floresorchestia mkomani (SRR23898670 under PRJNA938803)71, Cochinorchestia sp. (SRR23898667 under PRJNA938803)71, Platorchestia pacifica (SRR23898669 under PRJNA938803)71. We analyzed all predicted proteins in the genome using HMMER v3.1b297 with an E-value threshold of 1e−5.

To assess the molecular evolution of the GH gene family, we identified GH7, GH9, and GH30_1 protein sequences in all eight talitrids. To further verify the protein sequence variations, we predicted the protein structure based on the all GH7 amino acid sequences on the SWISS-MODEL online website (https://swissmodel.expasy.org/). Moreover, we obtained these protein sequences from other organisms by a BLASTp search of the GenBank non-redundant protein database using the sequences identified in the talitrid species as query sequences. We constructed a phylogenetic tree with these three proteins separately with the PROTGAMMALG substitution model using RAxML.

To compare the genomic difference between rafting and non-rafting species, we used OrthoFinder v2.5.4 to define gene families based on protein-coding genes98. Gene family expansion and contraction analyses were performed using CAFE v4.299.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript \"Genomics of rafting crustaceans in tropical oceans: adapting to marine climate change\" through ocean current dispersal simulation, population genomics, molecular evolution and salinity stress experiments proposes to evaluate the rafting dynamics on a direct-development species distributed across the Indo-Australian Archipelago. Using climatic and dispersal models, the authors predict that the species will be able to adapt to changes in salinity levels driven by shifting oceanic patterns.\nThis is a well-written and thoroughly presented study, with analyses that are mostly soundly performed. I recommend it for publication, provided the following points are clarified and addressed:\nLine 394: The manuscript does not specify how linkage disequilibrium was accounted for during the filtering process. This is a crucial detail, as linkage disequilibrium can significantly impact demographic and selection estimates. Please clarify this aspect.\nLine 494: The MAF cut-off threshold appears to be very conservative, potentially excluding important rare variants in coding regions. Could you provide a rationale for this choice and discuss its implications for your results?\nLine 516 and throughout the manuscript: Replace the term \"adaptive SNPs\" with alternatives such as \"SNPs putatively under selection\" or \"putative adaptive SNPs\". The association results do not necessarily confirm that these SNPs are adaptive. Additionally, it is unclear how \"adaptive SNPs\" were defined—did you consider the union of results from LFMM and RDA, or only the 1,180 SNPs identified as common to both methods? The latter, being more conservative, would likely be the more robust approach. Please clarify and justify this choice.\nLastly, the bioproject/biosample SRA numbers should be publicly available.", "reviewer_2": "Given my limited expertise in the broader scope of this study, my evaluations focus primarily on the application of Lagrangian particle tracking.\nOverall evaluation:\nThis study employs OpenDrift, a widely used particle tracking package, to simulate the ocean current dispersal trajectory and indicate the possible drift trajectory of T. martensii, which rafts on driftwood. The approach provides valuable insights into the potential rafting pathways of this species. However, I have one major comment that I believe warrants further discussion in the manuscript:\nCould the authors comment on the accuracy of using numerical particles without mass or volume to simulate the trajectory of driftwood, which may have varying sizes and shapes? Are there previous studies or evidence supporting the use of a 2% wind drift factor for simulating rafting driftwood with numerical particle tracking? A discussion addressing this question would strengthen the manuscript.\nBelow, I have outlined several specific comments and suggestions to improve the clarity and robustness of the study.\n(1) Line 490:\nThe sentence “…recorded the trajectories every 6 h, with a final output time step of 12 h” is somewhat unclear. Does “6 h” refer to the timestep used to move particles from current location to the next? What is the output frequency of ocean current from GLORYSV12 – is it daily mean or something else? Please clarify.\n(2) What are the definitions of particle “active” and “stranded”? In Fig. 2d and supplementary Fig. 6 & 7, particles become stranded in both the ocean and on land. Do the ocean-stranded particles represent the final locations of particles after the eight-month tracking period? Do particles strand immediately upon reaching the coastline?\n(3) In Fig. 2d, e:\n• Consider using a different color to distinguish the initial particle locations from others. In the current figures, some initial particles are obscured by the orange-colored stranded particles.\n• How were the 500 initial particles distributed? Were they uniformly seeded across a specific spatial range?\n• Do “Summer” and “Winter” refer to the averaged trajectories of particle released in April and October from 2014 to 2016?\n• The figure caption states, “particles were released hourly for eight months”, but this seems inconsistent with the main text (line 488), which indicates that particles were released only on 1 April and 1 October. Please clarify the discrepancy.\n• Additionally, please add the units of ∆D, V, and Transit time.\n(4) Supplementary Figs. 6&7:\n• In the captions, are the stranding particles shown in “Red” rather than “Orange”?\n• The sub-titles mention April and October — what year do these represent?\n• In supplementary Fig. 6d, why are the stranded particles so close to their initial release locations in the southwest of South Sulawesi? The trajectories of particles released from the south coastline of South Sulawesi appear to be missing – did these particles disperse away from their release locations, or did they remain nearby?\n• In Supplementary Fig. 7e, should “November” in the title be corrected to “October”?\n• The figure resolution could be improved in the final publication.\n(5) Least-cost dispersal path:\nI am not familiar with the methodology behind “least-cost dispersal path,” so my comments may be less informed. However, I am curious about:\n• Transit time could also be estimated from particle tracking by recording the time for particles to travel from their source to destination. How does the transit time estimated in this study differ from that derived from particle tracking?\n• Particle tracking could provide an estimate of the percentage of particles released from one location that arrive at another location. For example, the percentage of particles traveling from East Borneo to Sulawesi relative to the initially released particles could be quantified. Could this percentage be used to evaluate genetic connectivity quantitatively?\n(6) Thermal and salinity shifts:\nThe authors note that rafting species must cope with thermal and salinity shifts during their journeys. I think OpenDrift has a functionality that allows temperature and salinity values from ocean models to be interpolated to particle locations. Could the authors use this feature to provide additional information on how these environmental variables vary along the dispersal trajectories of T. martensii?", "reviewer_3": "This manuscript provides a thorough investigation of the genomics of rafting crustaceans, detailing how ocean dispersal impacts gene flow, what genes are associated with adaptation, how those genes are expressed, and how this might translate to future maladaptation. The strength of this study is the interdisciplinary and comprehensive set of methods across multiple fields to understand the genomics of a rafting species, T. martensii. The authors integrate ecological and population genetics, dispersal simulation, ocean current modeling, and salinity stress experiments. The most robust findings include the recovery of 5 populations of T. martensii in the Indo-Australian Archipelago, determination of divergence times and effective population size for this species, and identification of wood-feeding genes required for energy intake during long rafting journeys.\nWhile there is evidence supporting the above mentioned conclusions and claims in the manuscript, there are a number of other claims that could be problematic. In general, this study highlights the importance of plasticity for the species in variable salinity environments, but the genomic maladaptation framework assumes no plasticity in the prediction. This contradiction calls into question the validity of the genomic offset analysis. The salinity stress experiment and subsequent plasticity analyses also use a terrestrial species rather than the marine rafting species, calling into question the applicability of these results to rafting species. Finally, the dispersal simulation results are unclear. More detail on these issues is provided below.\nLack of framework: The claim in the abstract that this paper provides a framework for studying adaptive shifts and predictability of adaptation seems a little strong; while analyses are comprehensive, no actual framework is proposed within the text.\nConcerns with GEA analysis: The genotype-environment associations (GEAs) could also benefit from using the more up-to-date multivariate implementation of LFMM from the LEA3 package, and the authors should also note that recent work has shown RDA is not the most reliable method for GEA analysis, though much of the field still uses these methods (see Lotterhos 2023). GEA additionally assumes that each population occupies only a single environment, which is not true for highly migratory rafting species experiencing large environmental fluctuations.\nLimitations of Genomic offsets analysis: Genomic offset modeling assumes no plasticity, so the fact that significant plasticity was found under experimental conditions undermines the genomic offset results and should be discussed. Furthermore, much like GEA, genomic offset assumes that populations occupy and are locally adapted to a single environment. Again, as a rafting species, T. martensii experience large environmental fluctuations and local adaptation may be limited in this species. Finally, while forward and reverse genetic offset methods are appropriate for species that can disperse, how dispersal and genomic offset modeling are integrated using ∆D is unclear from the methods. In general, genomic offset predictions have only been empirically tested in a few systems (mostly trees), and these tests have yielded results spanning from inaccurate predictions to relatively accurate predictions. Some empirical validation is needed before using genomic offsets in this system to predict maladaptation.\nInterpretability of gene flow and dispersal modeling: The dispersal and current modeling analyses are unclear, specifically what ∆D actually means - while the authors explain that this is a dispersal direction deviation, it is unclear why dispersal direction would be different from current direction if this species truly lacks autonomous dispersal and instead disperses passively. Figure 2 shows the quantified model and particle trajectory simulations, but does not clarify how ∆D is calculated - is it averaged across all modeled cells, or specific to just the final cell? Furthermore, when doing the gene flow model testing, the authors use K=3 despite recovering K=5 in previous analyses.\nLimitations of Plasticity analysis: The salinity stress experiments use a terrestrial species from a different genus than T. martensii, which is the focal species of the rest of the manuscript. While the plasticity analyses themselves seem sound, this calls into question the applicability of these results to marine rafting species, particularly T. martensii. At a minimum the discussion should address how applicable the salinity stress experiments are to the focal species of the manuscript. Additionally, there is significant evidence that plasticity can promote acclimation, but theoretical work suggests that it can also be harmful in the long term by precluding the need for species to genetically adapt to stress (see Ghalambor et. al. 2007). The authors should caveat these results in relevant sections of the discussion.\nMissing details and code: Because the manuscript covers so many topics and fields, detail is sometimes lacking in the main text as to what some of the terms and results actually mean. For example, the distinction between forward and reverse genetic offset is unclear in the main text, and this section would benefit from definitions of these terms that are biologically relevant. This was just one example, but other result and discussion sections would also benefit from additional detail. While the text is already long, clarifying terminology and what each analysis can actually do is important given that it is unlikely that any one reader will be an expert in every field covered in the manuscript. The written methods section is detailed and thorough. However, the provided code does not include all of the analyses; scripts for the DEG analysis, DAPC, selective sweep analysis, quantification of ocean current dynamics, gene ontology enrichment analysis, LFMM, and RDA were not in the zip file - additionally, the demographic analysis and ocean current drift modeling have no provided scripts, though again, the methods are detailed in the text. Inclusion of code used to conduct these analyses is necessary to ensure the reproducibility of the work and to ensure that the analysis is sound.\nOverall: Overall, the work integrates concepts across the fields of oceanography and ecological genomics to gain a more comprehensive understanding of how rafting species disperse, adapt, and respond to shifting climates, particularly ocean currents and salinity. While some of the conclusions of the work are well supported by the analyses conducted, many others are not, namely the genetic offsets, the salinity experiments and subsequent plasticity analyses, and the dispersal modeling. In addition to potential reanalysis of these data and separation of the experimental salinity work from the population genomic and dispersal/current modeling work, this manuscript would benefit from a more nuanced discussion of the caveats and limitations of the methods used and the results found, and from more explicit descriptions and definitions of relevant terminology and analyses in the main body of the text.\nDetailed comments:\nLine 81: “Based on these high-quality SNPs” is a bit unclear, do the authors mean that they used the high quality SNPs to do the analyses they mention next?\nLine 129: Given the claim of stepping stone dispersal, it would be nice to do an analysis confirming IBD\nLine 140: include a more definitive definition of genomic offset to better characterize what these results mean - i.e. what would a higher or lower offset mean? How does gradient forest define genomic offset?\nLine 146: Clarify the difference between local, forward, and reverse offset in biologically relevant terms\nLine 161: Why was the cutoff for dispersal direction deviation 90º - what is the biological relevance here?\nLine 173: Clarify whether identified SNPs were identified as those recovered by both methods (i.e. overlapping), or the combination of all SNPs identified, including those identified by just one method - the supplemental shows that only SNPs identified by both methods were kept, but this isn’t clear in the text\nLine 176: Give a number for how many genes were enriched in each of the GO categories\nLine 202: “and so on” is non-specific, it would be good to mention why the first two pathways were described but not others - was there some cutoff for why some pathways made the main text and others didn’t?\nLine 237: Why was Paryhale hawaiensis used as a reference?\nLine 260-261: It would be helpful to quickly outline what some of the expectations of Pleistocene sea-level fluctuations and Wallace’s Line are, either here or in the introduction\nLine 292: Careful about using the term “genetic vulnerability” or “vulnerability” when talking about genetic offset, the models can describe potential maladaptation, but how that actually translates to meaningful fitness changes and population vulnerability is still an area of active research in the field\nLine 315: Plasticity does not always facilitate long term population survival even if it does facilitate acclimation and survival in fluctuating environments\nLine 362: It is unclear exactly what downstream analyses the mitochondrial fragments are used for\nLine 373: What was the sequencing depth?\nLine 374: How are “low-quality read” defined?\nLine 379: What is the threshold of nucleotide diversity used?\nLine 400: What is meant by “concatenated all SNPs” - I assume this mean SNPs for all samples, but clarifying this would be helpful\nLine 426: Each model was run “multiple times”, but what is the exact number there?\nLine 438: The authors define groups here based on K=3, but used K=5 for all other analyses and recovered K=5 in their population structure analyses\nLine 501: It is unclear to me what is meant by “transforming multidimensional climatic variables into multidimensional genetic space” - I assume this is part of variable selection, but an actual description of this method would be helpful\nLine 501: It seems like the authors did not scale the climate variables before using them in GEA/genetic offset analysis, which would allow for a more accurate comparison of their influence on the response variables within the analysis\nLine 551: Defining the salinity of seawater surrounding each collection site as their salinity tolerance when there is migration between populations seems counterintuitive, can there really be a site-specific salinity tolerance when the authors also describe significant plasticity and tolerance to large salinity shifts during dispersal?", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "The study investigates how the rafting crustacean Talorchestia martensii adapts to climate change in the Indo-Australian Archipelago. The authors combined ecological modeling and genomics to explore the dispersal patterns and adaptive mechanisms of rafting species under rapid environmental changes. From 203 samples collected across the region, five geographically distinct T. martensii populations were identified. Simulations revealed stronger northward migration patterns to shape the populations compared to random or southward dispersal. Using a quantified model and particle trajectory simulations, the authors further confirmed asymmetric northward gene flow, particularly during summer monsoons, consistent with seasonal ocean currents.\nGenomic data from the identified populations revealed that salinity fluctuations acted as a major selective pressure during rafting dispersal. The authors utilized transcriptome data from salinity stress experiments to identify genes associated with salinity adaptation, such as GRIK2 and SLC6A13. However, the lack of statistical significance tests and functional validation raises concerns regarding the robustness of these findings.\nOverall, this is a well-executed study on the adaptation of rafting crustaceans, and the manuscript should be considered for publication after addressing the following comments:\nThe manuscript applies various bioinformatics tools, which are essential for obtaining reliable results. However, listing all software without explanation in the main text is distracting. I recommend moving the software to the Methods section and briefly explaining the conceptual rationale behind each analysis in the main text.\nLine 180: Could the authors provide a list of GO terms enriched in salinity-associated genes? It is unclear whether these genes represent false positive associations or are truly enriched in certain biological functions.\nLine 186: The duration of salinity stress experiments (24, 48, and 72 hours) seems too short to reflect long-term acclimation. Given that rafting voyages can last for months, a longer experimental period would better represent natural conditions.\nLines 190-192: The manuscript highlights a few selected differentially expressed genes (DEGs). However, what is the general expression pattern for all DEGs? Cherry-picking without functional validation makes it difficult to draw robust conclusions.\nLine 215: It would be interesting to see how gene expression changes under longer treatments. This could provide insights into whether plastic responses increase over time.\nLine 218: Are there significance tests for these genes? Including statistical validation would strengthen the interpretation of results.\nLines 332-333: This section seems a bit weak to me. The discussion on wood degradation genes seems narrowly focused on a few candidates. Whole-genome comparisons between rafting and non-rafting species would help identify a broader set of key genes involved in driftwood digestion.\nThe discussion is fairly weak, I would expect to see more discussion on how biotic interactions (e.g., competition) affect rafting success; comparisons between salinity adaptation mechanisms in T. martensii and other marine species; the evolutionary history of GRIK2 and its role in other species.\nLines 348-349: Could you quantify the expression dosage of duplicated genes in T. martensii and compare it to species with fewer gene copies?\nZhiqiang Ye, Professor\nCentral China Normal University"}, "author_response": "Response: Thank you for your wonderful suggestions. We have added more details and explanations in the revised manuscript according to your suggestions. Detailed explanations can be found below.\nLine 394: The manuscript does not specify how linkage disequilibrium was accounted for during the filtering process. This is a crucial detail, as linkage disequilibrium can significantly impact demographic and selection estimates. Please clarify this aspect.\nResponse: We further excluded highly correlated SNPs by performing a linkage disequilibrium based SNP pruning in PLINK v.1.9.0 with a pairwise r2 threshold 0.2 (Purcell et al. 2007).\nWe have added description in the “Materials and methods” section, line 518 to 520, as “We further excluded highly correlated SNPs by performing a linkage disequilibrium based SNP pruning in PLINK v.1.9.072 with a pairwise r2 threshold 0.2.”\nLine 494: The MAF cut-off threshold appears to be very conservative, potentially excluding important rare variants in coding regions. Could you provide a rationale for this choice and discuss its implications for your results?\nResponse: In the context of genetic studies, minor allele frequency (MAF) > 0.1 is often used as a threshold for genetic association studies. Alleles with a frequency higher than 10% are typically refers to alleles that are common enough to be reliably studied and statistically analyzed in large populations. Such alleles are often of greater interest in association studies because they are more likely to be involved in common diseases or traits. Rare variants (with MAF < 0.01 or 1%) are often excluded from analysis because they may not have enough statistical power or may not be reliably detected in smaller sample sizes.\nWe have added more explanation on this value “Materials and methods” section, line 623 to 625, “we retained SNPs with a minor allele frequency > 10%, which refers to alleles that are common enough to be reliably for the detection of traits associations53.”\nLine 516 and throughout the manuscript: Replace the term \"adaptive SNPs\" with alternatives such as \"SNPs putatively under selection\" or \"putative adaptive SNPs\". The association results do not necessarily confirm that these SNPs are adaptive.\nResponse: We have changed the term to \"putative adaptive SNPs\" throughout the revised manuscript.\nAdditionally, it is unclear how \"adaptive SNPs\" were defined—did you consider the union of results from LFMM and RDA, or only the 1,180 SNPs identified as common to both methods? The latter, being more conservative, would likely be the more robust approach. Please clarify and justify this choice.\nResponse: The \"putative adaptive SNPs\" represents the SNPs identified as common to LFMM and RDA methods. We identified adaptive variants as the genetic variants show higher environmental associations from the intersections of different GEA (genome-environmental analysis), which would provide a more conservative and robust results for us.\nWe have reanalyzed this part according to the reviewer’s suggestions and corrected the “Results” section as “Accordingly, we found sea surface salinity to be strongly associated with genomic variation in T. martensii, with a total of 575 significant salinity-associated SNPs, located within 459 genes, consistently identified across multivariate latent-factor mixed model (LFMM)33, univariate LFMM, and redundancy analysis(RDA)34(Fig. 4a and Supplementary Fig. 17a).”\nLastly, the bioproject/biosample SRA numbers should be publicly available.\nResponse: We have submitted the sequenced data to NCBI under the accession number PRJNA1297316.\nReviewer #1 (Remarks on code availability):\nThe codes for evaluating the ocean currents dynamics are presented, bur not for the DE analyses. The latter is partially described on methods section. I would suggest that the authors add this one to the figshare repository as well.\nResponse: We have submitted 11 R, Python or Shell scripts used in this study to figshare https://figshare.com/s/e164e7f678e789439b69. We also upload a file named readme.txt to introduce each script used in which analysis. For example, “1.Population_genetics.txt”, the script was used for population genetic analysis like SNP calling, filtering, TreeMix, Admixture, Dsuite, and Selective sweep analysis. “8.DEGs.R” is the R script used for DEG analysis.\nReviewer #2 (Remarks to the Author):\nGiven my limited expertise in the broader scope of this study, my evaluations focus primarily on the application of Lagrangian particle tracking.\nOverall evaluation:\nThis study employs OpenDrift, a widely used particle tracking package, to simulate the ocean current dispersal trajectory and indicate the possible drift trajectory of T. martensii, which rafts on driftwood. The approach provides valuable insights into the potential rafting pathways of this species. However, I have one major comment that I believe warrants further discussion in the manuscript:\nCould the authors comment on the accuracy of using numerical particles without mass or volume to simulate the trajectory of driftwood, which may have varying sizes and shapes? Are there previous studies or evidence supporting the use of a 2% wind drift factor for simulating rafting driftwood with numerical particle tracking? A discussion addressing this question would strengthen the manuscript.\nResponse: Thank you for your wonderful comments. OpenDrift is designed to model ocean currents, and by using particles without mass or volume, it simplifies the model and reduces the computational complexity by focusing purely on the movement and dispersal of particles driven by ocean currents. This allows for faster simulations, making it easier to predict broad trajectories of drift over large areas.\nTo answer the questions about driftwood mass and shape, and incorporate salinity and temperature variables in the dispersal simulations, we emailed the author of the software and got some potential methods to solve these problems.\nReal driftwood does not behave like a dimensionless particle. Driftwood traits (e.g., shape, size and density) differ strongly, e.g., density determines the position of a wood fragment relative to the water surface, and specific morphological attributes (e.g., a dorsal sail) may increase mechanical friction. According to Murphy et al., 2021, wood volume, density and water density would be the main variables influencing the buoyancy of driftwood, which all influence how it moves in the water. Larger pieces of driftwood tend to stay afloat for longer, because larger pieces are more buoyant.\nIn transport models for particles floating near the surface, the surface wave effects are often included by simply superimposing their Stokes drift (the small net drift induced by waves) upon wind‐driven flows and currents. The \"wind_drift_factor\" is the fraction of the wind need to add to the surface current, to match the observed drift. It's physical meaning is thus not a direct wind drag on the object (by some called windage), but rather a compensation factor for the inability of ocean models to resolve the very strong shear observed in the upper few centimeters of the ocean (Laxague, 2018). By experiments, the \"wind_drift_factor\" is found to be around 3-3.5% of the wind (Schwartberg, 1971, Jones, 2016). The surface Stokes drift is about 1.5% of the wind. Thus when Stokes drift is included in the simulations, wind_drift_factor should be ~2%. If Stokes drift is disabled or not available, wind_drift_factor should be increased to about 3%.\nTherefore, the \"wind_drift_factor\" in the OceanDrift model could also be used and would be determined as function of the driftwoods’ shape and mass. For example, for particles floating higher with large part in air, the wind_drift_coefficients would be higher. To compare these effects, we simulated the dispersal path difference under wind_drift_factor 1%, 1.5%, 2%, 2.5%, 3%.\nTo choose the best model of wind_drift_factor, we set two criteria: (i) northward dispersal from South Sulawesi to Vietnam/Thailand should be supported, (ii) more connections between populations are present, because this species exits significant gene flow rather than independent species.\nFrom these simulations, we found larger driftwood could be easy to be stranded near their source location. For example, in the summer simulations, the South Sulawesi source particle dispersed to Vietnam-Thailand would be decreased from 0.8% to 0% as the wind factor enlarges from 1% to 3%. When the wind drift factor was setted as 3%, most particles released from VT, SB, EB and NB were stranded near their source locations. Similarly, in the winter simulations, we found the VT populations could reach SS when the wind drift factors as 1%, but when the wind drift factor increased to 2% or more, their connection would disappear and all particles released from VT stranded in the Vietnam-Thailand area. Therefore, we thought the small driftwood could build a better connection between islands of IAA.\nWith the simulation of wind_drift_factor 1%, we found the summer connections stronger than winter connections. At summer simulations, most populations show connection with other populations except Vietnam-Thailand. At winter simulations, three geographical populations (East Borneo, South Sulawesi, South Borneo) just did not show connections with other populations, and stay at their original site. Therefore, combined with the northward gene flow direction calculated from population genetics analysis, we propose the summer ocean currents is the driving force of the population connections.\nSupplementary Fig. 10. Simulations with different wind drift factors 1%, 1.5%, 2%, 2.5%, and 3% based on summer ocean currents that released in 1 April 2016 and experienced eight months drift.\nSupplementary Table 7. The dispersal ratio of populations from source to sink under Summer ocean currents simulations.\nRatio Ratio Ratio Ratio Ratio\nSource Sink\nWind1% Wind1.5% Wind2% Wind2.5% Wind3%\nVietnam/Thailand Vietnam/Thailand 100 100 99.9 100 100\nVietnam/Thailand EastBorneo 0 0 0 0 0\nVietnam/Thailand SouthSulawesi 0 0 0 0 0\nVietnam/Thailand SouthBorneo 0 0 0 0 0\nVietnam/Thailand NorthBorneo 0 0 0.1 0 0\nEastBorneo Vietnam/Thailand 6 0.4 0 0 0\nEastBorneo EastBorneo 50.3 79.1 80.1 96.9 99.5\nEastBorneo SouthSulawesi 3.7 1.1 5.1 2.3 0.4\nEastBorneo SouthBorneo 35.3 17.1 14.3 0.9 0\nEastBorneo NorthBorneo 4.7 2.3 0.5 0 0\nSouthSulawesi Vietnam/Thailand 0.8 0.2 0 0 0\nSouthSulawesi EastBorneo 7.3 16.5 17.8 18.7 17.4\nSouthSulawesi SouthSulawesi 69 64.2 68.5 65.3 65.6\nSouthSulawesi SouthBorneo 20.6 18.6 13.5 15.9 16.6\nSouthSulawesi NorthBorneo 2.2 0.5 0.1 0.1 0.4\nSouthBorneo Vietnam/Thailand 0 0 0 0 0\nSouthBorneo EastBorneo 6.3 1.5 0 0 0\nSouthBorneo SouthSulawesi 0 0 0 0 0\nSouthBorneo SouthBorneo 93.7 98.5 100 100 100\nSouthBorneo NorthBorneo 0 0 0 0 0\nNorthBorneo Vietnam/Thailand 2.2 0.8 0 0 0\nNorthBorneo EastBorneo 14.6 18.9 4.6 2.8 6.1\nNorthBorneo SouthSulawesi 0.3 0 0 0 0\nNorthBorneo SouthBorneo 2.9 0.4 0.4 0 0\nNorthBorneo NorthBorneo 80 79.9 95 97.2 93.9\nSupplementary Fig. 11. Simulations with different wind drift factors 1%, 1.5%, 2%, 2.5%, and 3% based on winter ocean currents that released in 1 October 2016 and experienced eight months drift.\nSupplementary Table 8. The dispersal ratio of populations from source to sink under winter ocean currents simulations.\nRatio Ratio Ratio Ratio Ratio\nSource Sink\nWind1% Wind1.5% Wind2% Wind2.5% Wind3%\nVietnam/Thailand Vietnam/Thailand 92.8 97.8 99.4 99.9 100\nVietnam/Thailand EastBorneo 0 0 0 0 0\nVietnam/Thailand SouthSulawesi 1.1 0.9 0.1 0 0\nVietnam/Thailand SouthBorneo 5.8 0.8 0.3 0 0\nVietnam/Thailand NorthBorneo 0.4 0.6 0.2 0 0\nEastBorneo Vietnam/Thailand 0 0 0 0 0\nEastBorneo EastBorneo 100 100 100 100 100\nEastBorneo SouthSulawesi 0 0 0 0 0\nEastBorneo SouthBorneo 0 0 0 0 0\nEastBorneo NorthBorneo 0 0 0 0 0\nSouthSulawesi Vietnam/Thailand 0 0 0 0 0\nSouthSulawesi EastBorneo 0 0 0.1 0.4 0.2\nSouthSulawesi SouthSulawesi 88.7 87.5 85.3 92.9 90.3\nSouthSulawesi SouthBorneo 11.3 12.5 14.6 6.6 9.5\nSouthSulawesi NorthBorneo 0 0 0 0 0\nSouthBorneo Vietnam/Thailand 0 0 0 0 0\nSouthBorneo EastBorneo 0 0 0 0 0\nSouthBorneo SouthSulawesi 0 0 0 0 0\nSouthBorneo SouthBorneo 100 100 100 100 100\nSouthBorneo NorthBorneo 0 0 0 0 0\nNorthBorneo Vietnam/Thailand 18.3 5.1 3.3 3.9 1.9\nNorthBorneo EastBorneo 6.1 6.2 5.7 6.7 1.8\nNorthBorneo SouthSulawesi 3.2 0.6 2.5 0.3 5.7\nNorthBorneo SouthBorneo 9.4 11.5 10.2 4.4 5.7\nNorthBorneo NorthBorneo 62.9 76.5 78.4 84.7 85\nWe have added more discussion in line 311 to 319 as “Particle trajectory simulations provide a powerful tool in rafting dispersal studies45,46. By using drifting particles that lack mass or volume, such simulations carry the advantage of reducing computational complexity. However, the hydrodynamics of real driftwood fragments may be strongly influenced by parameters such as mass, volume, density and water density47. By comparing \"wind_drift_factors\" of 1% to 3%, we found that large driftwood pieces often strand near their source locations, whereas smaller driftwood fragments more readily disperse between islands of IAA(Supplementary Figs. 10 and 11). Moreover, by integrating salinity and temperature variables with these complex particle models, our simulations revealed a substantial increase in successful rafting between SS and VT (14.6%) than did comparable winter simulations (0.7%) (Fig. 2d, e).”\nBelow, I have outlined several specific comments and suggestions to improve the clarity and robustness of the study.\n(1) Line 490:\nThe sentence “…recorded the trajectories every 6 h, with a final output time step of 12 h” is somewhat unclear. Does “6 h” refer to the timestep used to move particles from current location to the next? What is the output frequency of ocean current from GLORYSV12 – is it daily mean or something else? Please clarify.\nResponse: Our initial analysis recorded the trajectories every 6 h. However, in the newly revised manuscript, we record trajectories hourly. The ocean current data from GLORYSV12 is daily mean. When we analyzed the dispersal trajectories, the software saved data only at a specified time interval by providing time_step_output, with default value of 3600 seconds (one hour). Therefore, we can record the trajectory hourly and plot the dispersal trajectory.\nWe have changed the sentence as: “We simulated each drifting process for a conservative eight-month period and recorded the trajectories hourly for the dispersal trajectory plot.”\n(2) What are the definitions of particle “active” and “stranded”? In Fig. 2d and supplementary Fig. 6 & 7, particles become stranded in both the ocean and on land. Do the ocean-stranded particles represent the final locations of particles after the eight-month tracking period? Do particles strand immediately upon reaching the coastline?\nResponse: “active” refers to particles that remain drifting on the sea surface, whereas “stranded” accounts for particles that have come into contact with shorelines.\nWe added one code as “o.set_config('general:coastline_action', 'stranding')”, so the particles would be deactivated if they hit the coastline and exhibit as stranded immediately upon reaching the coastline. In these supplementary Figures, they look like stranded in the ocean, but in fact, they were stranded on some small islands (Figures below, including Paulu Bangka, Paulu Belitung, Kepulauan Riau).\nWe have added the explanation in the legend of Fig. 2d and e, “Simulated T. martensii stranding events (orange circles represent particles that have come into contact with shorelines) and active (blue circles are particles that remain at sea) distribution across the IAA in summer(d) and winter (e) in 2016.”\n(3) In Fig. 2d, e:\n• Consider using a different color to distinguish the initial particle locations from others. In the current figures, some initial particles are obscured by the orange-colored stranded particles.\nResponse: Thank you for your suggestions, we have changed the particle colors in Fig. 2d, e. In the new figures, purple for initial particles, orange for stranded particles, and blue for active particles.\n• How were the 500 initial particles distributed? Were they uniformly seeded across a specific spatial range?\nResponse: We have explained these particles’ distribution in the “Materials and Methods” section as: “We ran general drift simulations using 500 particles randomly distributed within a 3-km radius around the central coordinates of each collection site. In total, 22,000 particles were released from 44 sites, representing the full extent of the current T. martensii distribution.”\n• Do “Summer” and “Winter” refer to the averaged trajectories of particle released in April and October from 2014 to 2016?\nResponse: The summer means the particles were released on 1 April 2016 and experienced eight months drift. Winter means the particles were released on 1 October 2016 and experienced eight months drift. The correlation analyses indicated that gene flow and ocean characters (direction, velocity) show similar trends from 2014 to 2016 (Fig. 2b). We chose 2016 to do particle release analysis.\nWe have changed the legend in Fig. 2 as “Simulated T. martensii stranding events (orange circles represent particles that have come into contact with shorelines) and active (blue circles are particles that remain at sea) distribution across the IAA in summer(d) and winter(e) in 2016.”\n• The figure caption states, “particles were released hourly for eight months”, but this seems inconsistent with the main text (line 488), which indicates that particles were released only on 1 April and 1 October. Please clarify the discrepancy.\nResponse: We have changed the sentence as: “Particles were recorded every 6 h for eight months after simulated releases from a total of 44 sites (22,000 particles) in 1 April (d) or 1 October (e).”\n• Additionally, please add the units of ∆D, V, and Transit time.\nResponse: the units of ∆D is degree, V is m/s, and transit time is 106 seconds. We have added the units in the figures.\n(4) Supplementary Figs. 6&7:\n• In the captions, are the stranding particles shown in “Red” rather than “Orange”?\nResponse: We have changed the color and corrected the captions.\n• The sub-titles mention April and October — what year do these represent?\nResponse: The correlation analyses indicated that gene flow and ocean characters (direction, velocity) show similar trends from 2014 to 2016 (Fig. 2b). We chose 2016 for the particle release analysis.\nWe have added explanation in each figure as “All the simulations were based on ocean currents data from 2016.”\nIn supplementary Fig. 6d, why are the stranded particles so close to their initial release locations in the southwest of South Sulawesi? The trajectories of particles released from the south coastline of South Sulawesi appear to be missing – did these particles disperse away from their release locations, or did they remain nearby?\nResponse: We reanalyzed the dispersal pathways of South Sulawesi population with ocean currents and environmental variables (salinity and temperature). This simulation started on 1 October and experienced eight months of rafts. We found most particles still stranded in South Sulawesi, and some particles will move southward and still remain active without stranding (below figure panel a).\nWe thought this was determined by the winter ocean currents, which were flowed toward the island or southward, so these particles were stranded or moved southward (below figure panel b). From this perspective, the northward dispersal of South Sulawesi population is very weak in the winter, which suggest winter ocean currents is not the driving force for the species dispersal. We added explanation in the legend as “Because the winter ocean currents flowed southward, particles released from South Sulawesi were moved southward and stranded around their source.”\nFig. Stranded particles in the southwest of South Sulawesi. a, Most particles still stranded (red points) in South Sulawesi, there are also some particles will move southward and remain active (blue). b, Winter ocean currents in this region lowed toward the island or southward. \n• In Supplementary Fig. 7e, should “November” in the title be corrected to “October”?\nResponse: In the original Supplementary Fig. 7d and e, we want to express that East Borneo population first connected the South Sulawesi, regardless of summer or winter drift simulations. Therefore, we show the summer ocean currents, the EB populations connected with the SS in three months that start in April (d). While for the drift in winter, we found the particle simulations that start in October, the EB population will strand near their source site in 14 days, so the simulation will stop on October 14th. Despite it, we found the EB populations first connected the SS within one month and ten days when the simulation started in November.\nTherefore, we reorganized this as Supplementary Fig. 9 and added figure legend as “d-f, East Borneo (4 sites, 2000 particles) population across the IAA, suggesting individuals from South Borneo were first connected with South Sulawesi regardless of summer (d) or winter (e, f) drift simulations. While for the drift in winter, we found that the winter ocean currents that start in October, the EB population will be stranded near their source site in 14 days (e). If the simulation starts from November, the EB populations first connected the SS within one month and ten days (f).\n• The figure resolution could be improved in the final publication.\nResponse: We have reorganized all the figures with higher resolution.\n(5) Least-cost dispersal path:\nI am not familiar with the methodology behind “least-cost dispersal path,” so my comments may be less informed. However, I am curious about:\nResponse: The least cost dispersal path refers to the shortest route between two sites in the ocean drifting. The methods were inspired by Hodel et al., 2018, using the shortest path on the sea between two sites. The ocean current variables on the shortest path could be viewed as the environmental resistance, which may have importance for explaining geographic patterns of genetic variation (Hodel et al., 2018). From the source site to the sink site, we identified every grid cell that intersected the least-cost path line. If the ocean current bearing of a pixel were identical to the bearing of the line, the pixel would be given a score of zero. If the bearing of the pixel were different from the bearing of the line, it would be given a score between zero and 180. This is the fundamental explanation of ∆D of each grid cell. Under such a designation, a ∆D value larger than 90° precludes dispersal along the least-cost path.\nWe further explain this approach in the Discussion section, “Furthermore, the quantified ocean dynamics simplify the species dispersal as the shortest path between two sites, and calculate the ocean currents variables on this path to explain their relationship with gene flow26. The mutual corroboration among analyses suggests that quantified ocean dynamics represent a promising tool for revealing correlations between ocean currents and rafting dispersal.”\n• Transit time could also be estimated from particle tracking by recording the time for particles to travel from their source to destination. How does the transit time estimated in this study differ from that derived from particle tracking?\nResponse: In our original calculation, we defined the transit time as distance/velocity. This time may differ from transit time in the OpenDrift simulations, because our estimations were calculated from the shortest dispersal distance, which may not be the real dispersal path.\nWe also used the dispersal time estimated from OpenDrift to estimate their relationship with gene flow. Because the winter simulations failed to connect most population pairs (13/25), we instead relied on summer simulations for the correlation test. The summer simulations in 2014–2016 years suggest a negative relationship between transit time and gene flow, which means longer dispersal time would lead to lower gene flow. Such results correspond well with our quantified ocean currents method, which shows that the summer ocean currents would be the major driving force for the gene flow. The mutual corroboration among analyses suggests that quantified ocean dynamics could be a useful method to estimate the rafting dispersal.\nWe have added Results as “Due to the difference between least-cost path and actual dispersal path, we simultaneously calculated transit times under two scenarios using equations in Fig. 2a and Lagrangian particle trajectory model27(Fig. 2d,e), and tested the correlation between transit time and gene flow. Both analyses yielded negative correlations(Fig. 2c and Supplementary Fig. 7).”\nSupplementary Fig. 7 The transit times estimated by opendrift and the correlation with gene flow.\n• Particle tracking could provide an estimate of the percentage of particles released from one location that arrive at another location. For example, the percentage of particles traveling from East Borneo to Sulawesi relative to the initially released particles could be quantified. Could this percentage be used to evaluate genetic connectivity quantitatively?\nResponse: Thank you for your suggestions. We exported the source location and sink location of each particle and quantified the dispersal percentage as the ratio between particles from source to sink and total particles from source (Supplementary Tables 6-8). From the estimations added salinity and temperature, we found the connection between East Borneo and Sulawesi (30.1%) was stronger than that with South Borneo (20.2%). Therefore, we think the stronger gene flow between East Borneo and Sulawesi is reasonable.\nWe have added these Results in the revised manuscripts “From the dispersal trajectories and quantification of dispersal percentage between populations, we found the EB population first connected with South Sulawesi (30.1% reached) and then reached SB (20.2% reached) (Supplementary Fig. 9d, e, Supplementary Table 6), consistent with the close phylogenetic relationship between EB and SS.”\n(6) Thermal and salinity shifts:\nThe authors note that rafting species must cope with thermal and salinity shifts during their journeys. I think OpenDrift has a functionality that allows temperature and salinity values from ocean models to be interpolated to particle locations. Could the authors use this feature to provide additional information on how these environmental variables vary along the dispersal trajectories of T. martensii?\nResponse: Originally, we used oceandrift model in the OpenDrift software, which only estimates the dispersal pathway with ocean currents. In the revised manuscripts, we used the PelagicEggDrift model, which can combine environmental variables such as salinity and temperature in addition to ocean currents for advection. According to results from models adding temperature and salinity values, we found the general dispersal directions and dispersal ratios between populations were similar to our original simulations (Supplementary Table 6).\nChanges have been made in the section “Results”, including “According to the results from models that integrated ocean currents and environmental temperature and salinity, intraspecific connectivity occurred with both winter and summer simulations. Nonetheless, the summer drift simulations (Fig. 2d) corresponded more closely to the contemporary distribution of T. martensii than did the winter simulations (Fig. 2e). For example, 14.9% of particles from SS reached VT in summer simulations, whereas only 0.7% reached SS from VT in winter simulations (Supplementary Table 6).”\nWe also added more changes in the “Discussion” section, “Moreover, by integrating salinity and temperature variables with these complex particle models, our simulations revealed a substantial increase in successful rafting between SS and VT (14.6%) than did comparable winter simulations (0.7%) (Fig. 2d,e).”\nIn the “Materials and methods” section, we added the methods description as “To investigate the influence of environmental variables on the rafting pathways, we used the PelagicEggDrift model and added the environmental temperature and salinity into the simulations. We recorded the stranded locations and calculated the percentage of particles to each region to evaluate populations connectivity quantitatively. We used modeled surface current fields and sea water temperature and salinity from Global Ocean Ensemble Physics Reanalysis (GOEPR), which was developed by the Copernicus Marine Environment Monitoring Service (CMEMS; https://marine.copernicus.eu/).”\nSupplementary Table 6. The dispersal ratio of populations from source to sink under simulations with ocean currents, and that combined with environmental variables. oc: ocean currents; oc+ts: ocean currents+salinity+temperature.\nSummer Winter\nSource Sink\nRatio_oc Ratio_oc+ts Ratio_oc Ratio_oc+ts\nVietnam/Thailand Vietnam/Thailand 98.4 98.7 75.3 75.6\nVietnam/Thailand EastBorneo 0 0 0 0\nVietnam/Thailand SouthSulawesi 0.2 0 1.4 0.7\nVietnam/Thailand SouthBorneo 0 0 12 4.8\nVietnam/Thailand NorthBorneo 1.4 1.2 11.2 19\nEastBorneo Vietnam/Thailand 0 0 0 0\nEastBorneo EastBorneo 46.5 49.7 100 100\nEastBorneo SouthSulawesi 23.9 30.1 0 0\nEastBorneo SouthBorneo 29.6 20.2 0 0\nEastBorneo NorthBorneo 0 0 0 0\nSouthSulawesi Vietnam/Thailand 2 14.6 0 0\nSouthSulawesi EastBorneo 0.3 0.7 0 0\nSouthSulawesi SouthSulawesi 31.1 22.6 90.3 86.6\nSouthSulawesi SouthBorneo 58.1 54.4 9.7 13.4\nSouthSulawesi NorthBorneo 8.5 7.7 0 0\nSouthBorneo Vietnam/Thailand 0 0 0 0\nSouthBorneo EastBorneo 0 0 0 0\nSouthBorneo SouthSulawesi 0 0 0 0.7\nSouthBorneo SouthBorneo 100 100 100 99.3\nSouthBorneo NorthBorneo 0 0 0 0\nNorthBorneo Vietnam/Thailand 6.8 9 0 0\nNorthBorneo EastBorneo 21.1 15.7 3.6 3.6\nNorthBorneo SouthSulawesi 1.4 0.8 0.4 0.4\nNorthBorneo SouthBorneo 0.1 0 0 0\nNorthBorneo NorthBorneo 70.6 74.5 96 96\nReviewer #3 (Remarks to the Author):\nThis manuscript provides a thorough investigation of the genomics of rafting crustaceans, detailing how ocean dispersal impacts gene flow, what genes are associated with adaptation, how those genes are expressed, and how this might translate to future maladaptation. The strength of this study is the interdisciplinary and comprehensive set of methods across multiple fields to understand the genomics of a rafting species, T. martensii. The authors integrate ecological and population genetics, dispersal simulation, ocean current modeling, and salinity stress experiments. The most robust findings include the recovery of 5 populations of T. martensii in the Indo-Australian Archipelago, determination of divergence times and effective population size for this species, and identification of wood-feeding genes required for energy intake during long rafting journeys.\nWhile there is evidence supporting the above mentioned conclusions and claims in the manuscript, there are a number of other claims that could be problematic. In general, this study highlights the importance of plasticity for the species in variable salinity environments, but the genomic maladaptation framework assumes no plasticity in the prediction. This contradiction calls into question the validity of the genomic offset analysis. The salinity stress experiment and subsequent plasticity analyses also use a terrestrial species rather than the marine rafting species, calling into question the applicability of these results to rafting species. Finally, the dispersal simulation results are unclear. More detail on these issues is provided below.\nResponse: Thank you very much for your constructive comments, which have helped us to improve and clarify our Manuscript:\n(1) We now note more clearly that plasticity and genetic offset analyses together yield complementary perspectives, providing insights into short term versus longer term change, respectively. Specifically, the plastic gene-expression shifts detected by our experimental studies represent immediate responses to rapid environmental variation, shifts which are essential for individuals undertaking transoceanic dispersal. By contrast, the strong signals of selection detected in our genomic data, the focus of offset analysis, highlight the importance of longer-term evolutionary responses to broader environmental variation.\nWhile offset analyses inherently focus on genomic change, rather than plasticity, the strong signatures of selection obtained from genomic data (e.g. outlier osmoregulatory loci linked to regional salinity gradients) emphasise the important role of adaptive change in this rafting species. We thus feel both short-term (plastic) and longer-term (evolutionary) change together provide crucial insights into how this species copes with rafting in a dynamic ocean.\nThis integration is now more clearly highlighted in the Discussion section, where we propose that plastic responses to short-term environmental shifts on the rafting journeys are complemented by longer-term adaptive processes, providing a crucial, broader perspective on how rafting species disperse and adapt changing conditions over a range of temporal and spatial scales.\n(2) We reran all the genetic offset analyses, with multivariate LFMM, univariate LFMM, and RDA methods based on robust environmental related SNPs. We obtained similar results of genetic offsets and found that the species have a northward dispersal trend.\n(3) For the revised salinity stress experiments, we exposed the coastal species Talorchestia cf. martensii to the same salinity treatments, and sequenced their RNA to quantify expression changes. Additionally, we have added two treatments in 20 ppt salinity for 15 days and 30 days to further test for salinity plasticity relevant to long-lived natural rafting voyages. The gene expression data have all been reanalyzed, yielding numerous genes and pathways in gills exhibiting plastic expression, with higher expression shifts corresponding to larger salinity shifts and longer rafting duration (Fig. 4f–i), providing a key parallel to the results previously obtained from the terrestrial species.\n(4) We described the dispersal simulation and ΔD, and provided more details in “Results” section. “The direction and velocity of ocean currents play important roles in explaining gene flow between populations26. We defined direction deviation (ΔD) as the angle between the direction of ocean currents and the least-cost path line. The ΔD was calculated as the average of each grid cell along the least-cost path line, and the larger value represents the mismatch between directions of ocean currents and potential dispersal direction. The dispersal velocity (V ) was the average value of ds grid’s effective dispersal velocity, which was calculated as the projection of the ocean current velocity (V ) onto the direction of the path. The transit time between sites was estimated by dividing the least-cost distance by dispersal velocity (Fig. 2a).”\n(5) To further infer detailed population migration and gene flow patterns, we performed genetic simulations twice with 25 simulation models in five populations (K=5), or 20 simulation models in three populations (K=3). Both simulations suggest asymmetric gene flow, with a northward dispersal model favoured over random and southward dispersal models (Fig. 1g–i).\nLack of framework: The claim in the abstract that this paper provides a framework for studying adaptive shifts and predictability of adaptation seems a little strong; while analyses are comprehensive, no actual framework is proposed within the text.\nResponse: Our highly integrative case-study provides a comprehensive investigation into the evolution, physiology, and ecology of marine rafting in a fast-changing ocean. To this end, we now clarify the Introduction: “Thus, there is an urgent need to address how rafting species disperse, adapt, and respond to shifting climates. Here we address these questions by integrating oceanographic, physiological, and ecological genomic analyses.”\nCrucially, to provide a broad and accessible framework, we have also produced a new schematic summary of the study: Figure 6.\nFor rafting dispersal, we first identified geographical population and used gene flow analysis and ocean currents modelling to prove that the summer ocean currents drive the species' northward gene flow (Fig. 6a).\nFor their adaptation and response to shifting environments, we explore the genetic changes of long-term adaptations to future climate change. We used genetic offset, genetic-environmental associations, and population genetic differentiation to analyze their adaptive strategies that come from local adaptation. Our results first support that northward rafting decreases their maladaptation, but climate warming would weaken their habitat tracking abilities. Then we found the species would have some osmoregulatory genetic changes to the future salinity changes (Fig. 6b).\nRegarding how the rafting species adapt to the dynamic shifts on their dispersal journey, which come from the short-term fluctuated sea water salinity and nutrient poor wood habitats. According to the salinity experiments, we found the species show a higher osmoregulatory plasticity to the fluctuated sea water salinity (Fig. 6c). To obtain enough energy from the driftwood, the wood-feeding genes duplication and expression are important for their successful living.\nFig. 6 Framework of how rafting species disperse, adapt, and respond to shifting environments. a, Schematic representation of rafting-mediated dispersal between populations. b, Facing future oceanographic changes (increasing temperature and decreasing salinity), we found northward dispersal and osmoregulatory genetic changes that may facilitate their adaptation to future climate change. c, We detected osmoregulatory gene expression plasticity that helps talitrids cope with short-term salinity fluctuation during rafting events, and key gene duplications of wood-degrading enzymes that facilitate energy acquisition during rafting.\nConcerns with GEA analysis: The genotype-environment associations (GEAs) could also benefit from using the more up-to-date multivariate implementation of LFMM from the LEA3 package, and the authors should also note that recent work has shown RDA is not the most reliable method for GEA analysis, though much of the field still uses these methods (see Lotterhos 2023). GEA additionally assumes that each population occupies only a single environment, which is not true for highly migratory rafting species experiencing large environmental fluctuations.\nResponse: Thank you for your wonderful comments, and we have reanalyzed our data and added more discussion about the GEA.\n(1) Multivariate lfmm2 supports arbitrary multivariate environmental matrices (e.g., climate variables, phenotypic traits, and their interactions), enabling joint testing of complex ecological gradients. For example, in genotype–environment association studies, users can include multiple climatic factors (e.g., temperature, precipitation) and demographic covariates (e.g., age, population structure) in the design matrix X, directly estimating locus-specific effects for each variable.\nWe have performed the GEA analysis using multivariate lfmm2, univariate lfmm1 and rda methods, and then we combined all the result to select the most conservative results. We finally obtained 575 SNPs distributed in 459 genes (Supplementary 17a). These genes were enriched in several pathway in the Supplementary Table 15.\nSupplementary Fig. 17. Gene numbers and function related to salinity adaptation. a, Venn diagram depicting overlaps between outlier single-nucleotide polymorphisms (SNPs) detected using multivariate LFMM, univariate LFMM, and RDA analysis. b, Salinity exposure comparison between different populations.\nSupplementary Table 15. The enriched pathway of 459 genes from GEA analysis using multivariate LFMM, univariate LFMM, and RDA analysis.\nInput Background Corrected\n#Term ID\nnumber number P-Value\nNeuronalSystem R-DME-112316 32 176 0.000214\nGeneexpression(Transcription) R-DME-74160 62 481 0.000252\nTransmissionacrossChemicalSynapses R-DME-112315 27 139 0.000366\nMetabolism R-DME-1430728 103 973 0.000476\nSphingolipidmetabolism R-DME-428157 13 48 0.003636\nRegulationofPTENgenetranscription R-DME-8943724 11 30 0.001655\nSignalTransduction R-DME-162582 79 725 0.001733\nHSF1-dependenttransactivation R-DME-3371571 9 19 0.001733\nMetabolismofcarbohydrates R-DME-71387 22 117 0.002316\nGenericTranscriptionPathway R-DME-212436 44 343 0.003629\nImmuneSystem R-DME-168256 79 749 0.003875\nRNAPolymeraseIITranscription R-DME-73857 52 437 0.004444\nSphingolipiddenovobiosynthesis R-DME-1660661 9 23 0.004444\nCytosolicsulfonationofsmallmolecules R-DME-156584 6 8 0.004621\nActivationofNa-permeablekainatereceptors R-DME-451307 7 13 0.005362\nN-glycanantennaeelongationinthemedial/trans-Golgi R-DME-975576 7 13 0.005362\nKeratansulfatebiosynthesis R-DME-2022854 5 5 0.006893\nSLC-mediatedtransmembranetransport R-DME-425407 25 159 0.006913\nLipophagy R-DME-9613354 7 14 0.006962\nKeratansulfate/keratinmetabolism R-DME-1638074 5 6 0.009966\nDigestionofdietarylipid R-DME-192456 5 6 0.009966\nChaperoneMediatedAutophagy R-DME-9613829 7 16 0.010967\nGABAsynthesis,release,reuptakeanddegradation R-DME-888590 8 22 0.011253\nDevelopmentalBiology R-DME-1266738 19 110 0.01134\nSPhase R-DME-69242 15 76 0.013336\nNeurotransmitterreceptorsandpostsynapticsignal\nR-DME-112314 15 77 0.014384\ntransmission\nTranscriptionalRegulationbyTP53 R-DME-3700989 21 133 0.015203\nAxonguidance R-DME-422475 14 70 0.016501\nMetabolismofamine-derivedhormones R-DME-209776 7 18 0.016501\nMultifunctionalanionexchangers R-DME-427601 6 13 0.018165\nPresynapticfunctionofKainatereceptors R-DME-500657 7 19 0.018576\nNeurotransmitterreleasecycle R-DME-112310 12 56 0.019751\nTransporttotheGolgiandsubsequentmodification R-DME-948021 18 110 0.020379\nThyroxinebiosynthesis R-DME-209968 6 14 0.021644\nMismatchRepair R-DME-5358508 3 9 0.180575\nPTENRegulation R-DME-6807070 11 50 0.024139\nMembraneTrafficking R-DME-199991 38 328 0.025242\nCellularresponsetoheatstress R-DME-3371556 12 59 0.025242\nmRNASplicing-MajorPathway R-DME-72163 21 144 0.025242\nPhaseII-Conjugationofcompounds R-DME-156580 8 28 0.025242\nAUF1(hnRNPD0)bindsanddestabilizesmRNA R-DME-450408 8 28 0.025242\nVesicle-mediatedtransport R-DME-5653656 38 330 0.025763\nmRNASplicing R-DME-72172 21 146 0.027334\nAntigenprocessing:Ubiquitination&Proteasome\nR-DME-983168 21 146 0.027334\ndegradation\nPIP3activatesAKTsignaling R-DME-1257604 14 79 0.029208\nMetabolismoflipids R-DME-556833 41 369 0.029208\nRegulationofmRNAstabilitybyproteinsthatbind\nR-DME-450531 11 53 0.029928\nAU-richelements\nClassIMHCmediatedantigenprocessing&presentation R-DME-983169 21 149 0.031942\nAdaptiveImmuneSystem R-DME-1280218 27 213 0.032105\nGlycosaminoglycanmetabolism R-DME-1630316 9 38 0.0334\nRoleofphospholipidsinphagocytosis R-DME-2029485 4 6 0.03548\nPost-translationalproteinmodification R-DME-597592 67 695 0.036341\nSynthesisofDNA R-DME-69239 12 65 0.039851\nAutophagy R-DME-9612973 13 75 0.042572\nMicroautophagy R-DME-9615710 7 25 0.042572\nActivationofCa-permeableKainateReceptor R-DME-451308 7 25 0.042572\nIonotropicactivityofkainatereceptors R-DME-451306 7 25 0.042572\nO-linkedglycosylationofmucins R-DME-913709 8 33 0.045582\nReactionsspecifictothecomplexN-glycansynthesis\nR-DME-975578 4 7 0.046684\npathway\nRegulationofHSF1-mediatedheatshockresponse R-DME-3371453 10 50 0.047671\nVitaminB5(pantothenate)metabolism R-DME-199220 6 19 0.047671\nRNApolymeraseIItranscribessnRNAgenes R-DME-6807505 11 59 0.048624\nO-linkedglycosylation R-DME-5173105 9 42 0.049576\n(2) As far as GEA assumes each population occupies only a single environment, which is not true for highly migratory rafting species. According to the population genomic analysis, we found clear geographical population structure within the species, which points to potential local adaptation to the differential environments. We also found many large-scale dispersal organisms showing local adaptation to their breeding ranges, like migratory bird (Bay et al., 2018, Science; Obiol et al., 2025), migratory Arctic fish (Layton et al., 2021).\nWe have added Discussion as “Despite the importance of rafting processes in marine ecosystems, few previous studies have combined rafting and genomic forecasting in the evaluation of future dispersal potential48. Genetic offset analyses can provide a powerful means to explore the alignment between locally adapted genotypes and the genetic requirements of future marine environments, as highlighted by recent analysis of marine invertebrates49, migratory birds50,51, and fishes52 to evaluate the climate maladaptation. While such offset analyses are yet to fully account for organismal plasticity, the substantial population genomic structure detected here in T. martensii highlights the importance of local adaptation in this rafting species.”\nLimitations of Genomic offsets analysis: Genomic offset modeling assumes no plasticity, so the fact that significant plasticity was found under experimental conditions undermines the genomic offset results and should be discussed. Furthermore, much like GEA, genomic offset assumes that populations occupy and are locally adapted to a single environment. Again, as a rafting species, T. martensii experience large environmental fluctuations and local adaptation may be limited in this species. Finally, while forward and reverse genetic offset methods are appropriate for species that can disperse, how dispersal and genomic offset modeling are integrated using ∆D is unclear from the methods. In general, genomic offset predictions have only been empirically tested in a few systems (mostly trees), and these tests have yielded results spanning from inaccurate predictions to relatively accurate predictions. Some empirical validation is needed before using genomic offsets in this system to predict maladaptation.\nResponse: (1) While gene-expression plasticity is clearly important for rafting individuals experiencing rapid salinity fluctuations, our genomic analysis also reveals strong outlier loci associated with regional salinity variation, revealing an unequivocal signature of local adaptation. In this context, genetic offsets quantify the alignment between this locally adapted genetic composition and the genetic requirements of future environments, which represents the genetic variations required for long-term adaptation of future climate change.\nWhile genomic offset analyses inherently focus on long-term genetic adaptation and do not explicitly incorporate plasticity, our study proposed plasticity as a complementary mechanism acting on distinct temporal scales (short-term dispersal vs. long-term settlement). This integration strengthens our understanding of their adaptive strategies crucial across different spatial and temporal scales. Rather than undermining the validity of genomic offset analyses (which have been undertaken for many marine species), we feel the consideration of both plasticity and genomic offset is valuable for forecasting responses to future environmental scenarios.\n(2) As far as the dispersal and genomic offset modelling, we used forward genetic offset to estimate the target location of one species disperse to. Then, we calculated the ∆D between the source and target location to estimate their dispersal possibility, and test whether ocean currents support T. martensii rafting to track suitable habitats. The direction deviation (ΔD) was defined as the angle between the direction of ocean currents and the least-cost path line. The ΔD was calculated as the average of each grid cell along the least-cost path line, and the larger value represents the mismatch between directions of ocean currents and potential dispersal direction. The dispersal velocity (V ) was the average value of ds grid’s effective dispersal velocity, which was calculated as the projection of the ocean current velocity (V ) onto the direction of the path. The transit time between sites was estimated by dividing the least-cost distance by dispersal velocity (Fig. 2a).\nWe have changed the sentence in the Materials and methods section as “As the GF algorithm and forward genetic offset methods can estimate the target location of species dispersal, we integrated the quantified ocean dynamics with forward offset prediction to test whether the ocean currents can support such dispersal.”\n(3) While it is conceivable that such genomic offset analyses could prove inaccurate, our study predicts northward dispersal to track suitable habitats. We note that a variety of marine species, from phytoplankton to fish and marine mammals, have shown direct evidence of such poleward dispersal under climate warming (Hastings et al., 2020). Indeed, several recent studies have combined genetic offset and genomic forecasting in the evaluation of future marine dispersal potential in the marine organisms, for example, ascidians invasion risks along Chinese coasts under climate change (Chen et al., 2024).\nWe added some Discussion in the revised manuscripts: “Under global change, rafting is likely to present a key mechanism for many marine species tracking their shifting habitats. Along these lines, northward dispersal of T. martensii could potentially mitigate the projected genomic maladaptation of this species under future oceanographic scenarios (Fig. 3b, Supplementary Fig. 13e-g). A wide range of marine species, from phytoplankton to fish and marine mammals, have shown direct evidence of such poleward dispersal under climate warming56.”\nWe also added some Discussion “Broadly, while our integrative analyses suggest that T. martensii has potential to disperse and adapt to the thermal and osmoregulatory challenges provided by global change (Fig. 6b), shifting oceanographic connectivity patterns may limit its ability to do so. Given the vast numbers of rafting marine species likely to be impacted by forecast marine climate change, there is also an urgent need for further data to help assess and forecast such biological responses. We particularly emphasize the need for additional common garden experiments and temporal genomic studies to further inform marine biodiversity change forecasts and potential management actions.”\nInterpretability of gene flow and dispersal modeling: The dispersal and current modeling analyses are unclear, specifically what ∆D actually means - while the authors explain that this is a dispersal direction deviation, it is unclear why dispersal direction would be different from current direction if this species truly lacks autonomous dispersal and instead disperses passively. Figure 2 shows the quantified model and particle trajectory simulations, but does not clarify how ∆D is calculated - is it averaged across all modeled cells, or specific to just the final cell? Furthermore, when doing the gene flow model testing, the authors use K=3 despite recovering K=5 in previous analyses.\nResponse: (1) The dispersal and currents modeling methods were inspired by the Hodel et al., 2018. It used the Euclidean distance of shortest path between two sites, which should be straight line theoretically. However, islands and lands exist between two sites, and the coastal species couldn’t disperse on the land, therefore, we used the least-cost path on the ocean to represent the dispersal way. The ocean currents variables on the shortest path could be measured its importance in explaining geographic patterns of genetic variation (Hodel et al., 2018). From the source site to the sink site, we identified every grid cell that intersected the least-cost path line. If the ocean currents bearing of a pixel were identical to the bearing of the line, the pixel would be given a score of zero. If the bearing of the pixel were different from the bearing of the line, it would be given a score between zero and 180°, where 180° implies that the ocean current bearing of the pixel is opposite to the line connecting the populations. The ΔD was calculated as the average of each grid cell along the least-cost path line, and the larger value represents the mismatch between directions of ocean currents and potential dispersal direction. Under such a designation, a ∆D value larger than 90° precludes dispersal along the least-cost path.\nWe have added more definition about the ∆D as: “To test relationships between gene flow and ocean currents, we developed a novel model to quantify ocean currents between two sites based on least-cost dispersal path (Fig. 2a). The direction and velocity of ocean currents play important roles in explaining gene flow between populations26. We defined direction deviation (ΔD) as the angle between the direction of ocean currents and the least-cost path line. The ΔD was calculated as the average of each grid cell along the least-cost path line, and the larger value represents the mismatch between directions of ocean currents and potential dispersal direction. The dispersal velocity (V ) was the average value of ds grid’s effective dispersal velocity, which was calculated as the projection of the ocean current velocity (V ) onto the direction of the path. The transit time between sites was estimated by dividing the least-cost distance by dispersal velocity (Fig. 2a).”\n(2) As far as gene flow model testing, we defined groups based on population structure (K=3) to simplify the model prior in the fastsimcoal analysis. Three groups are easy to design genetic models to test the directional or random gene flow. In the revised manuscript, we also tried to use five group (K=5) to perform the simulations, We listed 25 models and also found the northward model is the most likely dispersal models.\nWe have added all these methods in the Materials and methods section as “We first used five populations (K=5) to infer the potential dispersal hypothesis, and listed 25 models for simulation. The model fitness was evaluated by simulating SFS for each demographic model and comparing the simulated and observed likelihoods. The best demographic model was determined based on the maximum likelihood value. During this process, we found the dispersal hypothesis and model prior could be very complex.”\nWe also added results in the Results section and supplementary files. “To further infer detailed population migration and gene flow patterns, we performed genetic simulations twice using Fastsimcoal 2.725 with 25 simulation models in five populations (K=5), or 20 simulation models in three populations (K=3) (Supplementary Figs. 5 and 6, Supplementary Table 5). Both simulations suggest asymmetric gene flow, with a northward dispersal model favoured over random and southward dispersal models (Fig. 1g–i).”\nSupplementary Fig. 6. The tested models in fastsimcoal2 about dispersal and gene flow between five groups (K = 5). SS, EB, SB, NB, VT mean South Sulawesi, East Borneo, South Borneo, North Borneo, Vietnam-Thailand; arrows represent the dispersal direction. Numbers under each model were the MaxEstLhood. Light orange background shows the northward dispersal, light blue shows the southward dispersal, and light green shows the random dispersal. The red pentagram shows the best model with the maximum likelihood value.\nLimitations of Plasticity analysis: The salinity stress experiments use a terrestrial species from a different genus than T. martensii, which is the focal species of the rest of the manuscript. While the plasticity analyses themselves seem sound, this calls into question the applicability of these results to marine rafting species, particularly T. martensii. At a minimum the discussion should address how applicable the salinity stress experiments are to the focal species of the manuscript. Additionally, there is significant evidence that plasticity can promote acclimation, but theoretical work suggests that it can also be harmful in the long term by precluding the need for species to genetically adapt to stress (see Ghalambor et. al. 2007). The authors should caveat these results in relevant sections of the discussion.\nResponse: In the revised manuscript, we collected Talorchestia cf. martensii for salinity experiments and sequenced their transcriptomes, and used the same salinity and duration designation. Moreover, we used longer time treatments to further test for salinity plasticity relevant to long-lived natural rafting voyages. \n(1) In the coastal Talorchestia cf. martensii, we identified numerous differentially expressed genes and pathways among different salinity treatments, which are related to osmoregulation and are differentially expressed only in the gills, not in the legs (Supplementary Table 18), e.g., amino acid transmembrane transport (five genes, GO:0003333), potassium ion transport (four genes, GO:0006813), cell volume homeostasis (three genes, GO:0006884), chloride transport (three genes, GO:0006821). By comparing top10 GO pathways, we found the gills of Talorchestia martensii and Morinoia aosen show a similar function in our analysis.\nFig. 4. Genetic changes in T. martensii associated with osmoregulation promote rafting acclimations to salinity shifts. f, The top 10 GO enrichment pathways for the upregulated expression genes in gills. Values in the circle represent the gene number.\nTherefore, we added Results description as “According to the differentially expressed genes (DEGs) and gene enrichment analysis between 5 ppt and 35 ppt treatments, the top 10 GO pathways in gills are primarily related to osmoregulation (Fig. 4f, Supplementary Tables 17 and 18). Notably, pathways such as amino acid transmembrane transport (five genes, GO:0003333), potassium ion transport (four genes, GO:0006813), and cell volume homeostasis (three genes, GO:0006884) were differentially expressed exclusively in gills rather than legs. Thus, we propose that amphipod gills may be crucial osmoregulatory organs, with comparable differential expression pattern between tissues also detected in the terrestrial talitrid, M. aosen (Fig. 4f, Supplementary Fig. 19b).”\nWe also added Discussion as “We found increasing differential expression of genes in response to decreasing salinity, with evidence that gills function as the major osmoregulatory organs in maintaining homeostasis under fluctuating salinity58,59.”\nSupplementary Table 18. The significant enriched pathways based on the identified differentially expressed genes only in gills but not in legs.\nInput Gene Corrected\n#Term ID\nnumber Ratio P-Value\nstructuralconstituentofchitin-basedcuticle GO:0005214 10 0.3125 0.0000\ntransferaseactivity,transferringhexosylgroups GO:0016758 10 0.2500 0.0000\ncellularresponsetounfoldedprotein GO:0034620 7 0.5385 0.0000\nUDP-glycosyltransferaseactivity GO:0008194 9 0.2813 0.0000\nlysosome GO:0005764 11 0.1930 0.0000\nproteolysis GO:0006508 26 0.0714 0.0000\nextracellularspace GO:0005615 32 0.0588 0.0000\nintracellularmembrane-boundedorganelle GO:0043231 13 0.1238 0.0001\nsulfotransferaseactivity GO:0008146 5 0.6250 0.0001\nsmoothseptatejunction GO:0005920 4 0.5714 0.0014\nvacuole GO:0005773 4 0.5714 0.0014\napicolateralplasmamembrane GO:0016327 5 0.2941 0.0019\ncarbohydratebinding GO:0030246 8 0.1356 0.0020\ncell-matrixadhesion GO:0007160 4 0.4444 0.0026\nbodymorphogenesis GO:0010171 5 0.2500 0.0033\ncarbohydratehomeostasis GO:0033500 4 0.4000 0.0033\nresponsetoreactiveoxygenspecies GO:0000302 4 0.4000 0.0033\nnursecellapoptoticprocess GO:0045476 4 0.3636 0.0043\nsphingolipidmetabolicprocess GO:0006665 4 0.3636 0.0043\nresponsetounfoldedprotein GO:0006986 4 0.3333 0.0055\naminoacidtransmembranetransporteractivity GO:0015171 5 0.2083 0.0062\ntissuedevelopment GO:0009888 4 0.3077 0.0069\nserine-typecarboxypeptidaseactivity GO:0004185 3 0.6000 0.0081\ncysteine-typeendopeptidaseactivityinvolvedin\nGO:0097200 3 0.6000 0.0081\nexecutionphaseofapoptosis\nproteinfoldingchaperone GO:0044183 4 0.2857 0.0081\naminoacidtransmembranetransport GO:0003333 5 0.1786 0.0101\ncysteine-typeendopeptidaseactivityinvolvedin\nGO:0097153 3 0.5000 0.0105\napoptoticprocess\nsteroltransport GO:0015918 3 0.5000 0.0105\ncellvolumehomeostasis GO:0006884 3 0.5000 0.0105\nprogrammedcelldeathinvolvedincelldevelopment GO:0010623 3 0.5000 0.0105\nmetallocarboxypeptidaseactivity GO:0004181 4 0.2500 0.0109\nphosphatidylinositolbisphosphatebinding GO:1902936 4 0.2353 0.0129\npotassiumiontransport GO:0006813 4 0.2353 0.0129\ndopaminemetabolicprocess GO:0042417 3 0.4286 0.0140\ndefenseresponsetoGram-negativebacterium GO:0050829 8 0.0899 0.0159\ncationbinding GO:0043169 3 0.3750 0.0180\nimmuneresponse GO:0006955 5 0.1429 0.0203\noxidation-reductionprocess GO:0055114 15 0.0547 0.0216\ngermarium-derivedfemalegerm-linecystformation GO:0030727 3 0.3333 0.0223\nlipidtransporteractivity GO:0005319 3 0.3333 0.0223\nregulationofimaginaldisc-derivedwingsize GO:0044719 3 0.3333 0.0223\nmemory GO:0007613 4 0.1818 0.0247\napoptoticprocess GO:0006915 5 0.1316 0.0252\nbasalplasmamembrane GO:0009925 3 0.3000 0.0261\nnervoussystemdevelopment GO:0007399 7 0.0897 0.0265\nheatshockproteinbinding GO:0031072 4 0.1739 0.0267\ncysteine-typeendopeptidaseactivity GO:0004197 4 0.1667 0.0299\nprogrammedcelldeath GO:0012501 3 0.2727 0.0308\ndefenseresponsetobacterium GO:0042742 5 0.1220 0.0310\ngenesilencingbymiRNA GO:0035195 3 0.2500 0.0355\nchloridetransport GO:0006821 3 0.2500 0.0355\nalkalinephosphataseactivity GO:0004035 3 0.2500 0.0355\nscavengerreceptoractivity GO:0005044 3 0.2500 0.0355\nmuscleattachment GO:0016203 4 0.1538 0.0355\nchaperonecofactor-dependentproteinrefolding GO:0051085 4 0.1538 0.0355\ncellsurface GO:0009986 5 0.1136 0.0379\nlateralplasmamembrane GO:0016328 3 0.2143 0.0469\npotassiumionimportacrossplasmamembrane GO:1990573 3 0.2143 0.0469\n(2) To show the transcriptional plasticity, we also performed DAPC analysis to compare the magnitude of expression shift, which represents a quantitative measure of gene expression plasticity. We found salinity shift from 35 ppt to 10 ppt would lead to larger expression shifts than salinity changes from 20 ppt to 5 ppt (P < 0.001) MCMC (Fig. 4k). For the same salinity shifts (20 ppt to 5 ppt), larger gene expression shifts can be found in longer duration treatments (72 h) than shorter duration treatments (24 h; P < 0.001). We have added Results as “These dynamic expression changes MCMC highlight transcriptional plasticity36 (Fig. 4g–l, Supplementary Fig. 19c–f), with the stronger shifts linked to the larger salinity shifts (Fig. 4k) and longer treatments (Fig. 4l).”\nFig. 4 Genetic changes in T. martensii associated with osmoregulation promote rafting acclimations to salinity shifts. k-l, Discriminant analysis of principal components (DAPC) indicate gene expression plasticity in response to salinity shifts in same durations (k) and duration changes in same salinity shifts (l).\nAccording to comparing these results generated from salinity experiments of Talorchestia martensii with the former analysis in the terrestrial species Morinoia aosen, we obtained similar plasticity patterns (Supplementary Fig. 19c–f).\nSupplementary Fig. 19 The expression plasticity found in Morinoia aosen. a, A clear divergence in expression patterns between gills and legs by the cluster analysis. b, The top 10 GO enrichment pathways for the DUGs specifically in gills. Values in the circle represent the gene number enriched in each pathway. c, The number of DEGs increased along salinity gradient treatments between different treatments. d, Ternary plots depicting the proportion of each gene’s expression variance between 5, 20, and 35 ppt salinity treatments, and more genes show higher expression at 35 ppt. e, f, Discriminant analysis of principal components (DAPC) for variation in gene expression plasticity in response to salinity shifts (e) and duration changes (f).\n(3) To further test for salinity plasticity relevant to rafting voyages, we added long time treatments, such as 15 and 30 days in 20 ppt. Relatively strong expression shifts were detected at low salinities (Fig. 4h) and over the longest treatments (Fig. 4i). Moreover, the number of genes enriched in ion channel transport pathways increased as treatment time increased (Fig. 4j, Supplementary Table 19). Therefore, we thought the Talorchestia martensii also used more genes involved in ion transport to keep the ion balance on the rafting journey as the duration time increased.\nFig. 4 Genetic changes in T. martensii associated with osmoregulation promote rafting acclimations to salinity shifts. h, Ternary plots depicting the higher proportion of gene expression at 5 ppt than 10 and 20 ppt. i, Ternary plots depicting the proportion of gene expression increased along the duration from 72 h to 15 d and 30 d compared with 24 h and 48 h. j, The number of genes that enriched in ion channel transport pathways along the treatment duration. e.g., R-DME-983712: Ion channel transport; GO:1902600: proton transmembrane transport; details in Supplementary Table 19.\nSupplementary Table 19. The gene pathways where the number of genes increased with duration time.\n#Term Database ID 30d 15d 72h 48h 24h\nIonchanneltransport Reactome R-DME-983712 8 14 7 3 NA\nOxygen-dependentprolinehydroxylationof\nReactome R-DME-1234176 4 2 1 1 NA\nHypoxia-inducibleFactorAlpha\nDetoxificationofReactiveOxygenSpecies Reactome R-DME-3299685 4 3 3 2 1\nO2/CO2exchangeinerythrocytes Reactome R-DME-1480926 3 2 1 NA NA\nTransportofnucleosidesandfreepurineand\nReactome R-DME-83936 2 1 2 1 NA\npyrimidinebasesacrosstheplasmamembrane\nMembraneTrafficking Reactome R-DME-199991 15 16 9 11 1\nVasopressinregulatesrenalwaterhomeostasis\nReactome R-DME-432040 2 1 1 NA NA\nviaAquaporins\nResponsetoelevatedplateletcytosolicCa2+ Reactome R-DME-76005 4 5 6 4 1\nplasmamembraneproton-transportingV-type\nGO GO:0033181 6 8 6 1 NA\nATPasecomplex\nproton-transportingATPaseactivity,rotational\nGO GO:0046961 7 10 8 2 1\nmechanism\nbasementmembrane GO GO:0005604 5 8 3 NA NA\nprotontransmembranetransport GO GO:1902600 7 12 10 2 2\nvacuolarproton-transportingV-typeATPase,\nGO GO:0000221 3 3 2 1 NA\nV1domain\nmitochondrialcalciumionhomeostasis GO GO:0051560 2 1 1 1 NA\n(4) Specifically, the V-type H+ ATP genes related to proton transport showed decreased expression with increasing treatment duration (VATA, VATB, VATE, VATG, VATO). Notably, genes involved in the transport of small molecules (Tret1, ABCA3, HXT1) and inorganic ion (Sodium/hydrogen exchanger, Sodium/Potassium/Chloride Transporter and KCAB2) showed increased expression in these treatments (Fig. 4m, Supplementary Fig. 20, and Supplementary Table 20).\nFig. 4 Genetic changes in T. martensii associated with osmoregulation promote rafting acclimations to salinity shifts. m, Hypothesized model of ion transport across the gill based on genes with increased (red) and decreased (blue) expression over duration gradient treatments.\nSupplementary Table 20. The genes show gradient increased or decreased expression changes along the increased duration.\nGene 24h 48h 72h 15d 30d Function\nVATO 175.82 124.81 131.40 111.05 100.54 protontransmembranetransporter\nVATA 158.05 117.68 118.34 62.81 68.43 protontransmembranetransporter\nVATE 245.53 148.51 192.51 141.63 107.71 protontransmembranetransporter\nVATB 295.06 193.84 216.91 128.91 125.83 protontransmembranetransporter\nVATH 98.01 73.46 69.64 30.52 27.06 protontransmembranetransporter\nCLCA1 4.72 3.60 3.79 2.05 1.53 Calcium-DependentChlorideChannel\nAQP4 26.97 30.53 26.32 12.67 9.26 Aquaporin\nS15A2 32.43 23.41 5.85 3.30 2.17 H+/PeptideTransporter\nVATG 301.54 154.53 165.24 192.62 137.30 protontransmembranetransporter\nTret1 0.67 0.76 0.70 1.05 2.43 Trehalosetransporter\nSodium/hydrogen\n2.44 3.27 4.05 4.19 8.94 Sodium/hydrogenexchanger\nexchanger\nInsP3R 0.97 1.36 2.26 1.98 2.86 Ca2+channelactivatedbyinositoltrisphosphate\nRyanodine\n5.94 9.45 13.27 16.68 14.16 intracellularcalciumchannels\nreceptor\nDscam2 0.23 0.26 0.54 1.41 1.37 largesingle-passtransmembraneproteins\nKCAB2 2.61 2.37 3.43 5.65 7.99 PotassiumVoltage-GatedChannel\nABCA3 0.00 0.00 0.00 4.29 18.82 ATP-bindingcassette(ABC)transporters\nTMC7 1.88 3.26 1.97 4.90 14.05 TransmembraneChannel-LikeProtein7\nSLC6A5 2.59 4.44 4.30 4.15 13.78 sodium-andchloride-dependentglycinetransporter2\nSLC28A3 0.00 0.00 0.00 3.48 23.75 ConcentrativeNa(+)-NucleosideCotransporter\nHXT1 0.01 0.00 0.00 0.33 5.17 hexosetransport\nHXT1 0.00 0.00 0.29 8.82 103.07 hexosetransport\nHXT1 0.00 0.00 0.00 2.27 16.94 hexosetransport\nS12A2 1.08 0.19 1.39 1.54 10.73 Sodium/Potassium/ChlorideTransporter\n(5) As far as plasticity can promote acclimation or harmful in the long term, we added some discussion about plasticity in the Discussion section as: “We found increasing differential expression of genes in response to decreasing salinity, with evidence that gills function as the major osmoregulatory organs in maintaining homeostasis under fluctuating salinity58,59. Rapid gene expression shifts suggest plasticity may be a key prerequisite for successful rafting, at least in the short term60 (Fig. 6c). Nevertheless, theoretical studies suggest that, in the longer term, plasticity may slow or constrain populations’ adaptive genetic response to environmental stress61. Therefore, we propose that both short-term (plastic) responses to temporal environmental variation, and longer-lasting adaptive responses to local conditions, should be integrated to understand how rafting species may disperse and respond to future climates.”\nMissing details and code: Because the manuscript covers so many topics and fields, detail is sometimes lacking in the main text as to what some of the terms and results actually mean. For example, the distinction between forward and reverse genetic offset is unclear in the main text, and this section would benefit from definitions of these terms that are biologically relevant. This was just one example, but other result and discussion sections would also benefit from additional detail. While the text is already long, clarifying terminology and what each analysis can actually do is important given that it is unlikely that any one reader will be an expert in every field covered in the manuscript.\nResponse: In our revised manuscripts, we have added the purpose of each analysis, explained more terminology with biologically relevant languages, and streamlined the Discussion section.\nFor example, about genetic offsets, we added explanation as “To test how rafting crustaceans will respond to future oceanographic change (in the absence of plasticity), we calculated local, forward and reverse genetic offsets29. Specifically, local offset represents the mismatch between current and predicted future genomic variation in the absence of migration. By contrast, forward genetic offset represents the minimum genetic change a population would experience if migrating to a relatively favorable future location, with high values indicating strong mismatches between contemporary and future genetic compositions. Thirdly, reverse genetic offset seeks to assess whether the genotypes required to adapt to future environments already exist within contemporary populations. To this end, we extracted a total of 117,113 SNPs with minor allele frequency > 10%, to predict genetic offset using the gradientForest (GF) approach, a regression tree-based machine-learning algorithm30.”\nAs far as ΔD, we have explained as “To test relationships between gene flow and ocean currents, we developed a novel model to quantify ocean currents between two sites based on least-cost dispersal path (Fig. 2a). The direction and velocity of ocean currents play important roles in explaining gene flow between populations26. We defined direction deviation (ΔD) as the angle between the direction of ocean currents and the least-cost path line. The ΔD was calculated as the average of each grid cell along the least-cost path line, and the larger value represents the mismatch between directions of ocean currents and potential dispersal direction. The dispersal velocity (V ) was the average value of grid’s effective dispersal velocity, which was calculated as the projection of the ocean current velocity (V ) onto the direction of the path. The transit time between sites was estimated by dividing the least-cost distance by dispersal velocity (Fig. 2a).”\nThe written methods section is detailed and thorough. However, the provided code does not include all of the analyses; scripts for the DEG analysis, DAPC, selective sweep analysis, quantification of ocean current dynamics, gene ontology enrichment analysis, LFMM, and RDA were not in the zip file - additionally, the demographic analysis and ocean current drift modeling have no provided scripts, though again, the methods are detailed in the text. Inclusion of code used to conduct these analyses is necessary to ensure the reproducibility of the work and to ensure that the analysis is sound.\nResponse: We have submitted 11 R, Python or Shell scripts used in this study to figshare https://figshare.com/s/e164e7f678e789439b69. We also upload a file named readme.txt to introduce each script used in which analysis. For example, “1.Population_genetics.txt”, the script was used for population genetic analysis like SNP calling, filtering, TreeMix, Admixture, Dsuite, and Selective sweep analysis. “8.DEGs.R” is the R script used for DEG analysis.\nOverall: Overall, the work integrates concepts across the fields of oceanography and ecological genomics to gain a more comprehensive understanding of how rafting species disperse, adapt, and respond to shifting climates, particularly ocean currents and salinity. While some of the conclusions of the work are well supported by the analyses conducted, many others are not, namely the genetic offsets, the salinity experiments and subsequent plasticity analyses, and the dispersal modeling. In addition to potential reanalysis of these data and separation of the experimental salinity work from the population genomic and dispersal/current modeling work, this manuscript would benefit from a more nuanced discussion of the caveats and limitations of the methods used and the results found, and from more explicit descriptions and definitions of relevant terminology and analyses in the main body of the text.\nResponse: Thank you very much for your constructive comments, which greatly improved our manuscript. We have revised the manuscript thoroughly, including reanalyzed the genomic offset section, used the coastal species to perform salinity experiments and plasticity analysis. We also added several discussions about the plasticity and local adaptation, limitation of GEA, and added more biological descriptions to some terminology, and moved the descriptions of methods and software to the Methods section. We also streamlined the Discussion section to strengthen the framework in Figure 6.\nDetailed comments:\nLine 81: “Based on these high-quality SNPs” is a bit unclear, do the authors mean that they used the high quality SNPs to do the analyses they mention next?\nResponse: We used stringent quality filter methods to filter SNPs, including keep only biallelic SNPs with the following filtering criteria: (1) minQ > 30, (2) read depth > 5, (3) max-missing counts = 5, (4) SNPs at least 5 bp away from indel regions. we also performed a linkage disequilibrium based SNP pruning to keep high quality SNPs.\nTherefore, we changed our description as “Based on retained SNPs after stringent quality-filtering”.\nLine 129: Given the claim of stepping stone dispersal, it would be nice to do an analysis confirming IBD\nResponse: We have done Mantel tests to investigate the relationship between geographical distance and gene flow, but we didn’t find significant evidence of isolation by distance (Fig. 1g, Mantel statistic r=-0.06306, P=0.53333).\nHowever, when we used the ΔD and V that calculated from quantified ocean currents methods to plot with gene flow, we found the ΔD and V of ocean currents in June show a significant correlation with gene flow (Fig. 2c), and the mantel test also show significant relationship (ΔD: Mantel statistic r=0.6836, P=0.016667; V : Mantel statistic r=0.6617, P=0.016667). These suggest the species’ dispersal was not random but driven by the ocean currents.\nLine 140: include a more definitive definition of genomic offset to better characterize what these results mean - i.e. what would a higher or lower offset mean? How does gradient forest define genomic offset?\nResponse: (1) We have explained the terminology in the “Results” section as follows. “To test how rafting crustaceans will respond to future oceanographic change (in the absence of plasticity), we calculated local, forward and reverse genetic offsets29. Specifically, local offset represents the mismatch between current and predicted future genomic variation in the absence of migration. By contrast, forward genetic offset represents the minimum genetic change a population would experience if migrating to a relatively favorable future location, with high values indicating strong mismatches between contemporary and future genetic compositions. Thirdly, reverse genetic offset seeks to assess whether the genotypes required to adapt to future environments already exist within contemporary populations.”\n(2) In the “Materials and methods” section, we claimed how does gradient forest define genomic offset as “Based on the identified environment-associated SNPs, we used the gradient forest models to calculate the genetic offset and predict how allele frequencies would shift under future oceanographic change. We transformed the oceanographic variables from each of the four climate scenarios (2050 SSP1-2.6, 2100 SSP1-2.6, 2050 SSP5-8.5, and 2100 SSP5-8.5) into genetic importance and calculated the Euclidean distance between the present and future genomic importance values for each grid, which served as a metric of genetic offset.”\nLine 146: Clarify the difference between local, forward, and reverse offset in biologically relevant terms\nResponse: We have added more explanation in the “Results” section in our revised manuscript as: “Specifically, local offset represents the mismatch between current and predicted future genomic variation in the absence of migration. By contrast, forward genetic offset represents the minimum genetic change a population would experience if migrating to a relatively favorable future location, with high values indicating strong mismatches between contemporary and future genetic compositions. Thirdly, reverse genetic offset seeks to assess whether the genotypes required to adapt to future environments already exist within contemporary populations.”\nLine 161: Why was the cutoff for dispersal direction deviation 90º - what is the biological relevance here?\nResponse: We defined direction deviation (ΔD) as the angle between the direction of ocean currents and the least-cost path line. The ΔD was calculated as the average of each grid cell along the least-cost path line, which would be a score between zero and 180°, and the larger value represents the mismatch between directions of ocean currents and potential dispersal direction. As the following plot, ∆D equal to 90° means the two directions are perpendicular, while ∆D larger than 90° represents the ocean currents bearings opposite to the dispersal direction.\nWe have added explanations and changed the sentence as “Based on annual mean ocean current data under 2100 SSP5-8.5 scenarios, we found that 86.95% of populations will face dispersal-direction deviations (ΔD) exceeding 90°, preventing them from reaching suitable future habitats (Fig. 3d, Supplementary Fig. 14, and Supplementary Table 14).”\nLine 173: Clarify whether identified SNPs were identified as those recovered by both methods (i.e. overlapping), or the combination of all SNPs identified, including those identified by just one method - the supplemental shows that only SNPs identified by both methods were kept, but this isn’t clear in the text\nResponse: We used the overlapping results recovered by both methods. To clarify it, we have changed the sentence as: “Accordingly, we found sea surface salinity to be strongly associated with genomic variation in T. martensii, with a total of 575 significant salinity-associated SNPs, located within 459 genes, consistently identified across multivariate latent-factor mixed model (LFMM)33, univariate LFMM, and redundancy analysis(RDA)34(Fig. 4a and Supplementary Fig. 17a).”\nLine 176: Give a number for how many genes were enriched in each of the GO categories\nResponse: We have added the number of genes and changed the sentence as “Enriched genes included those involved in plasma membrane (77 genes, GO:0005886), ion transport (five genes, GO: 0006811) and SLC-mediated transmembrane transport (19 genes, R-DME-425407) (Supplementary Table 15).”\nLine 202: “and so on” is non-specific, it would be good to mention why the first two pathways were described but not others - was there some cutoff for why some pathways made the main text and others didn’t?\nResponse: In our analysis, we used KOBAS to perform the gene pathways enrichment and kept the pathways which have more than two genes enriched and P-value less than 0.05. We have updated all the gene pathways in our revised manuscript and attached them in the Supplementary Tables 15 to 18.\nLine 237: Why was Paryhale hawaiensis used as a reference?\nResponse: Because the Parhyale hawaiensis is the marine non-rafting amphipods, which is the closest marine relatives with Talorchestia martensii.\nIn the revised manuscripts, we expanded the comparative range to nine amphipods and focus on non-rafting and rafting species. We also found the rafting species show higher proportion of GH7 gene family and significantly expansion of GH7 gene family in the branch that contains only rafting species.\nLine 260-261: It would be helpful to quickly outline what some of the expectations of Pleistocene sea-level fluctuations and Wallace’s Line are, either here or in the introduction\nResponse: we added some discussion about Pleistocene sea-level fluctuations and Wallace’s Line in the “Discussion” section: “During the early Pleistocene, sea-level decline exposed the Sundaland landmass, connecting coastlines of Vietnam-Thailand and Borneo19, and isolating the VT+SB+NB clade from the EB and SS populations (Supplementary Fig. 4). Subsequent sea-level rise created new intertidal zones around Borneo, enabling coastal talitrids to expand their ranges; the NB and SB lineages likely diverged from VT around 1.78 Ma. A subsequent split between NB and SB at ~0.66 Ma aligns with sea-level fluctuations during the Naynayxungla Glaciation (0.78–0.50 Ma)40, when land bridges between the Malay Peninsula and Borneo likely drove regional isolation. In addition to historical sea-level fluctuations, Wallace’s line has long been regarded as a barrier for marine and terrestrial species in the IAA42,43, even for species with extended pelagic larval dispersal44. In our study, the relatively deep genetic split between SS+EB and VT+SB+NB implies that modern oceanographic circulation in IAA has failed to fully erase these historical biogeographic patterns. Nevertheless, our analyses provide evidence of ongoing rafting-mediated gene flow between populations across the region (Fig. 1e, f).”\nLine 292: Careful about using the term “genetic vulnerability” or “vulnerability” when talking about genetic offset, the models can describe potential maladaptation, but how that actually translates to meaningful fitness changes and population vulnerability is still an area of active research in the field\nResponse: We agreed with you and have changed the “genetic vulnerability” or “vulnerability” as “maladaptation” throughout the revised manuscript.\nLine 315: Plasticity does not always facilitate long term population survival even if it does facilitate acclimation and survival in fluctuating environments\nResponse: Thanks for your comments. We have reconstructed this part as: “Rapid gene expression shifts suggest plasticity may be a key prerequisite for successful rafting, at least in the short term60 (Fig. 6c). Nevertheless, theoretical studies suggest that, in the longer term, plasticity may slow or constrain populations’ adaptive genetic response to environmental stress61. Therefore, we propose that both short-term (plastic) responses to temporal environmental variation, and longer-lasting adaptive responses to local conditions, should be integrated to understand how rafting species may disperse and respond to future climates.”\nLine 362: It is unclear exactly what downstream analyses the mitochondrial fragments are used for.\nResponse: We used these mitochondrial fragments for phylogeny and species delimitation in Supplementrary Fig. 1. According to the delimitation, we distinguished the focal species T. martensii in the IAA.\nWe have added the description as “We first performed species delimitation followed Liu et al. (2023)65 and inferred phylogenetic relationships based on the concatenated dataset comprising the mitochondrial (COI, COII, 16S) and nuclear (18S rRNA, 28S rRNA, H3, NaK) gene regions.”\nLine 373: What was the sequencing depth?\nResponse: The data used for draft genome assembly is around 360×, and the data used for population resequencing is at least 15×.\nWe have changed the sentence to “A total of 239.88 Gb clean reads (~360×) was obtained after removing low-quality reads with poor sequencing accuracy and adapters with Trimmomatic v.0.3966.”\nLine 374: How are “low-quality read” defined?\nResponse: A low-quality read refers to a sequencing read with poor accuracy, characterized by low Phred quality scores, high levels of ambiguous bases, or contamination. We have corrected the sentence as “A total of 239.88 Gb clean reads (~360×) was obtained after removing low-quality reads with poor sequencing accuracy and adapters with Trimmomatic v.0.3966.”\nLine 379: What is the threshold of nucleotide diversity used?\nResponse: We found the regions of high heterozygosity were removed at the lateral redundans analysis.\nLine 400: What is meant by “concatenated all SNPs” - I assume this mean SNPs for all samples, but clarifying this would be helpful\nResponse: Yes, it is the SNP from all samples. We have corrected the sentence as “To characterize population structure, we first used SNP data from all samples and constructed phylogenies using RAxML v.8.2.973.”\nLine 426: Each model was run “multiple times”, but what is the exact number there?\nResponse: We have changed as “The Markov Chain was run for a total of 50,000 iterations, and sampled every 100 iterations. The MCMC diagnostics were performed in Tracer v1.6(http://beast.bio.ed.ac.uk/Tracer).”\nLine 438: The authors define groups here based on K=3, but used K=5 for all other analyses and recovered K=5 in their population structure analyses\nResponse: We used K=3 to simplify the prior in Fastsimcoal model simulations, because when we set K=5, there are too many models (> 100) to be simulated. According to the models which generated by K=3, we also set 25 models used K=5.\nBased on the model comparisons, we found the northward dispersal from South Sulawesi to Vietnam/Thailand with South Borneo as the stepping stone is also the best models, which is similar to the K=3.\nLine 501: It is unclear to me what is meant by “transforming multidimensional climatic variables into multidimensional genetic space” - I assume this is part of variable selection, but an actual description of this method would be helpful\nResponse: The GF (Ellis et al. 2012) methods are founded on the common idea of modelling compositional turnover using nonlinear functions of environmental gradients. These turnover functions provide a means to transform environmental variables to a common biological scale of compositional turnover (in this case, turnover in allele frequencies), thereby allowing conversion from multidimensional environmental space to multidimensional genetic space while selecting and weighting variables such that they best summarise genomic variation. By applying the associated turnover function to each mapped environmental variable, these functions can be used to map expected patterns of genomic variation, including under scenarios of future climate to assess potential impacts to locally adapted genetic variants in a spatially explicit manner.(Fitzpatrick and Keller, 2015).\nTherefore, we have changed the sentence as: “Then, we used turnover function in GF models to convert complex climate data into patterns of genetic variation (allele frequencies), subsequently to select and weight these variables to best summarize genomic variation. ”\nLine 501: It seems like the authors did not scale the climate variables before using them in GEA/genetic offset analysis, which would allow for a more accurate comparison of their influence on the response variables within the analysis\nResponse: We used Z-Score standardization methods X’= (X-μ) /δ to scale the climate variables, while X is original value, μ is mean of the original data, δ is standard deviation of the original data, and X’ is the scaled standardized value. In our revised manuscripts, we have reanalyzed the whole part of genetic offset and climate analysis based on the scaled variables and all the figures were replotted based on the new data.\nFig. Comparisons of genetic offsets between unscaled environmental variables (a-c) and scaled environmental variables (d-e).\nAccording to the above figure, we found the simulations without scaled variables (a-c) and that with scaled variables (d-f) were similar. We found that populations in Souh Sulawesi, south Borneo, the coast along the Gulf of Thailand, and many of the areas along the dispersal route exhibited higher local, lower forward, and higher reverse genetic offsets, except the northern Vietnam range.\nWe further performed correlation tests between these two datasets and found a higher correlation between scaled and unscaled genetic offsets (Figure below).\nFig. Correlation tests of genetic offsets estimated from unscaled environmental variables (x-axis) and scaled environmental variables (y-axis).\nLine 551: Defining the salinity of seawater surrounding each collection site as their salinity tolerance when there is migration between populations seems counterintuitive, can there really be a site-specific salinity tolerance when the authors also describe significant plasticity and tolerance to large salinity shifts during dispersal?\nResponse: To quantify the salinity exposure of each population, which exhibits a distinct geographical distribution, we used salinity values from the sampling sites. For each population, we had at least four collection sites, and we plotted the observed salinity values to characterize the environmental salinity exposure (Supplementary Fig. 17b).\nWe have changed the sentence as “We defined the salinity of seawater surrounding each collection site as their salinity exposure and chose two populations with the most dramatic salinity exposure differences to test the selection signals related to salinity adaptation.”\nReviewer #3 (Remarks on code availability):\nThe written methods section is detailed and thorough. However, the provided code does not include all of the analyses; scripts for the DEG analysis, DAPC, selective sweep analysis, quantification of ocean current dynamics, gene ontology enrichment analysis, LFMM, and RDA were not in the zip file - additionally, the demographic analysis and ocean current drift modeling have no provided scripts, though again, the methods are detailed in the text. Inclusion of code used to conduct these analyses is necessary to ensure the reproducibility of the work and to ensure that the analysis is sound.\nResponse: We have submitted 11 R, Python or Shell scripts used in this study to figshare https://figshare.com/s/e164e7f678e789439b69. We also upload a file named readme.txt to introduce each script used in which analysis. For example, “1.Population_genetics.txt”, the script was used for population genetic analysis like SNP calling, filtering, TreeMix, Admixture, Dsuite, and Selective sweep analysis. “8.DEGs.R” is the R script used for DEG analysis.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: Thank you very much for your suggestions. Please check the point-by-point response and the revised manuscript.\nReviewer #4 (Remarks on code availability):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: We have submitted 11 R, Python or Shell scripts used in this study to figshare https://figshare.com/s/e164e7f678e789439b69. We also upload a file named readme.txt to introduce each script used in which analysis. For example, “1.Population_genetics.txt”, the script was used for population genetic analysis like SNP calling, filtering, TreeMix, Admixture, Dsuite, and Selective sweep analysis. “8.DEGs.R” is the R script used for DEG analysis.\nReviewer #5 (Remarks to the Author):\nThe study investigates how the rafting crustacean Talorchestia martensii adapts to climate change in the Indo-Australian Archipelago. The authors combined ecological modeling and genomics to explore the dispersal patterns and adaptive mechanisms of rafting species under rapid environmental changes. From 203 samples collected across the region, five geographically distinct T. martensii populations were identified. Simulations revealed stronger northward migration patterns to shape the populations compared to random or southward dispersal. Using a quantified model and particle trajectory simulations, the authors further confirmed asymmetric northward gene flow, particularly during summer monsoons, consistent with seasonal ocean currents.\nGenomic data from the identified populations revealed that salinity fluctuations acted as a major selective pressure during rafting dispersal. The authors utilized transcriptome data from salinity stress experiments to identify genes associated with salinity adaptation, such as GRIK2 and SLC6A13. However, the lack of statistical significance tests and functional validation raises concerns regarding the robustness of these findings.\nResponse: Thank you very much for your valuable and insightful suggestions. We have made major revisions as below according to your suggestions:\n(1) We have added statistical significance tests to all the analyses needed and provided the P-values or adjusted P-values. For example, in the Supplementary Fig. 20, we have added the P-values for each gene that showed differential expression with increasing treatment duration.\n(2) We constructed GRIK2 dataset based on 12 genomes of amphipods, which occur in freshwater, marine and terrestrial habitats. Phylogenetic analysis revealed that these GRIK2 sequences clustered strongly by habitat type (freshwater; marine; terrestrial) rather than reflecting the underlying species phylogeny (Supplementary Fig. 18b), providing further support for this gene’s role in amphipod salinity adaptation.\nSupplementary Fig. 18. Genetic changes of salinity-related GRIK2 genes. a, genotype–salinity associations between the allele frequency of SNPs in GRIK2 and annual mean sea surface salinity (AMS). Colors on the map are based on variations of the relevant salinity across the distribution range. ALT: Alternate allele, REF: Reference allele. b, The GRIK2 phylogeny based on 1:1 orthologs in 12 amphipod genomes suggests that species were clustered by freshwater, marine, and terrestrial habitats.\n(3) In addition to the salinity stress experiments (35, 20, 10, 5 ppt) for coastal T. cf. martensii, we performed long-time salinity treatments for T. cf. martensii to 15 days and 30 days in 20 ppt to further test for salinity plasticity relevant to long-lived natural rafting voyages. The related results can be found in Fig. 4 and the following point-by-point response.\n(4) We added more comparative genomics analyses between rafting and non-rafting species to gain a broader set of key genes involved in driftwood digestion, including searching the whole of CAZyme repertoire and comparing the CAZyme family composition between six rafting and three non-rafting species, gene family expansion and extraction analysis (Fig. 5b,c).\nFig. 5 Horizontal gene transfer and duplication events of GH gene families. a, Three major steps of wood-degrading from cellulose to glucose and associated glycoside hydrolase (GH) gene families. b, The reconstructed phylogeny based on three non-rafting and six rafting amphipod genomes. c, Comparisons of proportion of GH gene families in all cazymes indicate that rafting talitrids have more GHs than non-rafting species (Wilcoxon signed-rank test). d, Two detected duplication modes of the GH gene family in talitrids. e, The inferred phylogenetic relationships of the GH9, GH30, and GH7 proteins from talitrids, non-talitrid crustaceans, insects, and fungi. Circles on the clade represent different duplication models present in (d). Black dots on the nodes correspond to different talitrid clades. The middle left aligned sequences show the amnio acid variation of GH7 clade II (C-II) from Parhyale hawaiensis (Phaw), Hyalella azteca (Hazt), Cochinorchestia sp. (Cochi), Talorchestia martensii (Tmar), and Morinoia aosen (Maos). The middle right schemes show a summary of the GH7 evolution, including horizontal gene transfer (HGT) from fungi to crustacean ancestors and gene duplication in talitrids. f, the expression patterns of GH7, GH9 and GH30 between three species and five tissues. Each line represents one gene copy, and circle size indicates expression level. The black brace and numbers correspond to the clade in Fig. 5e.\n(5) We collected five tissues (antennae, brain, gut, leg, and gill) for RNA-sequencing and compared the expression of wood-feeding genes. We found that the GH genes showed a higher expression in the gut tissues. We also quantified the expression of GH genes between different duplication events and found that all the duplications show effective expression in the gut, with comparable expression levels in rafting species despite their greater gene copy numbers relative to non-rafting species (Fig. 5f).\nOverall, this is a well-executed study on the adaptation of rafting crustaceans, and the manuscript should be considered for publication after addressing the following comments:\nThe manuscript applies various bioinformatics tools, which are essential for obtaining reliable results. However, listing all software without explanation in the main text is distracting. I recommend moving the software to the Methods section and briefly explaining the conceptual rationale behind each analysis in the main text.\nResponse: Thanks for your suggestions. We have moved most software to the Method section and have added the purpose of each analysis, explained more terminology with biologically relevant languages, and streamlined the Discussion section.\nLine 180: Could the authors provide a list of GO terms enriched in salinity-associated genes? It is unclear whether these genes represent false positive associations or are truly enriched in certain biological functions.\nResponse: We used KOBAS to perform the gene pathways enrichment and kept the pathways that have more than two genes enriched and a P-value less than 0.05.\nWe provided the salinity-associated genes enriched pathways in the supplementary tables and added descriptions in the Results section as “Many of these genes are similarly implicated in osmoregulatory pathways, e.g., plasma membrane (10 genes, GO:0005886); transmembrane transport (four genes, GO:0055085); potassium ion transport (two genes, GO:0006813)(Supplementary Table 16).”\nSupplementary Table 16. The enriched pathway of 179 genes with selective signals between the VT group and SS group based on the top 5% Fst and θπ ratio.\nInput Background\n#Term GOterms P-Value\nnumber number\npreribosome,largesubunit\nGO:0030687 3 19 0.000212785\nprecursor\nregulationofheartrate GO:0002027 2 5 0.000585502\nribosomallargesubunitbinding GO:0043023 2 6 0.000777946\noxidativephosphorylation GO:0006119 2 6 0.000777946\nmaturationof5.8SrRNA GO:0000460 2 9 0.001512182\n3'-5'-exoribonucleaseactivity GO:0000175 2 11 0.002129639\nplasmamembrane GO:0005886 10 638 0.002447592\nmaturationofLSU-rRNA GO:0000470 2 12 0.002475932\nmRNAcatabolicprocess GO:0006402 2 14 0.003242292\nmitochondrialelectrontransport,\nGO:0006122 2 14 0.003242292\nubiquinoltocytochromec\nmitochondrialmatrix GO:0005759 3 53 0.003358149\nmitochondrialintermembrane\nGO:0005758 2 15 0.003661824\nspace\ntrachealoutgrowth,opentracheal\nGO:0007426 2 16 0.004105238\nsystem\nmitochondrialelectrontransport,\nGO:0006123 2 16 0.004105238\ncytochromectooxygen\npotassiumiontransport GO:0006813 2 17 0.004572271\npositiveregulationofTORC1\nGO:1904263 2 19 0.005576155\nsignaling\ncationchannelactivity GO:0005261 2 20 0.006112492\nkinaseactivity GO:0016301 2 20 0.006112492\nimaginaldiscdevelopment GO:0007444 2 22 0.007252685\nribosomebinding GO:0043022 2 23 0.007856039\ncytosol GO:0005829 8 547 0.009820574\ngermcellmigration GO:0008354 2 27 0.010485402\ncationtransport GO:0006812 2 27 0.010485402\nsynapticvesicleendocytosis GO:0048488 2 28 0.011195514\nregulationoftranscription,\nGO:0006355 5 249 0.011802892\nDNA-templated\nelectrontransferactivity GO:0009055 2 29 0.011926255\ncalmodulinbinding GO:0005516 2 34 0.015881178\nsleep GO:0030431 2 34 0.015881178\ncytoplasm GO:0005737 14 1386 0.016764745\ntranslationalinitiation GO:0006413 2 36 0.017599416\nnegativeregulationoftranslation GO:0017148 2 36 0.017599416\nchaetadevelopment GO:0022416 2 37 0.018486831\nchitinbinding GO:0008061 3 104 0.019779583\nnegativeregulationof\nGO:0045892 3 105 0.020267449\ntranscription,DNA-templated\nintegralcomponentofmembrane GO:0016021 9 760 0.021975131\ntranslationinitiationfactoractivity GO:0003743 2 41 0.022219975\ntransmembranetransport GO:0055085 4 198 0.023268374\nflavinadeninedinucleotidebinding GO:0050660 2 45 0.026236332\nearlyendosome GO:0005769 2 47 0.028346492\nhemebinding GO:0020037 3 120 0.028379122\nRNApolymeraseIIcis-regulatory\nregionsequence-specificDNA GO:0000978 3 125 0.031412918\nbinding\nLine 186: The duration of salinity stress experiments (24, 48, and 72 hours) seems too short to reflect long-term acclimation. Given that rafting voyages can last for months, a longer experimental period would better represent natural conditions.\nResponse: Thanks for your wonderful suggestions. We added salinity stress experiments for coastal species T. cf. martensii. The salinity gradients included 35, 20, 10, 5 ppt. We first used short-term salinity stress experiments (24, 48, and 72 hours) to show the rapid salinity shifts during rafting process. We describe the plastic expression based on these experiments as “These dynamic expression changes highlight transcriptional plasticity36 (Fig. 4g–l, Supplementary Fig. 19c–f), with the stronger shifts linked to the larger salinity shifts (Fig. 4k) and longer treatments (Fig. 4l).”\nIn addition to salinity stress experiments on 24, 48, 72h, we have performed salinity treatments under 15 days and 30 days in 20 ppt for coastal T. cf. Martensii in our revised manuscripts. Major results as follows:\n(1) Relatively strong expression shifts were detected at low salinities (Fig. 4h) and over the longest treatments (Fig. 4i). Moreover, the number of genes enriched in ion channel transport pathways increased as treatment time increased (Fig. 4j, Supplementary Table 19).\nFig. 4 Genetic changes in T. martensii associated with osmoregulation promote rafting acclimations to salinity shifts. h, Ternary plots depicting the higher proportion of gene expression at 5 ppt than 10 and 20 ppt. i, Ternary plots depicting the proportion of gene expression increased along the duration from 72 h to 15 d and 30 d compared with 24 h and 48 h. j, The number of genes that enriched in ion channel transport pathways along the treatment duration. e.g., R-DME-983712: Ion channel transport; GO:1902600: proton transmembrane transport; details in Supplementary Table 19.\n(2) These dynamic expression changes highlight transcriptional plasticity36 (Fig. 4g–l, Supplementary Fig. 19c–f), with the stronger shifts linked to the larger salinity shifts (Fig. 4k) and longer treatments (Fig. 4l).\nSpecifically, the V-type H+ ATP genes related to proton transport showed decreased expression with increasing treatment duration (VATA, VATB, VATE, VATG, VATO). Notably, genes involved in the transport of small molecules (Tret1, ABCA3, HXT1) and inorganic ion (Sodium/hydrogen exchanger, Sodium/Potassium/Chloride Transporter and KCAB2) showed increased expression in these treatments (Fig. 4l, Supplementary Fig. 20, and Supplementary Table 20).\nFig. 4 Genetic changes in T. martensii associated with osmoregulation promote rafting acclimations to salinity shifts. k and l, Discriminant analysis of principal components (DAPC) for variation in gene expression plasticity in response to salinity shifts (k) and duration changes (l). m, Hypothesized model of ion transport across the gill based on several increased (red) and decreased (blue) expression genes as duration increased.\n(3) Based on these results, we added Discussion as “We found increasing differential expression of genes in response to decreasing salinity, with evidence that gills function as the major osmoregulatory organs in maintaining homeostasis under fluctuating salinity58,59. Rapid gene expression shifts suggest plasticity may be a key prerequisite for successful rafting, at least in the short term60 (Fig. 6c). Nevertheless, theoretical studies suggest that, in the longer term, plasticity may slow or constrain populations’ adaptive genetic response to environmental stress61. Therefore, we propose that both short-term (plastic) responses to temporal environmental variation, and longer-lasting adaptive responses to local conditions, should be integrated to understand how rafting species may disperse and respond to future climates.”\nLines 190-192: The manuscript highlights a few selected differentially expressed genes (DEGs). However, what is the general expression pattern for all DEGs? Cherry-picking without functional validation makes it difficult to draw robust conclusions.\nResponse: In our revised manuscript, we reanalyzed the GEA analysis and found 459 genes showed significant correlation with the salinity. Among these genes, we found 32 genes showed differential expression between high salinity (35ppt) and low salinity (two genes in 5 ppt, 30 genes in 10 ppt). Among these 32 genes, 18 genes (56.25%) show differential expression between the high and low salinity gradients in the terrestrial species M. aosen. Among these 18 genes, 16 genes show the same upregulated or downregulated trends in both species.\nTherefore, according to the expression changes related to salinity changes in two species, we thought talitrids would use a similar genetic strategy to overcome the salinity changes.\nTable. The expression level of 32 genes estimated from GEA analysis show differential expression between high salinity(35ppt) and low salinity in two species.\nCoastal Terrestrial\nTalorchestia martensii Morinoia aosen\ngeneID DEGs GeneFunction\nLow High P Low High P\nSalinity Salinity value Salinity Salinity value\nmartensii_DN644_\n0.33 341.55 0.01 31.69 214.68 0.01 Both Transmembraneprotein45B\nc3_g1_i5_p1\nmartensii_DN3764\n2.67 257.33 0.00 16.22 1597.58 0.00 Both cuticleprotein\n_c0_g2_i1_p1\nmartensii_DN3297\n0.33 59.25 0.01 42.25 255.71 0.00 Both Beta-galactosidase\n_c0_g1_i2_p1\nmartensii_DN3276\n13.33 93.00 0.03 59.74 485.71 0.00 Both PolycombproteinAsx\n_c1_g3_i1_p1\nmartensii_DN1189 Coiled-coildomain-containing\n6.17 33.66 0.01 51.64 277.87 0.00 Both\n4_c0_g1_i45_p1 protein28B\nmartensii_DN1009\n5.33 31.67 0.02 22.17 244.75 0.00 Both mucin-5AC-like\n3_c0_g1_i1_p1\nmartensii_DN3655 AdhesionG-proteincoupled\n116.67 55.67 0.02 509.56 123.94 0.00 Both\n_c0_g1_i1_p1 receptorD1\nmartensii_DN3416\n24.33 77.00 0.04 426.25 969.52 0.00 Both Cyclin-T2\n_c0_g1_i1_p1\nmartensii_DN3234\n0.00 51.71 0.00 109.72 386.27 0.00 Both NA\n_c0_g1_i5_p1\nmartensii_DN515_ 29552.9\n0.00 684.26 0.01 3289.80 0.01 Both Phenoloxidase2\nc0_g1_i3_p1 1\nmartensii_DN4391\n18.78 98.26 0.02 172.21 526.09 0.01 Both Ran-bindingprotein9\n_c0_g1_i1_p2\nmartensii_DN9098 Tectoninbeta-propeller\n51.50 22.17 0.03 200.89 104.83 0.01 Both\n_c0_g1_i1_p1 repeat-containingprotein2\nmartensii_DN451_ TranscriptionfactorwithAP2\n45.71 170.18 0.02 36.46 339.07 0.01 Both\nc6_g1_i3_p1 domain\nmartensii_DN7713\n69.36 269.67 0.00 2718.81 964.09 0.01 Both Kinesin-likeproteinKIF13A\n_c0_g2_i2_p1\nmartensii_DN1135\n22.00 103.33 0.03 94.28 240.09 0.01 Both Dipeptidylpeptidase8\n_c1_g3_i2_p1\nmartensii_DN9608\n25.33 806.67 0.01 73.58 1997.90 0.02 Both Stabilin-2\n_c0_g1_i2_p1\nmartensii_DN1195 SUNdomain-containing\n17.11 63.78 0.05 715.38 2574.65 0.03 Both\n2_c0_g1_i1_p1 ossificationfactor\nmartensii_DN6174\n45.50 183.84 0.00 1682.21 3791.52 0.04 Both Nipped-Bprotein\n_c0_g2_i1_p1\nmartensii_DN8771 Only_ LLGLscribblecellpolarity\n8.44 44.55 0.02 489.01 764.12 0.06\n_c0_g1_i2_p1 martensii complexcomponent2\nmartensii_DN7668 Only_ Probablethiopurine\n421.83 151.59 0.00 10.92 1.94 0.08\n_c0_g1_i1_p1 martensii S-methyltransferase\nmartensii_DN2284 Only_ TELO2-interactingprotein1\n19.22 71.37 0.05 696.03 481.15 0.10\n_c4_g1_i1_p1 martensii homolog\nmartensii_DN3189 Only_ Low-densitylipoprotein\n14.67 91.33 0.00 368.33 602.26 0.16\n2_c1_g1_i1_p1 martensii receptor-relatedprotein5\nmartensii_DN9107 Only_\n251.33 6.00 0.00 9.62 27.19 0.18 NA\n4_c0_g1_i1_p1 martensii\nmartensii_DN1277 Only_ WDrepeatandHMG-box\n24.81 169.54 0.04 209.84 380.96 0.21\n9_c0_g1_i2_p1 martensii DNA-bindingprotein1\nmartensii_DN1430 Only_\n16.00 71.00 0.00 212.75 316.28 0.24 PHDfingerprotein12\n5_c0_g1_i1_p1 martensii\nmartensii_DN361_ Only_\n11.24 162.56 0.00 505.47 779.55 0.41 Focaladhesionkinase1\nc0_g1_i1_p1 martensii\nmartensii_DN1082 Only_\n25.27 8.57 0.00 126.37 173.20 0.43 CarboxylesterasetypeB\n8_c0_g1_i3_p1 martensii\nmartensii_DN3012 Only_\n0.00 19.33 0.04 9.63 19.40 0.51 C-typelectinfold\n_c0_g1_i2_p1 martensii\nmartensii_DN5691 Only_\n47.98 21.42 0.00 24.73 19.64 0.68 Dopaminereceptor2\n_c0_g1_i1_p1 martensii\nmartensii_DN3042 Only_\n33.70 5.57 0.03 164.61 177.24 0.86 Dopaminebeta-hydroxylase\n_c0_g1_i3_p2 martensii\nmartensii_DN1769 Only_\n100.00 311.00 0.04 1030.91 1020.86 0.99 ZincfingerCCCH-type\n5_c0_g1_i1_p1 martensii\nmartensii_DN2484 Only_ plateletbindingprotein\n0.67 366.83 0.04 6.43 10997.77 1.00\n_c0_g1_i2_p1 martensii GspB-like\nLine 215: It would be interesting to see how gene expression changes under longer treatments. This could provide insights into whether plastic responses increase over time.\nResponse: In addition to salinity stress experiments on 24, 48, 72h, we have performed salinity treatments under 15 days and 30 days in 20 ppt for coastal T. cf. Martensii in our revised manuscripts.\nWe have added Results as “To further test for salinity plasticity relevant to long-lived natural rafting voyages, we undertook extended stress experiments lasting for 15 and 30 days. Relatively strong expression shifts were detected at low salinities (Fig. 4h) and over the longest treatments (Fig. 4i). Moreover, the number of genes enriched in ion channel transport pathways increased as treatment time increased (Fig. 4j, Supplementary Table 19). Specifically, the V-type H+ ATP genes related to proton transport showed decreased expression with increasing treatment duration (VATA, VATB, VATE, VATG, VATO). Notably, genes involved in the transport of small molecules (Tret1, ABCA3, HXT1) and inorganic ion (Sodium/hydrogen exchanger, Sodium/Potassium/Chloride Transporter and KCAB2) showed increased expression in these treatments (Fig. 4l, Supplementary Fig. 20, and Supplementary Table 20).”\nWe also added Discussion as “We found increasing differential expression of genes in response to decreasing salinity, with evidence that gills function as the major osmoregulatory organs in maintaining homeostasis under fluctuating salinity58,59. Rapid gene expression shifts suggest plasticity may be a key prerequisite for successful rafting, at least in the short term60 (Fig. 6c). Nevertheless, theoretical studies suggest that, in the longer term, plasticity may slow or constrain populations’ adaptive genetic response to environmental stress61. Therefore, we propose that both short-term (plastic) responses to temporal environmental variation, and longer-lasting adaptive responses to local conditions, should be integrated to understand how rafting species may disperse and respond to future climates.”\nLine 218: Are there significance tests for these genes? Including statistical validation would strengthen the interpretation of results.\nResponse: We reported these genes as differentially expressed genes, which were identified using “DESeqDataSetFromMatrix” and “DESeq” functions in the R package DESeq2. The threshold value of these DEGs was set to P-value < 0.05 and absolute of fold change ≥ 2 (log2|FC| ≥1). Log2(FC) ≥ 1 indicates upregulation, and log2(FC) ≤ -1 indicates downregulation. In the revised manuscript, we used new genes to show gradient changes along the duration gradients and added P-value on the graph.\nSupplementary Fig. 20. Expression changes of salinity-related genes. Gradually up-regulated expression genes (red) involved in small molecular exchange along increased duration. Gradually down-regulated expression of V-type H+-ATPases encoding genes (blue) along with increased duration treatments. Gene expression level was scaled use log (FPKM+1). The P-value means significant expression 2 changes between the short term treatments (24 hours) and long term treatments (30 days).\nLines 332-333: This section seems a bit weak to me. The discussion on wood degradation genes seems narrowly focused on a few candidates. Whole-genome comparisons between rafting and non-rafting species would help identify a broader set of key genes involved in driftwood digestion.\nResponse: (1) To gain a comparison of a broader set of key genes involved in driftwood digestion, we searched the six rafting species (T. martensii, Trinorchestia longiramus, Florersorchestia mkomani, Cochinorchestia sp., Platorchestia pacifica, and Morinoia aosen) and three non-rafting species genomes (Hyalella azteca, Parhyale hawaiensis, and Hirondellea gigas) for Carbohydrate Active Enzymes (CAZymes) belonging to glycosyl hydrolases (GHs), carbohydrate esterases, auxiliary activities, along with associated carbohydrate-binding modules and glycosyl transferases. The CAZyme family composition was largely conserved between rafting and non-rafting species; however, the copy number of CAZyme genes varied notably (Supplementary Fig. 21).\nSupplementary Fig. 21. The proportion and detailed number of CAZymes between rafting and non-rafting species.\n(2) We also compared the number of CAZyme repertoire between rafting and non-rafting species. We found there are 16 GHs, three carbohydrate esterases, and two auxiliary activities exist in all six rafting and three non-rafting species. In particular, rafting species exhibited significantly higher proportions of GH genes compared to non-rafting species (Fig. 5c, Supplementary Table 21).\nFig. 5 Horizontal gene transfer and gene duplication of GH gene families. b, The reconstructed species tree based on three non-rafting and six rafting talitrid genomes. c, Comparisons of abundances of GH gene families in all cazymes suggest rafting talitrids have more GHs than the non-rafting species (Wilcoxon signed-rank test).\n(4) Based on the genome sequence of these rafting and non-rafting species, we also performed gene family expansion analysis, and we found the 89 gene families that were significantly expanded in the rafting species (P<0.05), and GH7 gene family is one of them. Therefore, we thought it is reasonable to propose that GH7 is one of the key genes in the wood feeding of rafting species. We also found a lot of genes participate in the carbon metabolism (22 genes) and the citric acid (TCA) cycle (six genes), which suggest more key genes would be involved in the energy generation and usage.\nSupplementary Table 22. The expanded gene families in the rafting species that estimated using Orthofinder and CAFE based on the protein sequence of nine rafting and non-rafting species.\nNon-rafting Rafting\nID GENE_function Ph Ha Hgi Flo PP T28 Trilo marte tgx\naw zt ga re F 15 ng nsii ia\nputative\nOG0002018 peptidyl-alpha-hydroxyglycine-alpha-ami 1 1 1 3 3 2 4 0 1\ndatinglyaseI\nOG0000006 hypotheticalproteinFHG87_002758 2 0 77 6 4 77 76 16 7\nhypotheticalprotein\nOG0000099 2 0 1 8 28 26 2 3 5\nHAZT_HAZT009215\nhypotheticalprotein\nOG0001831 1 1 1 3 1 1 6 2 1\nFOCC_FOCC013914\nhypotheticalproteinFHG87_024609,\nOG0000212 0 0 5 0 0 20 21 3 1\npartial\nOG0000614 ATPsynthaseF1complexbetasubunit 1 1 1 3 1 11 2 4 4\nOG0000169 solutesymporter 0 2 2 10 9 9 5 9 10\nOG0000036 hypotheticalproteinFHG87_019591 1 0 64 2 0 48 10 0 1\nOG0000072 hypotheticalproteinFHG87_019441 0 0 1 0 0 46 33 0 9\nOG0000294 PDZdomain 3 1 5 16 5 0 6 4 2\nOG0000124 AlcoholdehydrogenaseclassIII 0 8 1 11 8 17 9 5 7\nPyridoxalphosphate-dependent\nOG0000371 1 2 1 4 2 9 3 6 9\ntransferase\nhypotheticalproteinFHG87_018343,\nOG0000032 23 0 48 19 11 21 11 5 1\npartial\nOG0001222 hypotheticalproteinFHG87_018293 0 0 1 2 1 4 2 0 10\nOG0001207 RUNdomain 1 1 1 1 4 2 7 2 1\nOG0000011 WD40-repeat-containingdomain 9 0 10 6 4 71 92 1 1\nOG0000136 reductaseSDR 13 1 1 8 14 11 10 3 1\nOG0000079 hypotheticalproteinFHG87_013911 1 1 4 4 5 19 37 0 15\nOG0000146 hypotheticalproteinFHG87_013667 2 1 2 17 9 3 5 10 12\nOG0000037 FtsZGTPasedomain 10 12 34 13 11 10 12 11 13\nOG0000193 Chitinbindingdomain 2 0 2 16 8 8 6 5 5\nOG0002743 hypotheticalproteinFHG87_012374 1 0 1 2 1 5 3 0 1\nOG0000493 Pacifastindomain 1 2 1 16 1 2 5 1 2\nOG0001087 REJ-likedomain 1 1 1 5 6 1 3 2 1\nOG0000269 CytochromeP450 2 1 2 11 13 10 0 1 4\nOG0001246 RNArecognitionmotifdomain 1 1 1 7 4 2 2 1 1\nOG0000336 polyaminetransporterI 3 3 3 8 8 2 3 6 3\nOG0000263 ZincfingerC2H2-type,partial 1 1 1 18 16 1 4 2 1\nPiggyBactransposableelement-derived\nOG0000137 7 2 2 2 5 10 13 8 13\nprotein\nHistoneacetyltransferasedomain\nOG0001536 1 1 1 5 4 2 2 1 1\nMYST-type\nhypotheticalproteinFHG87_003766,\nOG0000348 2 2 3 7 8 4 2 4 6\npartial\nOG0000089 UDP-glucosyltransferase 4 9 1 16 15 6 8 9 13\nGprotein-coupledreceptor\nOG0001320 1 2 1 3 1 1 3 4 4\nrhodopsin-like\nOG0000050 hypotheticalproteinFHG87_002434 0 0 15 1 1 6 17 59 9\nOG0000187 ImmunoglobulinV-setdomain 1 4 7 1 6 14 5 8 7\nIonotropicglutamatereceptor\nOG0002263 1 0 1 0 1 3 3 3 4\nL-glutamateandglycine-bindingdomain\nOG0000191 Majorfacilitatorsuperfamily 4 1 1 2 13 6 4 4 17\nhypotheticalproteinA3Q56_08581,\nOG0000101 0 0 35 6 0 17 8 2 7\npartial\nOG0000008 proteinSZT2-like 25 6 23 0 2 42 57 50 1\nOG0000009 TransposaseTc1-like 5 0 26 31 3 78 48 9 2\nOG0002104 uncharacterizedproteinLOC108664535 1 1 1 5 1 0 3 3 1\nOG0001353 kinesin-likeproteinKIF21B 1 1 1 2 7 2 3 1 1\npiggyBactransposableelement-derived\nOG0000283 3 3 1 1 2 8 9 4 12\nprotein4-likeisoformX3\nOG0000623 uncharacterizedproteinLOC108670669 1 1 3 12 2 2 3 3 1\nOG0000280 sorbitoldehydrogenase-like 4 3 1 8 8 5 5 5 4\nOG0000416 probablecytochromeP450 2 2 2 6 6 6 3 6 2\nmyosinlightchainkinase,smooth\nOG0000473 2 4 1 5 6 1 3 8 2\nmuscle-like,partial\nOG0002207 cadherin-23-like 1 1 1 1 5 0 3 3 1\nOG0001473 uncharacterizedproteinLOC108675586 1 1 1 1 8 1 2 2 2\nglycoprotein\nOG0000145 0 6 3 20 7 8 6 3 8\n3-alpha-L-fucosyltransferaseA-like\nOG0002531 prickle-likeprotein2 1 1 1 2 4 0 3 2 1\nserine-richadhesinforplatelets-like,\nOG0000586 1 2 3 10 3 1 3 3 3\npartial\nOG0000118 alpha-2-macroglobulin-likeisoformX7 1 9 1 12 14 2 12 3 14\nOG0000458 Glycosidehydrolasefamily7 2 2 1 6 6 3 1 3 8\nOG0000609 proteinMON2homolog 2 1 2 9 2 2 5 4 1\ncentrosomalproteinof89kDaisoform\nOG0000237 5 1 1 1 3 22 8 1 5\nX4\nOG0000186 aldosereductase-like 4 3 1 12 10 2 6 6 9\nOG0002373 uncharacterizedproteinLOC113811614 1 1 1 3 4 0 2 2 1\nOG0000003 Zincfingerprotein 63 9 36 38 39 43 28 23 28\nOG0000012 phosphatase 5 0 16 8 1 62 85 1 0\nOG0000020 Reversetranscriptasedomain 4 4 5 8 0 42 99 0 0\nOG0000021 hypotheticalprotein 7 2 84 20 2 28 9 3 0\nOG0000024 Reversetranscriptasedomain 27 0 19 23 22 44 9 8 0\nGprotein-coupledreceptor\nOG0000030 3 0 4 0 8 93 29 4 0\nrhodopsin-like\nOG0000042 ZincfingerCCHC-type 2 2 3 0 5 14 88 0 0\nOG0000043 hypotheticalprotein 2 0 1 4 1 49 51 5 0\nOG0000056 hypotheticalprotein 0 0 1 4 2 19 73 1 0\nOG0000060 Reversetranscriptasedomain 2 0 2 6 10 23 53 2 0\nOG0000069 Homeoboxdomain-like 0 0 3 10 3 55 11 7 0\nOG0000070 Reversetranscriptasedomain 2 0 16 2 0 39 29 1 0\nReversetranscriptaseRNA-dependent\nOG0000071 4 0 2 2 4 18 59 0 0\nDNApolymerase\nOG0000078 hypotheticalprotein 1 0 10 12 12 17 1 33 0\nOG0000082 hypotheticalprotein 1 0 10 6 2 46 18 1 0\nOG0000083 Integrasecatalyticcore 5 0 3 6 3 17 49 0 0\nOG0000174 RibonucleaseH-likedomain 0 1 10 1 2 13 28 0 0\nOG0000178 RibonucleaseH-likedomain 0 1 4 2 1 13 32 1 0\nOG0000230 Transposasetype1 7 0 1 13 4 8 12 2 0\nOG0000255 hypotheticalprotein 0 0 12 0 11 17 5 1 0\nOG0000264 hypotheticalprotein 0 2 20 4 5 10 4 0 0\nOG0000301 hypotheticalprotein 0 0 3 0 4 24 9 2 0\nOG0000331 hypotheticalprotein 1 0 1 0 4 24 10 0 0\nOG0000525 Actinfamily 3 0 1 2 5 12 4 3 0\nOG0000701 Homeoboxdomain-like 0 0 1 0 1 14 11 0 0\nOG0002607 Cadherin 1 1 1 4 3 2 1 2 0\nOG0002690 PEP-utilizingenzymemobiledomain 1 1 1 2 2 0 3 4 0\nOG0002909 Neurexindomain 1 0 1 4 3 0 2 3 0\nOG0002953 ABCtransporter-like 1 1 1 5 2 0 2 2 0\nOG0002963 FAM91C-terminaldomain 1 1 1 5 2 0 2 2 0\nOG0003733 hypotheticalprotein 1 1 1 3 3 0 2 2 0\nThe discussion is fairly weak, I would expect to see more discussion on how biotic interactions (e.g., competition) affect rafting success; comparisons between salinity adaptation mechanisms in T. martensii and other marine species; the evolutionary history of GRIK2 and its role in other species.\nResponse: We have reconstructed and streamlined the whole discussion. According to your suggestions, major changes are as follows:\n(1) As far as how biotic interactions (e.g., competition) affect rafting success, we have added Discussion as: “In addition to physiological plasticity, competition for space and food on the driftwood also influences rafting success62. The slender morphology of T. martensii (< 15mm) likely facilitates their crypsis on and in driftwood, minimizing potential predation by birds or fish63. Furthermore, their fecundity (10-20 eggs per female) and rapid generation time likely enhance the survival of rafting populations. Crucially, the ability of this species to feed on driftwood, enables lengthy rafting voyages64. In particular, the high proportion of key GH genes that are highly expressed in the gut emphasizes the key role of wood feeding in this rafting system (Fig. 6c).”\n(2) As far as comparisons between salinity adaptation mechanisms in T. martensii and other marine species, we first compared the salinity adaptation gene changes of T. martensii with M. aosen, “The original habitat of T. martensii is marine and high salinity environments (ca. 35 ppt). Compared with the control samples, we found more genes would show differential expression as the salinity decreased to 20, 10 and 5 ppt, it means an increase of environment stress induces broader transcriptional responses (Fig. 4g, h). A similar pattern was found in the Morinoia aosen, a terrestrial species, with more gene expression changes in response to salinity shifts from low to high salinity stress (Supplementary Fig. 19). The salinity experiments from two species provide robust evidence that talitrids would use some plastic osmoregulatory mechanism to cope with the salinity stress.”\nThen, we added Discussion of gills as the osmoregulatory organ similar to other crustaceans. “Thus, we propose that amphipod gills may be crucial osmoregulatory organs, with comparable differential expression pattern between tissues also detected in the terrestrial talitrid, M. aosen (Supplementary Fig. 19b).”\nWe also referenced other marine organisms (Atlantic killifish, Fundulus heteroclitus) that use plastic strategies to adapt to salinity changes. “Rapid gene expression shifts suggest plasticity may be a key prerequisite for successful rafting, at least in the short term60 (Fig. 6c). Nevertheless, theoretical studies suggest that, in the longer term, plasticity may slow or constrain populations’ adaptive genetic response to environmental stress61.”\n(3) As far as the evolutionary history of GRIK2 and its role in other species, we used blastp method to search the 1:1 orthologs in another 11 amphipod genomes, which belong to freshwater, marine, and terrestrial habitats and could be a perfect model for us to compare the gene's role in salinity adaptation. We aligned the protein and reconstructed the phylogenetic relationship between these species.\nWe have added the phylogeny to the supplementary figures and added results description. “We used BLASTp to detect 1:1 orthologs in the genomes of 11 additional amphipod taxa, encompassing a variety of freshwater, marine and terrestrial habitats. Phylogenetic analysis revealed that these GRIK2 sequences clustered strongly by habitat type (freshwater; marine; terrestrial) rather than reflecting the underlying species phylogeny (Supplementary Fig. 18b), providing further support for this gene’s role in amphipod salinity adaptation.”\nLines 348-349: Could the authors quantify the expression dosage of duplicated genes in T. martensii and compare it to species with fewer gene copies?\nResponse: To quantify the expression levels of duplicated genes, we sequenced the transcriptomes of five tissues (gut, gills, legs, brains, and antennae) from three species, Parhyale hawaiensis, Talorchestia martensii, and Morinoia aosen. We compared the GH7 expression between different tissues, and found the GH7 genes show higher expression in their gut system than the gills, legs, brains, and antennae. Besides, we detected similar expression levels between different copies for non-rafting and rafting species in the same duplication events in the same duplication events (Fig. 5f, Supplementary Fig. 22). Considering the larger number of GH7 copies and their expression in the gut, we proposed these duplication events increased gene dosage.\nWe have added Results description as “Moreover, transcriptome analysis revealed that GH7, GH9, and GH30_1 gene copies are highly expressed in the gut, with comparable expression levels in rafting species despite their greater gene copy numbers relative to non-rafting species (Fig. 5f, Supplementary Fig. 22).”\nSupplementary Fig. 22. Gene expression and duplication events within three GH gene families. Pink dots represent duplication events that happened within rafting talitrids, and blue dots represent duplication events that happened at the talitrid ancestor, including marine aquatic talitrids.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nThe study \"Genomics of rafting crustaceans in tropical oceans: adapting to marine climate change\" uses an interdisciplinary approach combining ecogenomics, population genetics, dispersal and oceanography simulations and modeling, and salinity stress experiments. The authors have satisfactorily addressed all comments.\nMy only concern is that the LD threshold used appears too low, and I would expect a justification for this choice in the Methods section. Statistically, an r² ≥ 0.5 is considered sufficient to avoid LD, as values above 0.5 occur in less than 14% of the allelic frequency space (VanLieren & Rosenberg 2008. Theor. Popul. Biol. 74:130). Using a lower threshold can only be justified in specific cases, such as very large sampled populations. A brief explanation of this choice should be added to the text. That aside, I consider the manuscript acceptable for publication.\nResponse: Thank you for your wonderful suggestions. In the Plink software, --ld-window-r2 controls the r2 threshold for a variant-pair to be included in the report. The default setting is 0.2. Although this value might be considered relatively low, when we used a window-sliding based method in 50bp, using an r² threshold of 0.2 helps to reduce the impact of physical linkage while preserving enough independent markers for the following population structure inference. We also employed an r² threshold of 0.5 for population structure analysis, and found that the inferred structure (K=5) was nearly identical to that derived from the two comparison values.\nWe have changed the sentence in the Materials and methods section as “We further excluded highly correlated SNPs by performing linkage disequilibrium based SNP pruning in PLINK v.1.9.078 with a default pairwise r2 threshold 0.2. This threshold reduces LD-induced bias while preserving enough independent genetic markers to accurately infer population structure, an approach previously used in marine genomic studies79.”\nFigure. Admixture analysis based on the SNPs after linkage disequilibrium pruning using threshold 0.2 (a) and 0.5 (b).\nReviewer #1 (Remarks on code availability):\nThe code seems fine. There is a README file with brief explanations of each file.\nResponse: Thank you very much.\nReviewer #2 (Remarks to the Author):\nI am very pleased with the revisions. Thank you, and congratulations on all your hard work.\nResponse: Thank you for your wonderful comments.\nReviewer #2 (Remarks on code availability):\nI looked at the code, but I'm afraid I don't know how to run R code.\nResponse: Thank you very much. The R code could be run after packages are installed, and more information could be obtained by contacting the first author, liuhongguang@ioz.ac.cn.\nReviewer #3 (Remarks to the Author):\nThe revised manuscript addresses many of the concerns raised in the previous review, and is a robust and comprehensive exploration of the genomics of rafting crustaceans. Overall, the integration of interdisciplinary data in this manuscript is a clear strength. Putting all this data together was clearly a lot of work, and the authors are to be commended for the synthesis. With some further clarification of the limitations of the offset analysis, particularly emphasizing how their interpretation is debated, this will make a nice contribution to the literature.\nResponse: Thank you for your wonderful suggestions. We have made corrections according to your suggestions.\nComments:\nFramework: The new figure (Fig. 6) describing the study’s framework better synthesizes the study’s design, making it more easily reproduced in future studies on rafting species. Additionally, the revised sentence in the abstract better reflects what is presented in the manuscript. This resolves the previous comment that the study does not put forward a framework in the text.\nResponse: Thank you for your wonderful suggestions.\nGEA analysis: The updated GEA results are robust and comprehensive, and the added discussion better reflects the limitations of these analyses in rafting species. The GEA results presented here are consistent with the results from other analyses in this study and justified given the recovery of population structure. The authors could draw a parallel between rafting species and other high-dispersal marine species (e.g. Strongylocentrotus purpuratus, Chlorostoma funebralis, & Semibalanus balanoides), which have been shown to exhibit local adaptation despite high dispersal/gene flow.\nResponse: We have added a parallel between rafting species and other high-dispersal marine species in the Discussion section “These results also revealed signatures of local adaptation in rafting species despite their high dispersal potential, which align with previous studies of other high-dispersal marine species, such as Strongylocentrotus purpuratus56, Chlorostoma funebralis57, and Semibalanus balanoides58.”\nGenomic Offset analysis: The revised text on the different types of genomic offsets is helpful and appropriately clarifies the difference between forward, reverse, and local genetic offsets. The addition to the methods on how ocean current and genomic offset analyses are linked is helpful in clarifying how the two were integrated, however the manuscript would benefit from including a brief summary of the information that is only provided in the response to reviewers.\nResponse: We have added a brief summary as “Based on the GF results, we integrated forward genetic offset with ocean current modelling to calculate directional deviation (ΔD) between ocean currents and least-cost dispersal paths, to assess rafting potential of T. martensi to suitable habitats. In the absence of land barriers, the least-cost path was calculated as the Euclidean distance between two sites. By contrast, when this path was interrupted by land barriers, the least-cost path was inferred from oceanographic connectivity.”\nThe new discussion points provide important caveats about these results and I appreciate the call to further validate these methods in marine systems, although a number of claims by the authors need revision before this can be published. First, offsets do not measure or predict “long-term evolutionary responses to broader environmental variation”, nor do they measure “genetic variations required for long-term adaptation of future climate change”. Rather, they are a measure of an environmental distance, weighted by environments that correlate most strongly with genetic structure, and the interpretation is debated (see papers by Matt Fitzpatrick, papers by Lotterhos, especially Lotterhos 2024 Interpretation issues with “genomic vulnerability” arise from conceptual issues in local adaptation and maladaptation, see Rellstab et al 2021: Prospects and limitations of genomic offset in conservation management, see Hoffman et al: Opportunities and challenges in assessing climate change vulnerability through genomics). In addition, whether the finding of strong outlier loci justifies the use of genomic offsets is debatable: whether offsets correlate with fitness changes in new environments depends on patterns of local adaptation in the metapopulation, the fitness landscape, and other biological responses such as plasticity. In addition, despite many other studies predicting genomic offset to future climates, it is important to note that these analyses do not account for adaptive evolution, instead assuming an instantaneous environmental change. Thus, genomic offsets are not necessarily a reflection of long-term genetic adaptation. The manuscript would benefit from some further acknowledgement of these limitations in using genomic offsets to predict responses to future climates.\nResponse: Thank you very much for your constructive comments. We now note more clearly that genetic offsets are a measure of environmental distance weighted by environments that correlate most strongly with genetic structure. It shows the instantaneous degree of maladaptation of a genome when moved from the environment to which it is adapted into a new environment.\nWe have changed the former explanations about “long-term evolutionary responses to broader environmental variation”, “genetic variations required for long-term adaptation of future climate change”, and added some sentences in the Discussion section as “Genetic offsets are measures of environmental ‘distance’ that correlate most strongly with genetic structure, as highlighted by recent analysis of marine invertebrates49, migratory birds50, 51, and fishes52 to evaluate the climate maladaptation. Genetic offsets show the instantaneous degree of maladaptation of a genome when moved to new environments53.”\nLine 362: This would be a good place to also briefly note the limitations of genomic offsets predictions to future climates.\nResponse: we have added some limitations about the genomic offsets as “Genetic offset analysis could be a powerful approach for applied conservation, although their correlation with fitness changes in new environments depends on patterns of local adaptation in the metapopulation, the fitness landscape, and other biological responses such as plasticity61, 62.”\nGene flow and dispersal modeling: The revised methods and discussion clarify the definition of ∆D and clarify that both K = 3 and K = 5 recover the same results in gene flow analyses. The response to reviewers was also helpful in clarifying why dispersal direction would be different from current direction - I suggest adding part of the explanation about islands and land masses changing the least-cost path on the ocean in comparison to the current’s direction into the manuscript. Figure 3d is helpful in clarifying this point, but more clarity in the main text would be helpful, potentially with reference to this figure.\nResponse: Thank you for your constructive comments. We have added more explanation in the Results section “In the absence of land barriers, the least-cost path was calculated as the Euclidean distance between two sites. By contrast, when this path was interrupted by land barriers, the least-cost path was inferred from oceanographic connectivity.”\nLine 181: “Using the GF approach” is unclear. Clarify that you are testing whether the ocean current modelling shows that organisms are/are not able to reach suitable habitat as defined using forward offset.\nResponse: we have changed the sentence as “Based on the GF results, we integrated forward genetic offset with ocean current modelling to calculate directional deviation (ΔD) between ocean currents and least-cost dispersal paths, to assess rafting potential of T. martensi to suitable habitats.”\nPlasticity analysis: The new experiments resolve previous concerns about the applicability of experimental results. Additionally, the addition to the discussion about plasticity’s short-term effects clarifies the caveats of plasticity analyses.\nResponse: Thank you for your comments.\nCode availability: The requested code is now available.\nResponse: Thank you very much.\nReviewer #3 (Remarks on code availability):\nThe requested code is now available.\nResponse: Thank you very much.\nReviewer #4 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: Thank you for your wonderful suggestions. Please check the point-by-point response and the revised manuscript.\nReviewer #4 (Remarks on code availability):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: Thank you for your wonderful suggestions.\nReviewer #5 (Remarks to the Author):\nAll of my comments have been addressed, and I recommend this manuscript for publication.\nResponse: Thank you very much."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The study \"Genomics of rafting crustaceans in tropical oceans: adapting to marine climate change\" uses an interdisciplinary approach combining ecogenomics, population genetics, dispersal and oceanography simulations and modeling, and salinity stress experiments. The authors have satisfactorily addressed all comments. My only concern is that the LD threshold used appears too low, and I would expect a justification for this choice in the Methods section. Statistically, an r² ≥ 0.5 is considered sufficient to avoid LD, as values above 0.5 occur in less than 14% of the allelic frequency space (VanLieren & Rosenberg 2008. Theor. Popul. Biol. 74:130). Using a lower threshold can only be justified in specific cases, such as very large sampled populations. A brief explanation of this choice should be added to the text. That aside, I consider the manuscript acceptable for publication.", "reviewer_2": "I am very pleased with the revisions. Thank you, and congradualations on all your hard work.", "reviewer_3": "The revised manuscript addresses many of the concerns raised in the previous review, and is a robust and comprehensive exploration of the genomics of rafting crustaceans. Overall, the integration of interdisciplinary data in this manuscript is a clear strength. Putting all this data together was clearly a lot of work, and the authors are to be commended for the synthesis. With some further clarification of the limitations of the offset analysis, particularly emphasizing how their interpretation is debated, this will make a nice contribution to the literature.\nComments:\nFramework: The new figure (Fig. 6) describing the study’s framework better synthesizes the study’s design, making it more easily reproduced in future studies on rafting species. Additionally, the revised sentence in the abstract better reflects what is presented in the manuscript. This resolves the previous comment that the study does not put forward a framework in the text.\nGEA analysis: The updated GEA results are robust and comprehensive, and the added discussion better reflects the limitations of these analyses in rafting species. The GEA results presented here are consistent with the results from other analyses in this study and justified given the recovery of population structure. The authors could draw a parallel between rafting species and other high-dispersal marine species (e.g. Strongylocentrotus purpuratus, Chlorostoma funebralis, & Semibalanus balanoides), which have been shown to exhibit local adaptation despite high dispersal/gene flow.\nGenomic Offset analysis: The revised text on the different types of genomic offsets is helpful and appropriately clarifies the difference between forward, reverse, and local genetic offsets. The addition to the methods on how ocean current and genomic offset analyses are linked is helpful in clarifying how the two were integrated, however the manuscript would benefit from including a brief summary of the information that is only provided in the response to reviewers.\nThe new discussion points provide important caveats about these results and I appreciate the call to further validate these methods in marine systems, although a number of claims by the authors need revision before this can be published. First, offsets do not measure or predict “long-term evolutionary responses to broader environmental variation”, nor do they measure “genetic variations required for long-term adaptation of future climate change”. Rather, they are a measure of an environmental distance, weighted by environments that correlate most strongly with genetic structure, and the interpretation is debated (see papers by Matt Fitzpatrick, papers by Lotterhos, especially Lotterhos 2024 Interpretation issues with “genomic vulnerability” arise from conceptual issues in local adaptation and maladaptation, see Rellstab et al 2021: Prospects and limitations of genomic offset in conservation management, see Hoffman et al: Opportunities and challenges in assessing climate change vulnerability through genomics). In addition, whether the finding of strong outlier loci justifies the use of genomic offsets is debatable: whether offsets correlate with fitness changes in new environments depends on patterns of local adaptation in the metapopulation, the fitness landscape, and other biological responses such as plasticity.\nIn addition, despite many other studies predicting genomic offset to future climates, it is important to note that these analyses do not account for adaptive evolution, instead assuming an instantaneous environmental change. Thus, genomic offsets are not necessarily a reflection of long-term genetic adaptation. The manuscript would benefit from some further acknowledgement of these limitations in using genomic offsets to predict responses to future climates.\nLine 362: This would be a good place to also briefly note the limitations of genomic offsets predictions to future climates.\nGene flow and dispersal modeling: The revised methods and discussion clarify the definition of ∆D and clarify that both K = 3 and K = 5 recover the same results in gene flow analyses. The response to reviewers was also helpful in clarifying why dispersal direction would be different from current direction - I suggest adding part of the explanation about islands and land masses changing the least-cost path on the ocean in comparison to the current’s direction into the manuscript. Figure 3d is helpful in clarifying this point, but more clarity in the main text would be helpful, potentially with reference to this figure.\nLine 181: “Using the GF approach” is unclear. Clarify that you are testing whether the ocean current modelling shows that organisms are/are not able to reach suitable habitat as defined using forward offset.\nPlasticity analysis: The new experiments resolve previous concerns about the applicability of experimental results.\nAdditionally, the addition to the discussion about plasticity’s short-term effects clarifies the caveats of plasticity analyses.\nCode availability: The requested code is now available.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "All of my comments have been addressed, and I recommend this manuscript for publication."}, "author_response": "Response: Thank you very much. Please check the point-by-point response and the revised manuscript."}]

Decision letter