Retrospective Genomics Suggests the Disappearance of a Tiger Shark (Galeocerdo Cuvier) Population Off South-Eastern Australia

Over the last century, many populations of sharks have been reduced in numbers by overexploitation or attempts to mitigate human-shark interactions. Still, there is a general perception that populations of large ocean predators cover wide areas and therefore their diversity is less susceptible to local anthropogenic disturbance. Here we report retrospective genomic analyses of DNA using archived and contemporary samples of tiger shark (Galeocerdo cuvier) from eastern Australia. Using SNP loci, we documented a signicant overall change in genetic composition of tiger sharks born over the last century. The change was most likely due to a shift over time in the relative contribution of two well differentiated, but hitherto cryptic populations. Our data strongly indicate a dramatic shift in relative contribution of the two populations to the overall tiger shark abundance of the east coast of Australia, possibly associated with differences in direct or indirect exploitation rates.


Introduction
The tiger shark (Galeocerdo cuvier) is one of the world's largest sharks, with a circumglobal distribution in tropical and sub-tropical waters 16 . Satellite-tag tracking studies have revealed extensive movements of several thousand kilometres both within the Atlantic 17 and Indo-Paci c 18 and, possibly, across the Indian Ocean 19 . This apparent connectivity within basins implies that genetic structuring over these distances may be limited. This has been corroborated by microsatellite genetic studies that show a general lack of intra-basin structuring, but a clear inter-basin genetic split between the Atlantic and Indo-Paci c Oceans [20][21][22] . Accordingly, tiger sharks in Australia are thought to be part of a large Indo-Paci c population. However, population structure of tiger sharks at a local scale remains largely unknown. On the Australian east coast they are found from tropical Queensland to temperate Victoria 23,24 with some satellite-tagged individuals seen to move as far as New Caledonia and Papua New Guinea 24,25 .
Throughout this distribution, the species is exploited through various shing operations, from target and bycatch in commercial, artisanal, and recreational sheries 26,27 , to shark control operations to improve bather safety 28,29 , with clear indications of population decline and reduced average size of individuals [29][30][31] . Tiger sharks are currently listed as Near Threatened on the International Union for the Conservation of Nature's Red List of Threatened Species IUCN (IUCN) due to a suspected decline by ~30% over the past 60 years from exploitation in commercial, recreational, and unregulated sheries, as well as shark control programs 32 . However, further regional depletions have been suggested. For example, tiger sharks in the Arabian Seas region have been assessed as Vulnerable based on a suspected decline of at least 30% over the past three generations 33 . Off Queensland (Australia), tiger shark catch-per-unit-effort has dropped by 74% over the past 25 years, while the average size declined by 21% 31 . Due to its threatened status, potential for over-exploitation and risk to human lives, the tiger shark, like other large sharks, has a high pro le in public awareness, and is subject to much debate over its management and conservation actions 34 .
Here, we investigated a possible link between this apparent decline in tiger shark numbers and its genetic diversity off the east coast of Australia. In order to collect genetic data over the widest spatial and temporal range we extracted DNA from tissue samples taken from shark jaws archived in museums and private collections, or retained as trophies from shing competitions, and performed retrospective genomic analyses 35 . We expected to nd a single panmictic population in eastern Australia as previously reported 20 , but we were alert to the possibility of signi cant changes to the population over the last century as a consequence of documented population reductions.

