Large effective population size masks population genetic structure in Hirondellea amphipods within the deepest marine ecosystem, the Mariana Trench

The examination of genetic structure in the deep‐ocean hadal zone has focused on divergence between tectonic trenches to understand how environment and geography may drive species divergence and promote endemism. There has been little attempt to examine localized genetic structure within trenches, partly because of logistical challenges associated with sampling at an appropriate scale, and the large effective population sizes of species that can be sampled adequately may mask underlying genetic structure. Here we examine genetic structure in the superabundant amphipod Hirondellea gigas in the Mariana Trench at depths of 8126–10,545 m. RAD sequencing was used to identify 3182 loci containing 43,408 single nucleotide polymorphisms (SNPs) across individuals after stringent pruning of loci to prevent paralogous multicopy genomic regions being erroneously merged. Principal components analysis of SNP genotypes resolved no genetic structure between sampling locations, consistent with a signature of panmixia. However, discriminant analysis of principal components identified divergence between all sites driven by 301 outlier SNPs in 169 loci and significantly associated with latitude and depth. Functional annotation of loci identified differences between singleton loci used in analysis and paralogous loci pruned from the data set and also between outlier and nonoutlier loci, all consistent with hypotheses explaining the role of transposable elements driving genome dynamics. This study challenges the traditional perspective that highly abundant amphipods within a trench form a single panmictic population. We discuss the findings in relation to eco‐evolutionary and ontogenetic processes operating in the deep sea, and highlight key challenges associated with population genetic analysis in nonmodel systems with inherent large effective population sizes and genomes.


