Mixed-marker approach suggests maternal philopatry and sex-biased behaviours of narrow sawfish Anoxypristis cuspidata

Sawfishes are a unique group of benthic rays consisting of 5 species separated into 2 genera: Pristis and Anoxypristis (Faria et al. 2013). Characterised by a large toothed rostrum, sawfish are commonly distributed in shallow coastal and estuarine environments of tropical and subtropical regions (Last & Stevens 2009). All sawfish species have undergone significant declines in geographic range and abun-


INTRODUCTION
Sawfishes are a unique group of benthic rays consisting of 5 species separated into 2 genera: Pristis and Anoxypristis (Faria et al. 2013).Characterised by a large toothed rostrum, sawfish are commonly distributed in shallow coastal and estuarine environments of tropical and subtropical regions (Last & Stevens 2009).All sawfish species have undergone significant declines in geographic range and abun-dance due to anthropogenic pressures including fisheries exploitation, pollution (e.g.mining activity) and habitat loss (Simpfendorfer 2000, Seitz & Poulakis 2006, Dulvy et al. 2016).As a result, sawfish are considered the most endangered family within the class Chondrichthyes (Dulvy et al. 2014).The declining global status of sawfish populations has led to all species being listed as either Endangered or Critically Endangered under the International Union for the Conservation of Nature Red List of Threatened Species (Carlson et al. 2013, D'Anastasi et al. 2013, Kyne et al. 2013a,b, Simpfendorfer 2013), included in Appendix I of the Convention on Trade in Endangered Species (Faria et al. 2013) and included in Appendix I and II of the Convention on Migratory Species (CMS) (COP11, report Annex VII; CMS 2014).
Sawfishes are protected nationally in Australia: 3 Pristis species are listed as Vulnerable under the Environmental Protection and Biodiversity Conservation Act (C'th) (1999) and 'no take' restrictions are in place across all northern states.Australia is considered one of the last strongholds for sawfish (Morgan et al. 2011, Dulvy et al. 2016), with 4 of the 5 species found in northern waters (Last & Stevens 2009).Recent reports from Papua New Guinea (PNG), a neighbouring nation to Australia, have indicated the presence of 4 sawfish species (Anoxypristis cuspidata, Pristis pristis, P. zjisron and P. clavata) within its waters (White et al. 2017).It is thought PNG may provide suitable habitat given the unexpected persistence of sawfish in the region.Due to a lack of historic information, the importance of PNG for sawfish has been under-assessed.Investigations into the current status of sawfish in the region are urgently required, given that no protective measures are currently in place.
The narrow sawfish A. cuspidata is considered the most abundant member of the sawfish family in Australian waters (Peverell 2005) despite undergoing sub stantial global declines of > 50% (Dulvy et al. 2014).Their relatively high abundance is often attributed to their life-history characteristics, which include higher growth rates, increased fecundity (sexual maturity at 2−3 yr, annual breeding) and shorter life expectancy (9 yr) than other sawfish species (Dulvy et al. 2016).Narrow sawfish are commonly found inshore, at depths of 10−40 m (Peverell 2005, Last & Stevens 2009).Ontogenetic, sex-biased and seasonal behaviours drive A. cuspidata to occupy varying ecological niches across its lifespan (Peverell 2005, Tobin et al. 2014).Adults are commonly found in offshore waters at depths of ~40 m, while juveniles and pregnant females are found in shallow inshore and estuarine habitats at depths of less than 15 m (Peverell 2005).Commercial net and trawl fishing often pose a risk as the frequency, locality and types of fishing gear used inadvertently catch all sawfish size classes (as by catch) (Peverell 2005, Devitt et al. 2015).Moreover, juvenile and pregnant females are disproportionately vulnerable to habitat loss as a result of coastal anthropogenic pressures, due to a strong reliance on shallow coastal environments as nurseries (Devitt et al. 2015, Adkins et al. 2016, Dulvy et al. 2016).
In recent years, our understanding of sawfish ecology in the Indo-Pacific has advanced considerably (Peverell 2005, Thorburn et al. 2007, D'Anastasi 2010, Morgan et al. 2011, Phillips et al. 2011, Tobin et al. 2014, White et al. 2017).However, there still remains limited information regarding A. cuspidata population structure, sex-biased dispersal behaviours and use of nursery grounds.An absence of conventional tagging studies (due to high post-release mortality associated with tagging) (Devitt et al. 2015) highlights the importance of population genetic analyses to delineate population structure (Simpfendorfer et al. 2016).Available ecological information is primarily for juvenile and sub-adult classes, as they are most commonly captured in fisheries (Peverell 2005, Tobin et al. 2014).However, information on movement and connectivity of both juvenile and adult A. cuspidata remains unknown.Earlier assessments of female-mediated population structure in A. cuspidata in northern Australia, using mitochondrial DNA (mtDNA) control region (CR) sequences, revealed significant structure between the east coast of Australia and both the Gulf of Carpentaria and northwestern Australia (analysis of molecular variance [AMOVA], Φ ST = 0.016, p = < 0.01) (D'Anastasi et al. 2010, B.R.D. unpubl. data), consistent with findings for other sawfish species using the same marker (Phillips et al. 2011).To date, no ecological or genetic studies have included A. cuspidata populations inhabiting nearby PNG waters.
Establishing effective management approaches for sharks and rays can be facilitated through improved understanding of aspects of their biology, including behaviours specific to age and sex (Harry et al. 2011).Additionally, the identification of discrete areas in which specific life stages may occur (e.g.nursery areas) can also help prioritise management efforts (Heupel et al. 2007). For A. cuspidata, it is yet to be confirmed if dispersal behaviours exhibited by males and females differ.These behaviours can be evaluated using genetic markers with different modes of inheritance (mitochondrial versus nuclear DNA) when assessing population structure (Chapman et al. 2015).Furthermore, mitigating the effects of over exploitation posed by anthropogenic activities is para mount to ensuring the survival of sawfish.Re cently, there have been a number of sawfish recovery and conservation plans announced with both global and Australia-wide focus (Harrison & Dulvy 2014, Australian Department of the Environment 2015).Be cause of the endangered status of A. cuspidata, limited knowledge and call for strategic conservation efforts, this study characterised the genetic connectivity and stock structure of A. cuspidata in northern Australia and the Gulf of Papua, PNG.

Sampling
Fin clips from Anoxypristis cuspidata were sourced from scientific sampling efforts and commercial fisheries between 2000 and 2015 from northern Australia and PNG (Fig. 1).Due to the opportunistic nature of collection and the rarity of this species, sample sizes and measurements collected varied between locations.For some samples, size and sex measurements were not recorded and are therefore unavailable.

Laboratory procedures
Genomic DNA (gDNA) was extracted from fin clips using either a modified salting-out method (Sunnucks & Hales 1996) or a Wizard ® SV Genomic DNA purification system (Promega).A portion of the mtDNA CR was amplified following D'Anastasi (2010), while NADH dehydrogenase 4 (ND4) was amplified using primers ND4-F 5'-TGA CTA CCA AAA GCT CAT GTA GAA GC-3' and Leu-Scyliorhinus 5'-CAT AAC TCT TGC TTG GAG TTG CAC CA-3' (Naylor et al. 2005).PCR was performed using 5 µl 10X Buffer (Promega), 0.125 µl Taq polymerase (Promega), 1.0 mM MgCl 2 , 5 µM of each primer, 0.5 µM dNTP and 1 µl (10−30 ng µl −1 ) DNA in 25 µl volumes.PCR conditions for both fragments included an initial denaturation step at 94°C for 60 s, followed by a touchdown protocol with denaturation at 94°C for 60 s, annealing at 54°C for 30 s and extension at 72°C for 90 s across 5 cycles; the next 5 cycles were performed as stated above with an annealing temperature of 52°C and the final 25 cycles were as above, with an annealing temperature of 50°C.After a total of 35 cycles, a final extension step was performed at 72°C for 5 min.PCR reactions were undertaken in either a Bio-Rad C1000 Thermal Cycler (Bio-Rad) or an Applied Biosystem GeneAmp ® PCR system 9700 (Life Technologies, Thermo Fisher Scientific).PCR products were sequenced bi-directionally using BigDye ® Terminator v3.1 cycle sequencing chemistry (Life Technologies) as per the manufacturer's recommendations.Cycle sequenced products were purified using the CleanSEQ kit (Beckman Coulter) and run on an ABI 3130XL AutoDNA sequencer (Life Technologies).Sequences were edited, aligned (using MUSCLE ;Edgar 2004) and trimmed, with the CR and ND4 sequences concatenated using Geneious vR6.1 (Kearse et al. 2012).Aligned and edited partial mtDNA CR sequences for samples collected in northern Australia from the northwest coast, Gulf of Carpentaria and east coast were provided by B.R.D. (n = 164) (GenBank references JQ026198.1,JQ026199.1,JQ026200.1,JQ026201.1,JQ026202.1 and JQ026203.1)(D'Anastasi 2010).
We developed primers flanking core microsatellite repeats in the sequenced enriched sawfish DNA fragments using Primer3 (http://frodo.wi.mit.edu/ cgi-bin/primer3/primer3_www.cgi).Characterisation identified 7 putative microsatellite markers.Of these, forward primers were fluorescently labelled (FAM, NED, PET or VIC) and grouped into 2 multiplex reactions, with 4 and 3 markers per multiplex, respec-tively (Table 1).Optimization included PCR in 10 µl volumes containing 5 µl 2X TYPE-IT master mix (Qiagen), 1 µl multiplexed primers at final concentrations of 10 µM, and 1 µl DNA (15−40 ng).PCR thermal conditions included an initial denaturation step at 95°C for 5 min, a touchdown protocol with 5 cycles of 95°C for 30 s, 57°C for 90 s and 72°C for 30 s, followed by 5 cycles at 55°C for annealing and a final 22 cycles at 53°C for annealing.A final extension was performed at 60°C for 30 min.All PCR reactions were completed using a Bio-Rad C1000 Thermal Cycler.

Data analysis
For mtDNA sequence analyses, the number of haplotypes (H), haplotype diversities (h) and nucleotide diversities (π) were calculated in Arlequin v3.5 (Excoffier & Lischer 2010).A minimum spanning tree (MST) was constructed using output data from Arlequin v3.5 to visualise the spatial distribution of haplotypic diversity.Tests for past population expansions using Fu's F s were conducted in DnaSP v6 (Rozas et al. 2017).Global tests for differentiation, AMOVA calculations and pairwise Φ ST including accompanying p-values for regional locations were computed using the distance matrix model in Arlequin v3.5 (Excoffier & Lischer 2010).
Genetic diversity metrics for the microsatellite markers included mean number of alleles (N a ), allelic richness (A R ), observed and expected heterozygosities (H o , H e ) and inbreeding coefficients (F IS ), which were calculated in the R package diveRsity (Keenan  (Rousset 2008) using the Markov chain algorithm, with a dememorization of 100 000 over 20 batches and 10 000 iterations per batch in order to detect (and exclude) loci not conforming to the underlying population genetic assumptions.The presence of null alleles and effects of allele dropout were evaluated in FreeNA (Chapuis & Estoup 2007).For estimates of population differentiation, samples missing alleles at 3 or more loci were excluded so that each individual had at least 3 loci and each locus had < 5% missing data.AMOVA and estimates of global differentiation using F ST and accompanying p-values were calculated using the haplotype frequency model in Arlequin v3.5.For all pairwise analyses, the initial significance levels (p < 0.05) were adjusted for simultaneous pairwise comparisons using the false discovery rate described in Narum (2006) to reduce the likelihood of Type I errors.

RESULTS
A final mtDNA sequence alignment containing 757 bp (269 bp and 488 bp fragments of CR and ND4, respectively) with a GC content of 0.425 was analysed.Some samples from the northwest coast of Australia failed to amplify using CR and ND4 mtDNA primers (likely due to degradation of samples/DNA); therefore, sample size was reduced for mitochondrial results representing that region.The Gulf of Papua, east coast Australia and Gulf of Carpentaria populations all had greater than 40 individuals, while the northwest coast was represented by 19 individuals (Fig. 1, Table 2).There were 11 polymorphic sites contained in 13 individual haplotypes -4 of which were singletons (Fig. 2).Haplotype diversity ranged from 0.535 on the east coast to 0.789 on the northwest coast, with an average haplotype diversity of 0.658 across all sampled locations (Table 2).
Microsatellite results from the Gulf of Papua were unavailable, as loci did not amplify reliably.It is unknown if unreliable amplification in this population was an artefact of degraded primers.Therefore, micro satellite results presented here are for the 3 Australian locations exclusively.Six of the 7 microsatellite markers amplified reliably (Ancu152 did not amplify) and were used to genotype 123 A. cuspidata individuals from 3 regions in Australia.Of these 6 loci, 5 were suitable for population genetic analysis after an excess of null alleles was found for one locus (Ancu123) (Table 1).The average number of alleles per location (across the suitable 5 loci) was 11. ranged from 9.97 to 10.60 on the northwest coast and the Gulf of Carpentaria, respectively (Table 3).

Population structure and demographic history
Global tests for differentiation and AMOVA analyses, based on mtDNA haplotype frequencies, revealed significant population structuring (Φ ST = 0.082, p = < 0.001).Pairwise population comparisons revealed structuring between all population pairs, except the Gulf of Carpentaria and the northwest coast (Table 4).Past population expansion tests were significant for the east and northwest Australian coast populations (p < 0.03).Negative Fu's F s values were obtained for all populations (Table 2) indicating a greater than expected number of rare haplotypes within each population.
The east coast, Gulf of Carpentaria and northwest coast locations had average sample sizes of 35, 25.4 and 50.4,respectively, across the 5 microsatellite loci (Table 3).No significant population structure was identified using microsatellite loci in AMOVA analyses (F ST = 0.012, p = 1.000), and genetic diversity (H o ) was > 80%.Additionally, genetic variation within pu tative populations (> 98% of variation) exceeded inter-population variation.This suggests homogeneity across samples; thus, no further pairwise comparisons were considered.

DISCUSSION
This study marks the first assessment of female-and malemediated genetic structure of Anoxypristis cuspidata, using a mixed-marker ap proach (mtDNA and nuc lear markers).Our findings of significant mtDNA structuring between all population pairs, excluding the Gulf of Carpentaria and northwest coast, is likely driven by the effects of genetic drift in combination with female-mediated philopatry.In contrast, microsatellite results identified no significant population structuring, which may represent connectivity facilitated by male-biased dispersal across the Austra lian−Papuan region.

Mitochondrial DNA
Elasmobranch mtDNA evolutionary rates are low (Martin et al. 1992), and sawfish mitochondrial diversity is at the lower end of the range for elasmo- Using multiple markers to increase the amount of variability available in this scenario typically enhances power to detect structure (Ryman & Palm 2006).Our approach of analysing concatenated CR and ND4 sequences increased the number of parsimony-informative nucleotide sites available and confirmed findings of genetic structuring between regions based on CR alone (D'Anastasi 2010, B.R.D. unpubl.data), providing greater confidence in our results.
Haplotypes across the Australian−Papuan region sampled show a shallow evolutionary history for A. cuspidata, with low nucleotide diversity and few mutations between haplotypes.In terms of population differentiation, all populations have relatively little divergence and all share the 2 most common haplotypes.The northwest coast and the Gulf of Carpentaria appear to be a single population, while the east coast and the Gulf of Papua show significant differentiation from all populations.
The lack of genetic structuring between the northwest coast and the Gulf of Carpentaria may have multiple explanations, given that structure is detected at smaller spatial scales (PNG and east coast).Gene flow may be ongoing, but this is unlikely given the presence of numerous private haplotypes in each location.Alternatively, genetic structure may be masked due to a number of factors, including: the large geographic scales of regional grouping in the northwest coast sample (which spans biogeographic and genetic breaks -see DiBattista et al. 2017), infrequent migrations between regions over evolutionary time scales, and insufficient time for genetic diffe rentiation between samples with a shared evolutionary history, e.g.following a range expansion from northwest Australia into the previously landlocked Gulf of Carpentaria (see de Bruyn et al. 2004).Additional data collected from finer spatial scales and additional mtDNA markers are required to further assess this.
Our results indicate the presence of matrilineal structure, particularly between the east Australian coast and the Gulf of Papua.The observed structure, characterised by a few common haplotypes shared across regions, numerous private haplotypes and a shallow evolutionary history is likely the combined effect of genetic drift and a bottleneck/founder effect.The effect of genetic drift, where haplotype frequencies change with the random loss or retention of haplotypes through time, is enhanced in small or declining populations such as sawfish (Dulvy et al. 2014).Rare haplotypes are more likely to be lost through time, which can cause genetic diversity reductions (assuming selection, migration or mutations are not occurring) (Lacy 1987, Slatkin, 1987).Another possible factor defining matrilineal population structure in A. cuspidata populations is the presence of founder/bottleneck events that precede range expansions (e.g.into the Gulf of Carpentaria), which would account for observed high haplotype and low nucleotide diversities (e.g. as per Fauvelot et al. 2003).Similar genetic drift and bottlenecks were reported for Pristis sawfish populations from northern Australia and considered to reflect founder events that occurred during range expansions into new habitats during their evolutionary history in northern Australia (Phillips et al. 2017).Genetic drift/bottleneck events as a driver of population structure is supported by the finding of negative Fu's F s values, which indicate the presence of an excess of private mutations that would not be expected if populations had reached equilibrium (Fu 1997).
The behaviour of adult A. cuspidata may also enhance the effects of genetic drift for populations across the region.It is well documented that many female sharks and rays return to the same region to give birth, a phenomenon referred to as parturition site fidelity or maternal natal philopatry (see Chapman et al. 2015 for definitions).This behaviour is likely the result of selection favouring females that return to give birth in areas where previous generations have been most successful (Hueter et al. 2004, Portnoy et al. 2015).Maternal natal philopatry has been identified in numerous other elasmobranchs displaying sex-biased dispersal (Pardini et al. 2001, Keeney et al. 2005, Portnoy et al. 2010, Daly-Engel et al. 2012, Portnoy et al. 2015, Phillips et al. 2016), including P. pristis (based on comparable sample sizes and number of loci [6] to A. cuspidata) (Phillips et al. 2011(Phillips et al. , 2016)).The significant mtDNA values reported here in combination with characteristics such as demographic aggregation (by sex and age) (Peverell 2005, Tobin et al. 2014, Adkins et al. 2016) demonstrate that A. cuspidata have characteristics associated with maternal natal philopatry.However, additional sampling is required (including metadata of size and sex of individuals sampled) to enable separate female and juvenile testing for more accurate estimates of philopatry.

Microsatellites
In contrast to evidence of maternal natal philopatry derived from mtDNA data, microsatellite data failed to detect population structure for A. cuspidata.Together, these findings suggest male-biased dispersal in the presence of female natal philopatry.Whilst we interpret these results with some caution, given the large geographic scales encompassed in regional groups and the limited number of loci used, these findings are consistent with those for P. pristis, a sawfish that is highly reliant on freshwater reaches of estuaries for juvenile survival (Phillips et al. 2016).In contrast, sex-biased dispersal was not evident for 2 other sawfishes, P. zijsron and P. clavata, which, like A. cuspidata, are not reliant on freshwater for juvenile survival, but which had matching population structuring for both sexes based on both mtDNA and microsatellites (Phillips et al. 2011(Phillips et al. , 2016)).
Male-biased dispersal may provide a number of benefits for populations by reducing inbreeding (Pu sey 1987), dispersing local resource competition (Greenwood 1980) and reducing mate competition (Dobson 1982).Phillips et al. (2016Phillips et al. ( , 2017) ) hypothesised that differing patterns in habitat use may explain the observed differences in dispersal bias among Pristis sawfish.However, our findings suggest involvement of other drivers.In particular, different morphological characteristics between sawfish species displaying sex-biased dispersal and those without are noteworthy.A. cuspidata and P. pristis both possess a large ventral lobe on the caudal fin (Faria et al. 2013).In general, fins with a high aspect ratio and distinct ventral lobe, similar to A. cuspidata fins, and to a lesser extent P. pristis fins, are considered indicators of pelagic swimmers that spend limited time on the benthos (Thomson & Simanek 1977).Migration of distances greater than 1000 km is possible given the A. cuspidata caudal fin shape and availability of suitable habitat across northern Australia (Devitt et al. 2015).

CONCLUSIONS
This study has provided important population genetic information using mtDNA and microsatellite markers for Anoxypristis cuspidata; particularly, evidence of possible maternal natal philopatry and sexbiased male dispersal.Regional mtDNA genetic diversity may be prone to erosion through genetic drift (Phillips et al. 2011).However, more maternal and biparental genetic data are required to confirm structure at smaller geographic scales, in order to validate whether A. cuspidata exhibits maternal natal philo patry coupled with paternal dispersal.Based on current evidence, a decline in female numbers in a particular region is unlikely to be replen-ished by re gional migration, suggesting that a regional management approach is required for maternal stocks and their pups.Inter-regional connectivity, identified by microsatellites, has provided the first evidence for male-biased dispersal in A. cuspidata over spatial scales exceeding 1000 km.This requires confirmation with additional data (e.g. more microsatellites or single nucleotide polymorphisms, sampling over finer spatial scales) before management implications can be confidently resolved.If confirmed, cross-jurisdictional management may be required for male stocks.Lastly, assessing genetic connectivity represented in the nuclear genome between Australia and neighbouring countries, such as PNG and Indonesia, should be expanded to better understand A. cuspidata source−sink dynamics and how different conservation approaches between countries are affecting the conservation outcomes for this species.

Fig. 1 .
Fig. 1.Approximate sampling locations and sample sizes per genetic marker for Anoxypristis cuspidata in Australia and Papua New Guinea.mtDNA: mitochondrial DNA; msat: microsatellite

Table 1 .
Motif, sequence, fluorescent dye, annealing temperature (T a ), size range and multiplex for the Anoxypristis cuspidata microsatellite loci screened in population samples et al. 2013).Hardy-Weinberg equilibrium calculations and linkage disequilibrium (LD) was calculated in Genepop 4.0 a Post-optimisation, loci were not deemed suitable for use in population genetic analyses

Table 2 .
Anoxypristis cuspidata mtDNA results (based on concatenated control region [CR] and NADH dehydrogenase 4 [ND4] sequences) with averages across samples shown.n: number of samples, H: number of haplotypes, h: haplotype diversity ± SE; π: nucleotide diversity ± SE, F: Fu's F s .Asterisks indicate a significant p-value (p < 0.03) Fig. 2. Minimum spanning tree representing Anoxypristis cuspidata haplotype distributions across sampled Australian and Papua New Guinean regions.Size of circle represents frequency of individuals belonging to haplotype.Dashes represent the number of mutational changes between each haplotype, i.e. one dash is equal to a single base-pair change between haplotypes