Results
Bioinformatics pipeline and data ltering Sequencing yielded an average of 2,833,138 reads per individual, with slightly lower numbers for historical jaw samples (median 1,973,493) than contemporary tissue samples (median 2,301,001). After running the bioinformatics pipeline, 78% of the original sequences went into transcriptome mapping. All samples showed a very low percentage of contaminants, con rming the validity of the capture strategy for sequencing enrichment with shark DNA. An average of 0.12% of the cleaned reads was of mitochondrial origin and excluded from further analysis, except for the few cases when it was used for species con rmation. Out of the 20,000 catshark derived baits, 4,544 (22.72%) were post-hoc mapped back to the tiger shark transcriptome 36 covering 4,137 scaffolds. Of these, 4,143 had captured reads with a depth of coverage higher than non-target regions. Scaffolds with a bait had an average coverage of 68.5X, while scaffolds from off-target regions had an overall average depth of 40.7X (35.8X for historical and 43.4X for contemporary samples). Coverage was higher and less variable in contemporary than in historical samples. Single nucleotide polymorphisms (SNPs) were called from all transcriptomic sequences and we identi ed 35,061 raw SNP variants for 122 samples in 2,978 reference scaffolds. After ltering, 4,580 SNPs remained. SNPs with signi cant departures from Hardy-Weinberg Equilibrium (HWE) were removed to produce a nal dataset consisting of 1,840 SNP loci genotyped in 115 samples from G. cuvier specimens caught between 1939 and 2015. In addition, four samples were removed due to high levels of missing data (below 80% call rate) and four other samples were excluded due to either mislabelling of the samples (identi ed as a different species), or to a possible contamination with reads from other shark species. The nal dataset thus contained 107 samples ( Fig. 1, Table S1). Very low levels of DNA damage was observed, con rming DNA was well preserved in jaws 35 , at least for the relatively short time period (~80 years) compared to true ancient DNA (aDNA) studies. Thus, it is highly unlikely that the nal SNPs represent artefacts due to deamination or other DNA damage in the historical jaw samples.

Data analysis for spatial and temporal genomic variability
To examine the stability of patterns over time, we used a temporal genetic analysis based on the backcalculated decade of birth of tiger shark individuals which showed clear evidence of temporal genetic differentiation. This was apparent for the temporal dataset as a whole (Table 1), where pairwise F ST estimates increased with time and with the largest signi cant break occurring between the oldest samples in comparison to the 1990 and 2000 groups. Temporal differentiation was also apparent for the Tasman Sea samples alone, where the majority of the oldest historical samples originated (Table S2) individual (1 of 23) was found (1990). Likewise, only one 'cluster 2' individual (1 of 16) was found among the CRS 2000 samples. When pooling all putative 'cluster 1' and 'cluster 2' individuals across samples, the resulting mean F ST between the two population groups was 0.015, thus substantially higher than between any pair of spatiotemporal samples, supporting a mixed populations hypothesis. The distribution of F ST values across loci showed that a high number of loci (Fig. S1) contributed to the differentiation, suggesting that genetic differences between the two groups were not caused by technical artefacts or contemporary evolution at one or a few loci. However, while bayescan did not detect any outliers, pcadapt identi ed 39 possible outlier loci between the two clusters. Most of these loci (34 of 39) showed an overall higher allele frequency for 'cluster 2'. The mean level of heterozygosity in the two groups was noticeably different (Fig. 5), with 0.23 (± 0.0097) for 'cluster 1' and 0.26 (± 0.0104) for 'cluster 2'. Average missing data was lower for 'cluster 1' than 'cluster 2' (0.04% and 2.3%, respectively), likely re ecting the average age of samples. However, there was no correlation between mean individual heterozygosity and proportion of missing data using a linear model (p = 0.28; R 2 = 0.0021).