| INTRODUC TI ON
The deep-ocean hadal zone remains one of the most poorly understood and least explored ecosystems on Earth. It extends below 6000 m to full ocean depth at ~11,000 m in the Mariana trench, and whilst it accounts for the bottom 45% of the oceans bathymetric range it represents only 2% of the total seafloor area (Jamieson, 2015). Most of the hadal zone is made up of the deepocean trenches that form at subduction zones between converging tectonic plates, and most of these trenches are situated around the Pacific rim where they form a discontinuous chain of discrete topographic seafloor features separated by intervening abyssal plains (Stewart & Jamieson, 2018).
An enduring question has been to what extent the individual trenches are demographically and evolutionarily independent, each shaping the patterns of local intra-and interspecific genetic structure through the combined effects of extreme environmental conditions and their relative geographical isolation (Ritchie et al., 2017b;. The traditional narrative has been that the trenches promote a high level of species endemism for a specialized hadal fauna that is distinct from shallower deep-sea fauna on the abyssal plains by virtue of adaptation to extreme hydrostatic pressures and reduced food availability (Beliaev, 1989;Wolff, 1959Wolff, , 1970. Many taxa show a pattern of phylogeographical structure that is consistent with this view, such as the Lysianassoid amphipod genus Hirondellea where several species can have restricted geographical ranges associated with individual trenches (France, 1993;Lowry & Stoddart, 2010) while others do not .
Even in those taxa with a cosmopolitan distribution across ocean basins, such as the amphipods of the Eurythenes gryllus complex, there are monophyletic morphospecies found within individual trenches that indicate localized evolutionary divergence and potentially incipient speciation (Eustace et al., 2016;Havermans, 2016;Havermans et al., 2013;Havermans & Smetacek, 2018;Ritchie et al., 2015;Weston, Espinosa-Leal, et al., 2021;Weston, Peart, et al., 2021). A similar picture is seen in the genus Paralicella, where high levels of genetic divergence exist between neighbouring Kermadec and New Hebrides trenches based upon microsatellite DNA polymorphisms.
However, an apparent lack of structure for a different Paralicella species across the Pacific Ocean questions a total lack of gene flow among individual trenches (Ritchie et al., 2017b).
Most research on how trenches influence the spatial distribution of genetic diversity has focused on a broad geographical scale between trenches . However, there is also evidence for localized population and ontogenetic structure within trenches (Blankenship et al., 2006;Eustace et al., 2016;Lacey et al., 2016Lacey et al., , 2018Weston, Espinosa-Leal, et al., 2021;Weston, Peart, et al., 2021). Multiple different deep-sea amphipod species across different trenches show the same patterns of ontogenetic stratification, with juveniles found in the shallower depths within their individual species bathymetric distribution and only adults present in the deeper regions. It has been considered that at shallower depths the juveniles will experience reduced predation and intraspecific competition, will gain a nutritional benefit, and the reduced hydrostatic pressures will confer a physiological and metabolic advantage during growth (Blankenship et al., 2006;Hessler et al., 1978;Lacey et al., 2018). How such ontogenetic structure affects population genetic structure through the effects of differential selection and cohort formation across depth distributions remains unknown.
A major impediment to collating the body of data and evidence necessary to generalize about how the hadal trenches shape patterns of intraspecific genetic structure is the sheer logistical challenge of sampling the deep sea at an appropriate geographical scale and for an adequate number of individuals (Jamieson et al., 2010).
Indeed, many of the high-profile studies on deep-sea genetic structure work within a phylogeographical rather than a population genetic framework (Baco et al., 2016;Havermans et al., 2013;Taylor & Roterman, 2017), focus on topographic features such as hydrothermal vents with a proclivity for marked genetic structure driven by natural selection (Vrijenhoek, 2010;Xiao et al., 2020), or confirm hypotheses of wide-ranging dispersal in cosmopolitan species from relatively sparse sampling (France & Kocher, 1996). Advances in autonomous underwater vehicle technology have increased the capability for sampling the deep sea in a more systematic and less opportunistic manner. This can provide much better insight into how stochastic and deterministic microevolutionary forces shape the spatial distribution of genetic diversity over different geographical scales (Jamieson, 2015).
Here we exploit a sampling campaign to the Mariana Trench that yielded multiple samples of the amphipod Hirondellea gigas to provide the first reports of within-species population genetic structure over a small spatial geographical scale within a single hadal trench. H.
gigas (Birstein & Vinogradov, 1955) is the dominant scavenging amphipod species found deeper than 6000 m in the hadal trenches of the NW Pacific Ocean (Eustace et al., 2013;France, 1993). It is readily caught in large numbers using baited traps deployed on deep-sea autonomous lander vehicles, which simulates the natural occurrence of surface-derived carrion falls (Blankenship & Levin, 2007). Multiple individuals were sampled along a transect that encompasses both the overriding and subducting tectonic plates and the trench axis of the Mariana Trench. We utilize double digest restriction siteassociated DNA sequencing (ddRADseq) (Peterson et al., 2012) to genotype individuals at multiple thousand loci across the Hirondellea genome. The propensity for sampling H. gigas from within the NE Pacific Ocean hadal zones makes it an obvious target species to examine within-species population genetic structure in trenches.
However, this also presents two major challenges.
The first challenge is that the ease of attracting H. gigas to baited traps suggests an extremely large population size. Indeed, the number of individuals attracted to bait are typically hundreds to thousands within hours (Gallo et al., 2015;Hessler et al., 1978).
Such hyperabundancy can influence the detection and interpretation of population genetic divergence (Waples & Gaggiotti, 2006).
Population structure in neutral markers is driven primarily by the effects of genetic drift, with effects inversely proportional to population size (Hedrick et al., 1976;Waples, 1998). As such there will be very slow accumulation of genetic divergence between large subpopulations, even in the absence of gene flow. Such so-called "panmictic inertia" causes a disconnect between genetic differentiation as estimated from neutral polymorphisms and the actual demographic connectivity between locations. The use of genome-wide markers generated using RADseq in analysis means that focus can shift away from just examining purely neutral markers influenced by demographic processes, and instead allow examination of patterns of population structure at outlier loci that are highly differentiated relative to neutral expectations reflecting the influence of selection.
The spatial distribution of such adaptive genetic variation can be used to infer local adaptation and be particularly useful for identifying the footprint of genetic structure in marine systems where high levels of gene flow and large population sizes are common (Deagle et al., 2015).
The second challenge is that deep-sea amphipods have extremely large genome sizes (Ritchie et al., 2017a). These can range up to ~34 Gb in Alicellea gigantea, with Hirondellea dubia having an ~5-Gb genome with a large component of repetitive and paralogous DNA (Ritchie et al., 2017a). This makes it especially difficult to ensure that the relatively short DNA sequences generated using approaches such as RADseq are assembled into true orthologous loci from which single nucleotide polymorphism (SNP) genotypes and allele frequencies can be drawn to infer patterns of population genetic structure. The potential for confounding paralogous gene sequence variants has been highlighted in species that have undergone historical genome duplication events, such as the salmonids (Davidson et al., 2010;Waples et al., 2016), and also in species with large genomes (Deagle et al., 2015) especially where the focal organism lacks a reference genome and so RADseq relies on de novo locus assembly (Nadukkalam Ravindran et al., 2018). In the case of deep-sea amphipods, there is no evidence of ancestral genome duplication, but some indications that large genomes are a consequence of transposable element (TE) proliferation which will duplicate genomic regions as a by-product of transposition (Brown & Thatje, 2014;Ritchie et al., 2017a). Both genome duplication and the actions of repetitive elements increase the potential that sequence contigs from polymorphic multicopy genomic regions are erroneously collapsed into a single locus, artificially inflating heterozygosity and biasing allele frequency. Several bioinformatic pipelines are available for RADseq analysis to enrich for true analogues within and between samples by identifying erroneous paralogues based upon optimizing sequence similarity metrics, thus preventing over-splitting of true alleles into separate loci or conversely under-splitting and incorrectly merging paralogous sequences into a single locus (Harvey et al., 2015;McKinney et al., 2017;Nadukkalam Ravindran et al., 2018;O'Leary et al., 2018;Willis et al., 2017). Here we combine several of these approaches at high stringency to identify true independent loci from which proper biological inference can be made. This study therefore presents both new insight into the patterns of population genetic structure within hadal amphipods from the Mariana Trench, but also provides a cautionary tale for population genomics studies on nonmodel organisms with inherently large genomes. Mature male samples were identified according to Barnard and Ingram (1990) and Barnard and Karaman (1991). DNA was extracted according to Ritchie et al. (2017b). All samples were confirmed as from Hirondellea gigas from mitochondrial COI DNA barcodes obtained using LCO1490 and HCO12198 polymerase chain reaction (PCR) primers (Folmer et al., 1994) according to Ritchie et al. (2015).

| ddRAD sequencing, locus assembly and paralogue filtering
Individuals were genotyped using ddRADseq (Peterson et al., 2012) with the two 6-bp restriction enzymes PstI and EcoRI. Samples were multiplexed across four libraries and sequenced across two lanes on an Illumina HiSeq 2000/2500 instrument generating 125bp paired-end reads. Read pairs with evidence of at least 4 bp of 3′ read-through into the EcoRI restriction site plus Illumina adapter (forward reads) or the PstI restriction site plus 8-bp barcode and Illumina adapter (reverse reads) were removed using cutadapt 2.3 (Martin, 2011). The remaining reads were demultiplexed and qualityfiltered (-c -r -q options) using the process_radtags utility of stacks 2.58 (Rochette et al., 2019). Barcodes were rescued allowing one nucleotide mismatch, and read pairs containing residual Illumina adapters with at most two mismatches were discarded.
The sequence data were assembled using the de novo stacks pipeline with additional filtering at various points to enrich for unique nonparalogous loci (hereafter, singletons) ( Figure S1) (Nadukkalam Ravindran et al., 2018). Initial loci were assembled from the forward reads within each individual in ustacks, implementing a strict haplotyping approach to identifying paralogues, where no more than two haplotypes may be observed in an individual and allele merging must be optimized to maximize biallelic loci and minimize monoallelic loci per individual (Ilut et al., 2014;Willis et al., 2017). Alleles were built from primary stacks of at least three identical reads (m = 3) and merged into loci via ungapped alignments allowing up to seven mismatches (M = 7, equivalent to 94% sequence identity) and at most two primary stacks per locus. The deleveraging algorithm was disabled to ensure that overmerged multi-allelic loci originating from diverged paralogous genomic regions were blacklisted. We selected Loci from all individuals were then merged into a catalogue using cstacks in ungapped alignment mode. To prevent paralogues that were diverged between individuals from being merged into the same catalogue locus we required a 100% match (n = 0) between loci from different individuals. This strict approach ensures that paralogue variants are recorded as separate loci but does not preclude bona fide heterozygous catalogue loci to be formed from individuals with biallelic loci. Correspondence between catalogue loci and sample alleles was obtained using sstacks in ungapped alignment mode, and catalogue loci that linked to alleles from multiple loci within an individual were blacklisted since these originate from within-individual paralogues that were missed by the initial locus assembly. The catalogue locus sequences were then clustered by sequence similarity to identify diverged between-individual paralogues (Nadukkalam Ravindran et al., 2018). Sequence clustering was carried out at 40% identity in ungapped end-to-end alignments using vsearch 2.15.1 (Rognes et al., 2016), and all loci that clustered with other loci were blacklisted as putative paralogues. A parameter sweep indicated that the numbers of clusters reached a minimum at 40% (Document S1), suggesting that this degree of identity produces the most conservative set of singleton loci in this highly diverse data set.
Genotyping was completed by incorporating the reverse reads using tsv2bam and calling final genotypes from information across all individuals using gstacks with the Maruki-low model (alpha = 0.05) (Maruki & Lynch, 2017). The final assembled paired-end catalogue loci consensus sequences were examined for microbial and human contamination by querying against all bacterial, viral, fungal and protozoan genomes available in the NCBI REFSEQ database using blastn 2.9.0 (Camacho et al., 2009) with an e-value cut-off of 0.001, and classifying against the PlusPF-16GB database using kraken2 2.0.8 (Wood et al., 2019). Loci identified with at least one of these methods were blacklisted.
(between-individuals methods), plus those blacklisted by haplotyping (within-individuals method). Final genotypes at singleton loci that were common among all five populations and identified in at least 50% of individuals per population were extracted using populations, requiring a minimum minor allele frequency of 0.01 and a maximum observed heterozygosity of 0.6.

| Inference of genetic structure
Genetic structure among individuals and sampling sites was examined using principal components analysis (PCA) and discriminant analysis of principal components (DAPC) in adegenet 2.1.3 (Jombart, 2008). Missing data were replaced with the data setwide mean allele frequency at the corresponding SNP. PCA eigenvalues were used to determine the number of principal components retained for DAPC such that 80% of the variance was preserved. Genetic structure among sampling sites was further examined using pairwise F ST (Nei & Chesser, 1983) and populationspecific F ST (local F ST ) (Weir & Goudet, 2017), both calculated using finepop2 (Kitada et al., 2017(Kitada et al., , 2021. Relationships between individual-based scores on each of the first three linear discriminant functions (DAPC) and latitude or sampling depth were identified using Spearman's rank correlation. The pairwise F ST matrix was summarized using PCA and relationships between the first two principal components and latitude or depth were identified likewise. Isolation by latitude or depth based on the pairwise F ST matrix was tested using Mantel tests with 9999 permutations (vegan package) (Oksanen et al., 2013). Population-specific F ST was regressed onto latitude or depth variables using generalized leastsquares models, accounting for autocorrelation in the response (Kitada et al., 2021).
SNPs driving the observed divergence between sites were identified by hierarchically clustering DAPC variable contributions to the linear discriminant functions into "outlier" and "nonoutlier" SNPs using Ward's method (snpzip function in adegenet). F ST -based outliers were identified using outflank (Lotterhos & Whitlock, 2015) and outliers associated with individual-based structure were identified using multivariate Mahalanobis distances based on the first six principal components using pcadapt (Luu et al., 2017). The union of the three sets was retained as a final set of outliers.
Finally, genetic structure among individuals based on outliers and nonoutliers was inferred using two Bayesian methods. SNPbased inference was carried out by thinning SNPs by linkage disequilibrium using plink 1.9 (Chang et al., 2015) with an r 2 threshold of .5 and then using faststructure with the logistic prior model and inferring between two (K = 2) and five (K = 5) genetic clusters (Raj et al., 2014). Haplotype-based inference was carried out using radpainter and fineradstructure (Malinsky et al., 2018). The co-ancestry matrix inferred by radpainter was inspected using PCA and used to assign individuals to an optimal number of genetic clusters using fineradstructure with a burn-in of 10 6 iterations, 10 6 iterations for Markov chain Monte Carlo (MCMC) and thinning interval of 10 4 . The clustering dendrogram was obtained from the best posterior state identified from at most 10 5 search attempts.

| Functional annotation of loci
The assembled ddRAD loci consensus sequences were functionally annotated to identify functional enrichment among subsets of loci. Loci containing repetitive elements were identified using with an e-value cut-off of 0.001. Enrichment of GO terms among outlier loci against a background of nonoutlier loci was carried out using hypergeometric tests in clusterprofiler 3.14.2 (Yu et al., 2012) with a false discovery rate (FDR)-corrected significance threshold of 0.05. Inspection of locus-sharing patterns among individuals revealed systematic data missingness associated with the ddRAD library, which was partially confounded by sampling site. PCA on Jaccard distances and DAPC between libraries separated individuals from sites C or D in one particular library from individuals from the same site in other libraries ( Figure S2). This same separation was also apparent in SNP genotypes, even in loci shared across all individuals.

| ddRAD data assembly
The library in question had fewer loci overall (Document S1) and more overlap in read pairs (45% vs. 13%-25% pairs merged using flash 1.2.11; Magoč & Salzberg, 2011), suggesting a difference in size-selection regime. To account for this batch effect, loci or SNPs that were among the top 10% features driving the DAPC separation in locus sharing or SNP genotypes between individuals within sites C or D were removed (Deagle et al., 2015;O'Leary et al., 2018). This removed 70% of SNPs with heterozygote deficiency in sites C and D ( Figure S3) and left 824 loci containing 5147 SNPs with minimal residual discrimination between libraries and good discrimination between sites C and D ( Figure S2). Notwithstanding this, all further analyses were also confirmed with a reduced data set, comprising only individuals from sites C, D and E from the outlier library.

| Genetic structure
PCA on SNP genotypes indicated substantial genetic differentiation among individuals but not sampling sites, whereby the two axes accounted for 4.98% and 3.92% of the variance (Figure 2a). In contrast,  (Table S1). Based on DAPC variable loadings, this genetic structure was driven by 301 outlier SNPs. No significant Individual-based Bayesian inference of genetic structure failed to resolve any structure using nonoutlier SNPs and yielded genetic clusters inconsistent with geographical or bathymetric structure using outlier SNPs ( Figure S4). Similarly, the haplotype-based coancestry matrix resolved minimal structure using nonoutlier loci and failed to identify structure congruent with sampling site, geography or bathymetry using outlier loci ( Figure S5). Analysis of a reduced data set comprising only samples from a single library covering sites C, D and E corroborated genetic structure between these three sites ( Figure S6), and local F ST corroborated genetic similarity between the deep sites C and D compared to the shallow site E ( Figure S7). The various filters employed to enrich for singleton loci identified 25,402 loci as paralogues common between at least 50% of individuals from each sampling site-these contained significantly more repetitive elements (5.9% vs. 3.6%; 2 1 = 6.821; p = .009), significantly fewer transcripts (53.0% vs. 57.7%; χ 2 1 = 6.765; p = .009) and similar numbers of proteins (17.0% vs. 18.8%; 2 1 = 1.650; p = .199) compared to the singleton loci. Among these singletons, proportions of annotated loci did not differ significantly between 170 outlier and 654 nonoutlier loci, ranging from 2.4% to 4% for repetitive elements, 52.3% to 59.2% for transcripts and 18.2% to 21.2% for proteins ( Figure 5).

| DISCUSS ION
This study has resolved fine-scale population genetic structure We made several modifications to the standard ddRAD stacks pipeline to ensure that the final data set yielded robust inference of genetic structure. First, we enriched for nonparalogous loci within individuals during ustacks assembly by increasing the number of allowed mismatches between alleles, disabling the deleveraging algorithm and allowing at most two alleles per locus. These settings are different to standard recommendations (Paris et al., 2017), but are effective in blacklisting biologically implausible loci comprising multiple alleles originating from diverged paralogues. We empirically identified the ideal threshold in our data set to strike a balance between minimizing the numbers of monoallelic loci and maximizing the numbers of biallelic loci (Ilut et al., 2014;Willis et al., 2017).
Second, we observed high levels of between-individual loci similarity (Document S1), leading to a decrease in catalogue loci by >40% when allowing seven mismatches in cstacks as suggested by default recommendations (Paris et al., 2017). An effective way of avoiding overmerging of paralogues diverged between individuals into confounded catalogue loci is to undermerge loci in cstacks and blacklist catalogue loci that cluster at relatively low sequence similarity (Nadukkalam Ravindran et al., 2018). In our data set, even a similarity threshold of 90% reduced the number of unique loci across all individuals by more than 50%, and the number of locus clusters did not stabilize until 40% similarity, consistent with extensive locus paralogy (Document S1). Finally, we removed loci with uneven per-allele read depths, which may originate from identical paralogous genomic regions that our previous measures failed to identify (McKinney et al., 2017). Our strict modified RADseq pipeline may be overly conservative, possibly falsely blacklisting diverged singleton loci (Nadukkalam Ravindran et al., 2018), but is clearly essential when dealing with data with extensive paralogy.
Although our focus was not on identifying and characterizing TE content, for which more specialized pipelines exist (Chak & Rubenstein, 2019), comparing the functional annotations of the singleton vs. paralogous locus sets allowed us to gain some insight into F I G U R E 5 Functional annotation of ddRAD loci identified as paralogues, divergence outliers or singleton nonoutliers. Numbers and proportions are given for loci annotated as transcripts (Hirondellea gigas transcriptome), UniProt proteins (arthropoda), repetitive elements (repeatmasker) and the six most frequently observed GeneOntology (biological process) terms in each locus category. Percentages for GO terms are ratios based on total loci with GO annotations. with the singleton loci being associated with GO terms that included regulation of integration, transposition and methylation that would also be involved in policing the activity of TEs. Combined, these observations reinforce the significance of TEs in driving genome size in deep-sea amphipods (Naville et al., 2019;Ritchie et al., 2017a), and are consistent with hypotheses implicating TEs as drivers of evolutionary innovation necessary to allow the colonization of the deep ocean from shallower bathyal depths (Brown & Thatje, 2014). There is growing recognition and focus around how proliferation of TEs can generate genomic novelty through a variety of mechanisms that can lead to rapid adaptation across a broad taxonomic range and under varying ecological contexts (Bourque et al., 2018;Ricci et al., 2018;Wells & Feschotte, 2020). For deep-sea amphipods what is now required is the mapping of different TE families across host genomes to examine location relative to genes that could influence the adaptive potential to extreme hydrostatic pressure, especially for those situations where there is clear gene family expansion or signatures of positive selection (Kobayashi et al., 2018(Kobayashi et al., , 2019Lan et al., 2017).
A further inherent challenge for identifying patterns of population genetic structure in H. gigas is that the high local abundance of these amphipods could show, at most, only subtle levels of genetic divergence (Hedrick et al., 1976). Given that the effect of genetic drift on allele frequencies is inversely proportional to population size, even in the complete absence of gene flow among locations there will only be very slow accumulation of genetic differences at neutral loci, which makes inferring dispersal and population connectivity problematic (Waples & Gaggiotti, 2006). This genomic inertia can cause a disconnect between patterns of demographic and population genetic structure (Lowe & Allendorf, 2010) and can be a common occurrence in marine systems where species can have large effective population sizes and fewer physical barriers to gene flow (Baco et al., 2016;Deagle et al., 2015;Etter et al., 1999). Indeed, our initial data exploration using PCA suggested no genetic structure between sampling locations within the Mariana Trench, consistent with a signature of panmixia. However, a parallel DAPC analysis, which acts to maximize variance between sampling locations whilst minimizing variance within them, identified clear genetic structure across all sampling sites. This differentiation was driven by a relatively small number of SNPs and loci, which display higher levels of divergence among populations than the neutral genomic average and thus can be considered to reflect the effects of adaptive divergence. Markers under divergent selection will change in frequency more readily than neutral polymorphisms even in large populations, though the capacity to identify genomic regions under selection is difficult with low levels of linkage disequilibrium (Pritchard & Przeworski, 2001).
An intriguing finding was that the multivariate structure described by the first two discriminant functions was significantly correlated with both latitude and trench depth. Relationships with latitude may reflect a signature of isolation by distance, though given the structure was driven primarily by a limited set of outlier loci this would also suggest some gradual signature of divergence based upon some environmental gradient. These outlier loci represented GO terms associated with DNA integration, transposition and methylation, in line with the broader suite of singleton loci, again reflecting the potential importance of TE activity and policing in amphipod adaptation and evolution (Brown & Thatje, 2014;Ritchie et al., 2017a).
Critically, the association with trench depth challenges any oversimplistic perspective regarding patterns of dispersal, gene flow and genetic structure in the deep sea and the hadal zone. The traditional paradigm was centred around the depth-differentiation hypothesis (Rex & Etter, 2010) which postulated that there will be decreases in both environmental heterogeneity and barriers to gene flow with increasing depth and hence less structure. Whilst data from some deep-sea species do support the depth-differentiation hypothesis (Cowart et al., 2014;Quattrini et al., 2015;Ritchie et al., 2013) it is apparent that as more studies examine a broader taxonomic spread and range of habitats, the more exceptions there are to the general rule. Disjunct topographical features such as seamounts, fracture zones and ocean ridges do disrupt gene flow, and hydrothermal vents offer localized selection pressures that drive adaptive divergence (Baco et al., 2016). The hadal trenches, as a disconnected network of ultra-deep island-like habitats, are also recognized as drivers of genetic structure in the deep sea. This is evidenced both by the levels of endemism found in the hadal zone, and also from some population genetic studies that have compared across trenches in cosmopolitan amphipod species (Chan et al., 2020;Havermans et al., 2013;Ritchie et al., 2017b;. Patterns of genetic structure have been framed in the context of understanding how the intervening abyssal plains operate as barriers to gene flow and how differences in selection between abyssal and hadal zones may drive divergence.
The current study provides whole new insight into the microgeographical scales over which population genetic structure can accrue, rejecting any concept that highly abundant amphipods within a trench form a single panmictic population. Moreover, it highlights that environmental differences operating within a trench and over microgeographical scales do influence the spatial distribution of genetic variation.
Within-trench population structure has been hinted at before from patterns of ontogenetic variation in amphipods, with a consistent trend across species and trenches of juveniles residing in the shallower depths of the bathymetric range (Blankenship et al., 2006;Eustace et al., 2016;Lacey et al., 2018;Weston, Espinosa-Leal, et al., 2021;Weston, Peart, et al., 2021). Spatially this would mean that juveniles at shallower sites move in a downslope direction as they grow and mature, which could explain the patterns of isolation-by-depth observed in the current data given the environmental gradients associated with trench depth and the effects of genetic drift operating in a similar way to a founder effect as juveniles moved to deeper reaches of the trench. This would be especially pronounced if juveniles formed periodic cohorts. Chronobiological triggers that could entrain reproduction do occur even at full ocean depth (Taira et al., 2004), and temporal and seasonal shifts in surface-derived food supply with associated responses of deep-sea benthic fauna are well documented (Durden et al., 2020;Gooday & Lambshead, 1989;Kalogeropoulou et al., 2010). It is also known that Hirondellea amphipods do not move long distances through the water column, with video footage from the autonomous landers showing they do not move more than 1 m of so from the surface of the trench, and amphipod traps that are suspended above the trench floor rarely catching any individuals (Jamieson et al., 2009). This, coupled with a life cycle in Hirondellea that involves fertilization and juvenile development in a marsupium, would predict that longer distance dispersal is limited. Moreover, the trench topology could provide multiple localized barriers to dispersal, certainly from the trench axis onto the shallower aspects of the plates, and in the case of the Mariana Trench, there are five physically partitioned areas at hadal depths within this one trench .

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors have no conflicts of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw DNA sequence reads are deposited in NCBI SRA bioproject PRJNA834786. Sample genotypes and analysis scripts are available on the GitHub repository: https://github.com/wenze lm/ hgigas_ddRAD.

B EN EFIT-S H A R I N G S TATEM ENT
Benefits Generated: A research collaboration was developed with scientists from the countries providing genetic samples; the results of research have been shared with the provider communities and the broader scientific community. More broadly, our group is committed to international scientific partnerships, as well as institutional capacity building.