Discussion
This is the rst study to demonstrate that archived samples in combination with modern genomic tools can reveal temporal changes in biodiversity, which otherwise would have remained unnoticed. We identi ed genomic differentiation at a relatively local scale for a large migratory shark, the tiger shark However, this pattern was not apparent in contemporary samples. The data analysis of our spatiotemporal samples suggests that there are two distinct population groups of tiger sharks. One of the population groups was most abundant in the oldest and most southern samples, while it was almost completely absent from contemporary (including southern samples) and most northern samples. We propose that the most parsimonious explanation for this absence is either a shift in relative population abundance or population distribution of two previously cryptic populations of tiger sharks, or both. It is possible that the shift was associated with human activities mirroring the ndings from Brown and Roff 37 that reported a major decline in abundance of tiger sharks over three generations off the coast of Queensland with greater declines detected at the southern sites . However, we recognise a number of factors could have affected our observations and there may be other alternative causal explanations for this apparent shift in abundance, which we discuss below.
Our results showed evidence of a temporal change in the genetic composition of tiger sharks on the east coast of Australia over the last century, which are unlikely to be the result of technical issues associated with the use of historical DNA (hDNA). In general, capture sequencing of hDNA has a lower sequencing depth and coverage that leads to fewer mapped reads and thus more missing data with increasing sample age 38 which could affect downstream population genomic inferences [39][40][41] . Although our historical samples presented fewer reads than the contemporary samples, the amount of missing data was generally very low, even for the oldest samples. For example, in 'cluster 2' all except two individuals had less than 5% of missing data (i.e. missing information for 5% of the overall SNP loci). 'Cluster 2' individuals generally showed a higher level of heterozygosity than individuals from 'cluster 1'. This is the opposite of the expected pattern of reduced individual heterozygosity with low coverage caused by allelic drop-out 40 . In addition, the genetic differentiation that we observed between clusters over time was not caused by a few spurious high-differentiated loci (which could be the result of a technical issue in the genomic pipeline), but was found to be spread across the transcriptome. The observed pattern of differentiation was consistent with a scenario of random genetic drift accumulated at an evolutionary time scale in two semi-independent populations of tiger sharks. Thus, the relatively large differentiation between the two putative population clusters, the differences in their levels of heterozygosity, and the nding that not all historical Tasman Sea samples clustered together, renders an alternative hypothesis of a very strong short-term genetic drift in a single panmictic population distributed across the central Indo-Paci c less plausible. In summary, the temporal genetic differentiation was unlikely to be caused by artefacts in sequencing and in the bioinformatics pipeline, by historical genetic drift in a single panmictic population, or by a strong temporal genetic signatures of intra-population environmental selection.
The use of material from historical collections together with genomic-scale analyses provided a relatively large sample size and su cient statistical power for individual-based cluster analyses, presenting a unique window to explore population composition of tiger sharks in the past. The collection of historical and contemporary samples from the east coast of Australia comprised a huge logistic effort and, at least for the oldest specimens, represents the majority of high quality samples available and practically possible to sample in the community to date. However, our samples are unlikely to provide a full picture of the genetic variation across space and time, including population mixing. As sampling was opportunistic, variation in composition of samples with respect to age and sex, as well as sampling time and method, may have in uenced our results. As the oldest samples primarily consisted of large "trophy" sh selectively collected by game-shermen, they are likely to include a higher proportion of mature females than contemporary samples originating from shark control programmes (drumlines or nets) or research projects. Unfortunately, there is little information on sex for the historical samples. However, a large female bias in sex-ratio is apparent in drumline catches for sharks above approximately 280 cm 30 total length. Tagging studies have revealed a variety of behaviours ranging from large-scale movements to considerable site delity 42 ; and if tiger sharks show natal site delity to speci c parturition sites 15 , then samples of large, mature females could also explain a higher occurrence of the local ('cluster 2') population in the Tasman Sea. However, little is still known about the use of space by tiger sharks, and most of the movements and residency studied so far seem to be connected to feeding grounds or other extrinsic factors 15 .
Despite possible sampling-associated uncertainties, the observed patterns of temporal genetic divergence are remarkably clear. Firstly, our results based on contemporary samples are consistent with two recent studies that detected a single genetic population on the Australian east coast 20,21 , and are also compatible with the present understanding of population structure in the Indo-Paci c 21,22 . Secondly, the presence of the two identi ed cryptic groups was not random in space and time. Individuals from 'cluster 2' were mostly detected in the southern part of the species' distribution, most abundant in the oldest samples, and mostly absent from contemporary samples from the Tasman Sea. Consequently we hypothesize tiger sharks in east Australian waters consisted of at least two populations in the past, but likely comprises a single population now. This may sound counterintuitive in the light of satellite tagtracking studies that have shown evidence of individuals migrating over 1,000 km 25,30,43 . However, large migrations and local populations at ner geographical scales are not mutually exclusive, but could be caused by basic "triangle migrations" 15,44 of sh between parturition sites and juvenile and adult habitats 45 . Despite large migrations, tiger sharks have shown evidence of both site delity and residency 42,46 , i.e. they either stay in or return to speci c areas after extensive migrations. This is similar to the white shark, Carcharodon carcharias, which displays large-scale migrations, e.g. Southern Ocean 47 ; Paci c Ocean 48 ; Atlantic Ocean 49 , while also showing ne-scale population structure, e.g. within the Southern Ocean 50 . Little is known as to whether tiger sharks show natal or regional philopatry, as speci c pupping grounds or regions have not been identi ed 51 . For some areas, e.g. Hawaiian Islands, tiger sharks seem to have a different pattern of migration depending on their sex, where females show partial migrations possibly related to reproductive purposes while males do not display such behaviour 46 . Further, the species also displays life history traits that are remarkably different to other elasmobranchs, such as only producing single-sired litters while other species employ multiple paternity as a mating strategy 52 . Females producing litters consisting of pups sired by different males is widely documented to increase the genetic quality of offspring, maintain genetic variation in a population, or increase effective population size 53,54 . The lack of multiple paternity reported for this species may be indicative that tiger shark populations are more vulnerable to the loss of genetic diversity than polyandrous sharks, particularly where overexploitation has resulted in long-term population reductions.
The pattern of population mixing suggests a southern Australian distribution of one of the populations ('cluster 2'), while the apparently more abundant population ('cluster 1') is currently found throughout the entire distribution, with a previous limited intrusion to southern New South Wales (based on the small number of individuals in 'cluster 1' from the historical TAS samples). It is possible that the southerly population ('cluster 2') was a coastal and more resident ecotype and the other population is currently more widespread, more offshore and more migratory, as found in bony shes such as Atlantic cod (Gadus morhua) 55 and European anchovy (Engraulis encrasicolus) 56 and in marine mammals like the bottlenose dolphin (Tursiops truncatus) 57 . Intraspeci c differences in movement and residency patterns is also increasingly being reported in sharks 47 , including in large predatory species (e.g. 47 ). Importantly, our ndings suggest that the current abundance of the putative southerly population has declined compared to pre-1990s levels. The apparent local depletion of tiger sharks at the south eastern distribution of the species in Australia support previous studies showing a reduction in the abundance and mean size of tiger sharks caught 30,31,37 . Although the species is not commonly commercially targeted off eastern Australia, the species is caught as bycatch and in Queensland and New South Wales' shark control programs. Off Queensland, tiger shark catch-per-unit-effort has dropped by 74% over the past 25 years, while the average size declined by 21% 31 , with most of the decline occurring in the southern part of the state 30 . Catch-per-unit-effort also declined in the NSW beach meshing program 29 . The reduction in a tiger shark population, which coincides with the increase in lethal shark mitigation measures, is of concern and may imply ongoing lethal mitigation measures as possible drivers of population declines. Recent estimates of tiger shark catches obtained from commercial logbook records indicate between 5-10 t routinely caught in Queensland, and approximately 3 t per year in New South Wales. However New South Wales also accounts for approximately 10 t of tiger shark recreational catch 58,59 . Further, since 1936 the game shing shery in New South Wales has targeted larger tiger sharks for capture point scores, and continues to do so over several competitions annually 60 . Illegal foreign shing vessels, predominantly Taiwanese vessels that target large sharks for their ns, have been apprehended in the Tasman Sea region and north into the Coral Sea, with tiger shark found to comprise about 20% of the total biomass of shark on the Taiwanese vessels 61 : the largest tiger shark reported to be 442.1 cm total length 61 . Overall, these activities may have selectively removed the southern population ('cluster 2'), most likely through asserting a higher level of exploitation than on the more widespread northern ecotype. Ongoing climate change could also contribute to shifting distributions and abundance of marine sh species and populations 62 , particularly in Australia 63,64 . In particular, since south-eastern Australia has been identi ed as one of the global hotspot for ocean warming 65 . Thus, it is possible that increasing sea temperatures could have negatively affected the putative southern population or enhanced the northern population. Genetic analysis of more contemporary samples from along the east coast of Australia and the Paci c Ocean 20 using e.g. a targeted genotyping approach focusing on the most informative SNPs 66 , could help elucidating the apparent change of abundance of the two populations of tiger sharks in eastern Australian waters.
The apparent occurrence of localized cryptic populations of tiger sharks, at a ner geographical level than hitherto believed, raises a number of concerns regarding identi cation and monitoring of intraspeci c biodiversity in large sharks. Our study suggests that localized populations may be more common than anticipated from recent genetic studies using markers with lower resolution 20,67 .This highlights the importance of development of high-resolution genomic resources for elasmobranchs and other high gene ow marine organisms 68,69 which can provide information on high number of variable sites in the genome (e.g. SNPs). By genotyping many individuals for a high number of SNPs, it may be possible to identify putative populations in species with general low levels of genetic differentiation such as sharks.
More importantly, our work also points to signi cant challenges regarding the scale of current management and biodiversity protection schemes for large sharks. Sustainable management of local populations, through matching the scale of governance with population structure, is important for the protection of the evolutionary legacy of the species and the potential for adapting to future environmental changes 70 . It is also important for the maintenance of healthy marine ecosystems that could provide services to human society 71 . Accordingly, management focus will need to include localized protection measures, such as local seasonal closures or marine reserves to properly match the geographic scale of the population 72 . Speci cally, for the east coast of Australia it should be a priority to further con rm our ndings, elucidate the current abundance and distribution of the two populations and establish measures to protect the putative southern population component, which appears to have faced a signi cant historical decline, driven by either direct and indirect exploitation or environmental change.

Method Sample collection
Tiger shark specimens were caught over a time-span of close to 80 years (1939-2015). Samples originated from north-eastern and eastern Australia, extending from the Gulf of Carpentaria (GCA), through the Coral Sea (CRS) to the Tasman Sea (TAS) (Fig. 1a). Contemporary tissue samples (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) were obtained as n-clips from sharks caught in the Queensland Shark Control Program, the New South Wales Shark Meshing Program, commercial and recreational landings, and sharks caught for tagging and tracking research purposes 73 . Historical samples  comprised dried tiger shark jaws and vertebrae obtained from museum collections, shers, and other private or public collections. The initial dataset consisted of 115 unique sharks. As tiger sharks are long-lived and the sampled individuals were highly variable in age, we estimated the year of birth for each sample to allow for a more accurate temporal genetic comparison (Fig. 1b). For example, a large shark sampled in 2010 could have been born the same year as a smaller shark sampled in 1990. Individual year of birth was estimated using a locally derived relationship between total length (L T ) and age. Growth rate estimation (t) was based on vertebral aging for both females and males using the Von Bertalanffy growth function (VBGF) (1) as; (1) t = ln ((-L T + L ∞ ) / (L ∞ -L 0 )) / -k) 74 , where L 0 and L ∞ represent the length-at-birth and theoretical asymptotic length, respectively, and k represents the growth coe cient. We assumed different parameters for males and females (males: L ∞ = 441.1 cm, k = 0.08, and L 0 = 123.4 cm; females: L ∞ = 379.9 cm, k = 0.06, and L 0 = 116.8 cm), and a combined set where information on sex was not available (L ∞ = 433.7 cm, k = 0.06, and L 0 = 121.5 cm).
For individuals with only fork length (L F ) available, total length was calculated using the relationship L T = 22.607 + 1.096 L F 74 . For the 20 sharks without length information, total weight (W T ) was used to rst obtain L F using a regression equation parameterized for tiger sharks in the north-western Atlantic 75 , the closest population to our target species from which there is available data (L F = ((W T /2.5281) x 10 -6 ) (1/3.2603) ).

DNA extraction and target capture
Historical tissue material was collected following the protocol described in Nielsen et al. 35 and involved the collection of "bio-swarf" produced when drilling a 3.5-mm hole in the calci ed cartilage of jaws or vertebrae. Extraction of DNA from the bio-swarf and contemporary n tissue was performed with the Bioline ISOLATE II Genomic DNA kit according to the manufacturer's protocol, using 18 -37 mg (average 27 mg) of tissue per extraction. For genomic library preparation, DNA from contemporary samples was sheared to an average fragment size of 200 bp with a M220 focused ultrasonicator (Covaris, USA). DNA from historical material was fragmented due to degradation over time and therefore used directly. Genomic-capture libraries were prepared using the KAPA Hyper Prep Kit (Kapa Biosystems, USA) according to manufacturer's instructions. A total of 50 ng input DNA per sample was used with a one in ve dilution of the TruSeq DNA HT dual-index adaptors (Illumina, USA). Ten PCR cycles for library ampli cation were used for the contemporary samples and twelve for the historical. Selected regions of genomic DNA were captured with a MyBaits (MYcroarray) target enrichment kit, consisting of 20.000 biotinylated RNA baits (120 bp each) developed from pancreas, liver and brain derived transcriptome sequences 76 of the small-spotted catshark (Scyliorhinus canicula). At the time of bait development this was the most taxonomically similar species from which a large genomic resource was available for bait design that would likely capture tiger shark transcribed regions. For more details about DNA extraction, library preparation, and bait design see Nielsen et al. 35 . All DNA samples were captured individually, using 135 ng of DNA library as input, which had previously been treated with a 1x AMPure XP beads clean-up. Hybridization capture of tiger shark DNA was conducted for 24 h at 60 °C in solution for subsequent paired-end (2x125 bp) sequencing on an Illumina HiSeq2000 v4. Prior to sequencing, the captured libraries were ampli ed using thirteen PCR cycles and puri ed using 0.8x AMPure XP beads. Quality

Bioinformatics pipeline and data ltering
We customized a bioinformatics pipeline to ensure removal of potential contaminants, artefacts and low quality reads 77 before proceeding with the downstream analysis. Brie y, the de-multiplexed reads were controlled for quality using FastQC 78 . All adaptors were removed using AdapterRemoval 79 and reads were ltered by length and quality, with a minimum length of 30 bp and base quality of 28. Filtered reads were merged using FLASH 80 with default parameters, and checked for contaminants using Kraken2 81 . Both unpaired and concordantly merged reads were mapped against possible sources of contaminants (bacteria, fungi) using the Bowtie2 82 "sensitive" option. The cleaned reads were also mapped against the mitochondrial genome of tiger shark (NCBI Reference Sequence: NC_022193.1 83 ). Previous studies have shown that "off-target capture" is common, i.e. capture of genomic regions of both nuclear and mitochondrial origin not matching the baits. For example, in highly degraded samples target template DNA may not bind and amplify as well as in good quality samples and templates, resulting in more ampli cation of non-targeted regions of the genome 10,11 . Moreover, mtDNA sequences are commonly captured, or directly sequenced, due to the high copy number of mitochondrial DNA compared to nuclear DNA 69,84 . This phenomenon may even be desirable, as it allows assessment of mtDNA diversity. After removal of mitochondrial sequences, the reads were mapped against the transcriptome of tiger shark 36 , using the BWA-mem algorithm. This transcriptome includes 179,867 unique contigs greater than200 bp length (average 822 bp, maximum 15,892 bp). After mapping, PCR duplicates were removed using Picardtools (http://broadinstitute.github.io/picard/). We checked the patterns of DNA damage of the remaining reads, using Mapdamage2.0 85 . Coverage and depth of the target regions were estimated using Samtools 86 , and nally we called SNPs using Freebayes 87 with default parameters. The raw SNPs obtained were further ltered to keep only biallelic SNPs with quality above 30 and minimum allele count of three. Only SNPs with a maximum level of missing data of 20% were maintained. Additionally, we ltered for excess depth to reduce the possible presence of paralogs and multicopy loci. Linkage disequilibrium (LD) between SNPs within bins of 800 bp (maximum length of the merged reads plus 150 bp each side) was estimated using the prune function in bcftools (SAMtools package), by calculating the square correlation between alleles of each pair of loci, r 2 88 and keeping only SNPs with an r 2 < 0.25. To test the reliability of our nal SNPs, we compared genotypes for duplicate control samples and only maintained SNPs that matched more than 80% of the pairwise comparisons (to allow for missing data).
Finally, we ltered for signi cant departure from HWE (p < 0.05) to remove systematic genotyping errors.
All ltering steps but the LD pruning were done using VCFtools 89 .

Data analysis for spatial and temporal genomic variability
Back-calculated year of birth ranged between 1917 (the oldest) and 2012 (the youngest). For all downstream analysis, we grouped samples into four time periods based on their estimated decade of birth: 1910 -1960, 1970 -1980, 1990 94 . An initial PCA of all samples revealed two identical genotypes (two types of archived tissue from the same individual) and one of them subsequently removed. A PCA of the contemporary samples revealed two extreme outliers of which one of them was identi ed as another species (spinner shark; Carcharhinus brevipinna) based on mtDNA sequences. Both samples were removed from further analysis. In order to test for possible differences in individual heterozygosity among population groups, we used an ad-hoc R script, to calculate the proportion of heterozygous loci out of all genotyped loci for an individual, thereby accounting for differences in total number of genotyped loci (i.e. missing data) across individuals. The boxplot was realized using the ggplot2 package in R 95 to highlight the median of each cluster. The distribution of F ST across loci was plotted as a function of heterozygosity using ggplot2. Estimates of Weir and Cockerham 91 F ST and heterozygosity were obtained using VCFtools. The F ST was calculated between individuals belonging to different spatiotemporal sample groups. A test of selection was performed using pcadapt 96  Sampling locations and distribution through time and space. a. Sample distribution along the east-coast of Australia. Samples are grouped by decade of catch (1910-1960, 1970-1980, 1990