First evidence of multiple paternity and hybridisation in Australian sawsharks

Keywords: elasmobranch, genetics, hybrid, mating system, polyandry, pristiophorid, reproductive biology, single nucleotide polymorphisms.

The reproductive biology of cartilaginous fishes (chondrichthyans) is complex. Two general modes of reproduction are observed in this group: oviparity (egg-laying) and viviparity (live-bearing). However, there is much variation within these modes, including the length of time that eggs are retained within the female, whether sperm can be stored, and the level of maternal provisioning provided to the embryos during development (e.g. histotrophy, oophagy, intrauterine cannibalism, placentotrophy; Wourms 1977; Awruch 2015; Nakaya et al. 2020). Investigating the reproductive biology of this group is difficult, largely because their habitat and depth ranges often prevent direct observations of breeding behaviour. However, with the development of molecular techniques, it is now possible to gain new insights into the reproductive biology of chondrichthyans (Portnoy and Heist 2012). For example, molecular analysis has demonstrated the occurrence of parthenogenesis in both captive and wild individuals in a range of species (Chapman et al. 2007; Fields et al. 2015; Dudgeon et al. 2017; Feldheim et al. 2017).

Investigating the presence of multiple paternity is another area where molecular techniques are valuable. Also known as genetic polyandry, multiple paternity refers to the siring of a single litter by more than one male. Potential benefits of multiple paternity include avoiding inbreeding or genetically incompatible mates, increased fecundity, and creating genetically diverse offspring (Jennions and Petrie 2000; Neff and Pitcher 2005; Slatyer et al. 2012). In elasmobranchs, multiple paternity has been documented in both viviparous (Rossouw et al. 2016) and oviparous species (Chevolot et al. 2007; Griffiths et al. 2012; Hook et al. 2019) and appears to be common within the group. Indeed, multiple paternity is believed to be the ancestral condition of all elasmobranchs (Lamarca et al. 2020).

However, this generalisation may not apply to all species. Some elasmobranch species do not appear to exhibit multiple paternity even if its occurrence is widespread within the family. Within the family Carcharhinidae, for example, 11 of the 13 species examined thus far have been shown to exhibit multiple paternity (Bester-van der Merwe et al. 2019; Lamarca et al. 2020; Nash et al. 2021; Armada-Tapia et al. 2023), with the two exceptions being the tiger shark [Galeocerdo cuvier (Holmes et al. 2018; Pirog et al. 2020)] and Galapagos shark [Carcharhinus galapagensis (Daly-Engel et al. 2006)]. However, analysis of C. galapagensis comes with the caveat that only one litter was available for study (Daly-Engel et al. 2006). Furthermore, not all orders have been examined (Lamarca et al. 2020), with many groups such as rays (Torres et al. 2022) and deep-sea sharks (Duchatelet et al. 2020; Nehmens et al. 2020) being under-represented in the literature or having yet to be investigated. Addressing this knowledge gap in a wider variety of species will enable a more complete understanding of the ubiquity of multiple paternity among elasmobranchs.

Once the presence of multiple paternity has been established, molecular techniques also allow the frequency at which it occurs to be determined. When multiple litters have been analysed, results indicate variation in the level of multiple paternity within and among species. Some species exhibit low levels of polyandry (<35%), with the majority of litters being sired by a single male (genetic monogamy; e.g. Chapman et al. 2004; Daly-Engel et al. 2010; Boomer et al. 2013; Duchatelet et al. 2020), whereas in other species singly sired litters are relatively rare (>85% of litters sired by more than one male; e.g. Feldheim et al. 2004; Griffiths et al. 2012; Lyons et al. 2017). Furthermore, rates of multiple paternity have been shown to differ within species among years and geographic locations (Nosal et al. 2013; Chabot and Haggin 2014; Barker et al. 2019b). Therefore, the use of molecular tools to study multiple paternity can help uncover complexities in the mating systems in different species of elasmobranch.

Molecular tools can also provide insights into hybridisation between species. Hybridisation refers to the creation of viable offspring following the mating of two taxa and has been estimated to occur in at least 25% and 10% of plant and animal species, respectively (Mallet 2005). Although hybridisation has long been recognised and is common among marine and freshwater fishes (Scribner et al. 2001; Montanari et al. 2016), its discovery in chondrichthyans only occurred in the last decade. In fishes, hybrid individuals often display colour patterns or morphologies intermediate to the parental species and are therefore relatively easy to identify (Montanari et al. 2016; Tea et al. 2020). However, chondrichthyans often have conserved morphology between species and this makes the recognition of hybrids more difficult. It was only through the examination of the nuclear DNA of common (Carcharhinus limbatus) and Australian (Carcharhinus tilstoni) blacktip sharks, coupled with mismatches in morphology and mitochondrial DNA, that Morgan et al. (2012) provided the first evidence of hybridisation in chondrichthyans. Since then, molecular analysis has assisted in identifying hybridisation in several other species of sharks and rays (Table 1), and suggested the possibility of hybridisation in others, such as sleeper sharks [Somniosus microcephalus and Somniosus pacificus (Walter et al. 2017)], three species of shyshark [Haploblepharus fuscus, Haploblepharus pictus and Haploblepharus edwardsii (van Staden et al. 2020)] and skates in the Mediterranean [Raja polystigma and Raja montagui (Frodella et al. 2016)].

Table 1.  Genetically confirmed reports of hybridisation (first- and second-generation or backcross) in sharks and rays.

Sawsharks (family Pristiophoridae) are a distinctive group of demersal sharks readily identified by their tapered saw-like rostrum and pair of long barbels on the ventral surface of the rostrum (Nevatte and Williamson 2020). The group currently consists of 10 species that inhabit coastal to deep-water marine environments on the continental shelves and slopes throughout the world. Sawsharks have generally received very little attention in the scientific literature and, as such, much of the knowledge of their biology is limited (Ducatez 2019; Nevatte and Williamson 2020).

In Australia, two species of sawshark have overlapping distributions. The common sawshark (Pristiophorus cirratus) and southern sawshark (Pristiophorus nudipinnis) co-occur in the waters of the southern half of Australia, from Western Australia to southern New South Wales and Tasmania (Last and Stevens 2009; Bartes and Braccini 2021), with the distribution of P. cirratus extending up to southern Queensland on Australia’s eastern coast (Last and Stevens 2009; Nevatte et al. 2019). The two species have overlapping depth ranges (5–>600 m; Raoult et al. 2020) and are frequent by-catch in the commercial fisheries operating in the region (Walker et al. 2005; Braccini et al. 2012; Raoult et al. 2020). Similarities in morphology often lead to the two species being grouped together in fisheries records (Emery et al. 2019; Raoult et al. 2020). Harvesting of P. cirratus and P. nudipinnis is currently considered sustainable, with relatively stable catch rates over the past two decades (Raoult et al. 2020; Patterson et al. 2022), and both species are listed on the IUCN Red List as Least Concern (Walker 2016, 2021). Many aspects of the biology of these species remain poorly known, with detailed investigations having occurred only in recent years. These include studies of age and growth (Raoult et al. 2017; Burke et al. 2020b), trophic ecology and possible feeding behaviour (Raoult et al. 2015; Nevatte et al. 2017b; Burke and Williamson 2021), movement and population structure (Burke et al. 2020a; Nevatte et al. 2021), and sensory systems (Nevatte et al. 2017a; Wueringer et al. 2021).

The reproductive biology of P. cirratus and P. nudipinnis is another area where knowledge is limited. At present, information on reproductive biology is based on investigations of specimens caught during fisheries research. As with all sawsharks, P. cirratus and P. nudipinnis have been identified as aplacental viviparous (ovoviviparous), with a biennial reproductive cycle. Litter sizes can range from 3 to 22 pups for P. cirratus and from 7 to 14 pups for P. nudipinnis, with a gestation period of at least 12–16 months (Hudson et al. 2005; Walker and Hudson 2005; Last and Stevens 2009). The mating system of sawsharks has not been directly examined, but given the ubiquity of multiple paternity in elasmobranchs, it is possible that these sharks also exhibit this reproductive strategy. Incorporating information about the mating system of harvested species is important for developing appropriate management plans (Rowe and Hutchings 2003). For example, knowledge of whether a species exhibits monogamy or polyandry can be included in demographic models to assess population declines more accurately (Tsai et al. 2015; Liu et al. 2020). Thus, determining the mating system of P. cirratus and P. nudipinnis would provide information necessary for accurate assessment of their resilience to current levels of fishing pressure.

In this study, we provide a preliminary assessment of multiple paternity in P. cirratus and P. nudipinnis by genotyping pups from a single litter of each species with nuclear single-nucleotide polymorphisms (SNPs). Although a small sample size, multiple paternity has been successfully detected in other shark species where only one litter was examined (Daly-Engel et al. 2006; Larson et al. 2011; Corrigan et al. 2015). We also investigate the possibility of hybridisation between the two sawshark species by genotyping a P. nudipinnis sample that is suspected to have been misidentified. This sample was collected in north-eastern Tasmania and identified in the field as P. nudipinnis. However, sequencing of the mitochondrial DNA showed that this individual possessed a mitochondrial haplotype characteristic of P. cirratus (Nevatte et al. 2021).

Tissue samples (fin clips and muscle tissue stored in 70–90% ethanol) for assessment of multiple paternity and hybridisation were sourced from previous studies. This included samples of a pregnant P. cirratus caught off Wollongong, New South Wales in 2014, and her eight pups (early to mid-term embryos; Fig. 1a) (Nevatte et al. 2017a) and samples of a pregnant P. nudipinnis caught off Kangaroo Island, South Australia in 2017, and her two pups (mid- to late-term embryos; Fig. 1b). The mother and pup P. nudipinnis samples were donated by collaborating researchers as part of a population genetics study on sawsharks (Nevatte et al. 2021). The potentially misidentified P. nudipinnis sample was collected by the authors during research trawling off north-eastern Tasmania in 2011 (Raoult et al. 2015, 2017). Research trawling was conducted under an ethics permit issued by the University of Tasmania (see Raoult et al. (2015, 2017) for permit number). Ethics approval for the other studies (Nevatte et al. 2017a, 2021) was not required because the authors were not involved in the capture of the sawsharks. Tissue samples were collected from deceased sawsharks caught during routine commercial fishing operations. These sawsharks were either purchased directly from or donated by commercial fishers and fish cooperatives.

Fig. 1.  Developmental stages of the (a) common sawshark (Pristiophorus cirratus) and (b) southern sawshark (Pristiophorus nudipinnis) pups examined in this study. For scale, a 30 cm ruler is presented in (a) and a tape measure in centimetres (bottom half) is presented in (b). Photo credits: Ryan Nevatte (a); Matt McMillan (b).

Subsamples of tissue (~20 mg) from the females, their pups and the potentially misidentified individual were sent to Diversity Arrays Technology Pty Ltd (DArT; Canberra, ACT, Australia) for DNA extraction and SNP genotyping with the DArTSeq™ protocol (Sansaloni et al. 2011; Kilian et al. 2012; Cruz et al. 2013). These samples were included in the genotyping assay performed on a set of P. cirratus (n = 55) and P. nudipinnis (n = 26) samples sourced from Nevatte et al. (2021) as part of a concurrent population genetics study (R. J. Nevatte, M. R. Gillings and J. E. Williamson, unpubl. data). All samples were processed following the methodology described in Kilian et al. (2012) and Melville et al. (2017) using the restriction enzymes PstI and HpaII.

The dataset received from DArT was imported into R (https://www.r-project.org; R Core Team 2022) as a genlight object by using the dartR package (ver. 1.9.9.1, https://cran.r-project.org/web/packages/dartR; Gruber et al. 2018) for further filtering and analysis. The RStudio interface (https://posit.co/products/open-source/rstudio/; Posit Team 2023) was used for data management and visualisation of graphical output.

For the assessment of multiple paternity, mothers and their associated pups were filtered alongside the individuals of their respective species. The dataset was first reduced to the species of interest and then locus metrics were recalculated using the gl.recalc.metrics function in dartR. The filtering pipeline for both species then proceeded as follows: (1) removal of monomorphic loci; (2) removal of loci of <98% reproducibility; (3) retaining only one SNP per 69 bp fragment to minimise linkage (selected with the ‘best’ option in the gl.filter.secondaries function); (4) retaining loci with a read depth between 5 and 50 (to prevent poorly genotyped or paralogous sequences in the dataset); (5) removal of loci with a call rate of <90%; (6) retaining individuals with a call rate of >90% (and removing any subsequent monomorphic loci); and (7) removal of loci with a minor allele count of ≤3 [identified with the radiator (ver. 1.2.2 R package, https://thierrygosselin.github.io/radiator/; Gosselin 2020)].

An additional filtering step was implemented to identify and remove loci deviating from expectations of Hardy–Weinberg equilibrium (HWE). All individuals were pooled together and each locus was tested for HWE by using the R package pegas (ver. 1.1, https://cran.r-project.org/web/packages/pegas; Paradis 2010) with 1000 permutations. A Bonferroni correction was then applied to the obtained P-values for each locus by using the p.adjust function in R, and loci with adjusted P-values of <0.05 were removed from the dataset. Finally, loci with missing data were excluded from the dataset, so that all individuals were represented by the same number of loci. This resulted in a final dataset of 1720 and 493 loci for P. cirratus and P. nudipinnis, respectively.

The presence of multiple paternity in the litters of each species was assessed using two methods: pairwise relatedness and sibship analysis. Pairwise relatedness between the mothers and pups was calculated using COANCESTRY (ver. 1.0.1.10, https://www.zsl.org/about-zsl/resources/software/coancestry; Wang 2011). To determine which of the seven different estimators available in the program was most suitable for our data, we simulated 300 dyads each for several relationship categories, including unrelated (r = 0), parent–offspring (r = 0.5), full-siblings (r = 0.5) and half-siblings (r = 0.25), based on the allele frequencies of the final filtered datasets for the two sawshark species. Allele frequency information was extracted from the datasets by using related (ver. 1.0, https://github.com/timothyfrasier/related; Pew et al. 2015), an R implementation of COANCESTRY. The simulation included a conservative 0.01 genotyping error rate for each locus and no missing data or allelic dropout. Pairwise relatedness was then calculated for the simulated dyads with the seven estimators, with the default value (100) for the number of reference individuals used for the triadic likelihood estimator (TrioML). The best relatedness estimator was identified based on both its accuracy (closeness to the true value) and precision (variation around the estimated values) as suggested by Attard et al. (2018). This involved examining the Pearson’s correlation coefficient calculated for each estimator by COANCESTRY and using the simulated dyad data for each relationship category to calculate means and standard deviations, and to construct box plots to visually assess variance. The two maximum-likelihood estimators, TrioML and dyadic likelihood (DyadML), showed the highest correlation to the true value and were the least variable across relationship categories for both sawshark species (Supplementary material Table S1 and Figs S1, S2). Thus, these estimators were selected for the empirical analysis. Parameters for the empirical analysis were the same as for the simulation (i.e. allele frequencies from final filtered datasets, 0.01 genotyping error, no missing data or allelic dropout, 100 reference individuals for TrioML), with the genotyping error accounted for in the analysis. The genotype input file containing all individuals was generated with related, while any other required files were created following the instructions in the COANCESTRY manual.

To determine the most likely number of sires for each litter, a sibship analysis was performed using COLONY (ver. 2.0.6.8, https://www.zsl.org/about-zsl/resources/software/colony; Jones and Wang 2010). The program uses a full-likelihood approach to infer relationships among offspring (full siblings, half siblings or unrelated) and can reconstruct the genotypes of potential parents. For the analyses, COLONY assumes that the loci are in HWE and linkage equilibrium. While the former had been accounted for with a previous filtering step, it was not possible to guarantee complete linkage equilibrium of the markers used here because of the lack of a reference genome for any species of sawshark or other closely related shark species. However, the algorithm implemented in the program can accommodate some level of linkage in the markers (Wang and Santure 2009) and any tightly linked loci were removed during the filtering of the dataset (i.e. only a single SNP was retained when more than one SNP was present in a sequence fragment). Hence, we acknowledge that although there may be some linkage among the markers, this should not affect the results.

The genotypes of the pups and their mothers were included in the sibship analysis with COLONY, with this maternal relationship being specified in the project file. Each species was tested in separate analyses by using the following parameters. Males and females were set to a polygamous mating system, with no inbreeding and no clones present among the offspring. Because sharks are diploid organisms and dioecious (have distinct male and female individuals), these options were also selected. The full-likelihood method was chosen as the analysis method, with a medium likelihood precision, medium run length and a total of five runs. The starting random number seed for the first of these runs was set to 1234. No sibship prior was applied, sibship scaling was set to the default (yes) and allele frequencies were not updated. Markers were set to codominant, with an allelic dropout rate of 0 and a genotype error rate of 0.01 (as for the COANCESTRY analysis). Because the accuracy of sibship reconstruction in COLONY is not greatly affected by the genotype error rate (Ackerman et al. 2017), the error rate selected for the analysis was based on previous studies examining sibship in sharks (Barker et al. 2019b; Reid-Anderson et al. 2019). Allele frequencies were estimated during the analysis.

When determining full-sibling families, COLONY reports two likelihood probabilities to indicate whether the families may have been over- or under-split during the analysis: the inclusion probability and exclusion probability. The inclusion probability shows the probability that all individuals listed in a full-sibling family are true full-siblings and no further splitting into half-siblings is necessary, with a high value indicating that the family is not under-split. The exclusion probability shows the probability that all individuals listed in the full-sibling family are full-siblings and no other individuals should be considered as full-siblings to the individuals in this group. The family is considered not over-split with a high value for this metric. When these likelihood probabilities are used in conjunction, they can show a potential overestimation of the number of sires. For example, a full-sibling family with a high inclusion probability and low exclusion probability indicates that individuals assigned as full-siblings are likely to be real, but the family may have been over-split, and some individuals listed as half-siblings to this family may actually be full-siblings. Thus, the number of sires may be overestimated.

Finally, to assess the strength of the SNP markers in identifying individuals, the probability of identity (probability that two individuals possess an identical genotype) was calculated in GenAlEx (ver. 6.503, https://biology-assets.anu.edu.au/GenAlEx/Welcome.html; Peakall and Smouse 2006, 2012). Because of the presence of related individuals, the more conservative probability of identity with siblings (PISibs; Waits et al. 2001) was calculated for each species by using the same datasets for the COLONY analyses (i.e. only the genotypes of the mothers and her pups).

All P. cirratus (n = 56) and P. nudipinnis (n = 27) samples (excluding the pups) were filtered together with the potentially misidentified individual in dartR. After removing the pups from the dataset and recalculating locus metrics, the filtering pipeline consisted of the following steps: (1) removal of loci with <99% reproducibility; (2) retaining only one SNP per 69 bp fragment to minimise the inclusion of linked loci (selected with the ‘best’ option by using the gl.filter.secondaries function); (3) removal of loci with a read depth of <5 and >50 (to prevent retaining poorly genotyped loci or paralogous sequences); (4) removal of loci with a call rate of <95%; (5) retaining individuals with a call rate of >95%; (6) removing monomorphic loci; and (7) retaining loci with a minor allele frequency of >5%. These threshold values were selected following the inspection of the output provided by the report functions available for each filtering criteria in dartR. A total of 4277 SNPs remained post-filtering, with 82 individuals in the dataset [P. cirratus (n = 54); P. nudipinnis (n = 27); misidentified individual (n = 1)].

Following the initial filtering, a principal-component analysis (PCA) was performed in dartR to visualise the positioning of the potentially misidentified sample in relation to the samples of the two sawshark species. Inspection of the resulting plot (see Results) indicated that the potentially misidentified individual could be a hybrid, and so further analyses were conducted to further explore this hypothesis.

To identify SNPs that were fixed between the two sawshark species, and therefore diagnostic for each species, a dataset with the putative hybrid removed was created and F_ST values for each locus were calculated with the gl.basic.stats function in dartR. Loci that displayed F_ST values of 1 were then selected as diagnostic markers for the sawsharks, resulting in a total of 1348 SNPs (the majority of loci had F_ST values of >0.9; 3039 of 4227). The dataset containing the putative hybrid was then reduced to these SNPs. Finally, this dataset was filtered further to remove SNPs that contained missing data. Following this, a total of 759 SNPs remained.

An admixture analysis was performed on the dataset of 759 SNPs, by using the R package LEA (ver. 3.0.0, https://bioconductor.org/packages/release/bioc/html/LEA.html; Frichot et al. 2014; Frichot and François 2015) to estimate the admixture coefficients (Q) of the putative hybrid and the two sawshark species. The number of ancestral populations (K) that best explained the data was determined based on the calculation and plotting of the cross-entropy criterion, with 100 repetitions for each tested value of K (range 1–10). The cross-entropy plot showed a ‘knee’ at K = 2 (Fig. S3) and was thus selected as the value of K for the analysis. This value also made the most sense biologically because there were two recognised species present in the dataset and the putative hybrid had P. cirratus mitochondrial DNA. Bar plots of the best run for this value of K (i.e. the run with the lowest cross-entropy criterion) were then constructed to visualise admixture.

For a quantitative assessment of possible hybridisation, the SNP data were also analysed with NewHybrids (ver. 1.1, http://ib.berkeley.edu/labs/slatkin/eriq/software/software.htm#NewHybs; Anderson and Thompson 2002). The program utilises a Bayesian clustering algorithm to determine the posterior probability of an individual belonging to several different classes, including pure species, first- or second-generation hybrid and backcross with one species. Because of memory constraints, the program can only analyse approximately 200 loci; so, a random subset of 200 loci were selected from the 759-SNP dataset by using the gl.subset.loci function in dartR. This dataset was then exported from R as a STRUCTURE file and converted into a NewHybrids file with PGDSpider (ver. 2.1.1.5, http://cmpg.unibe.ch/software/PGDSpider; Lischer and Excoffier 2012). The NewHybrids analysis was run in parallel with EasyParallel (https://github.com/hzz0024/EasyParallel; Zhao et al. 2020) and consisted of five runs with Jeffreys-like priors, a burn-in period of 10 000 and 100 000 sweeps. Parental individuals were not specified a priori. All five runs were inspected for consistency of results.

Pairwise relatedness and sibship analysis showed the presence of multiple paternity in both sawshark species. As the relatedness estimates produced by COANCESTRY for the TrioML and DyadML estimators were near identical for both species, only the TrioML estimates are reported here. For P. cirratus, one pup in the litter (Pup 3) failed to produce a library during the SNP genotyping and thus only seven pups were analysed for this species. Nevertheless, pairwise relatedness showed that the seven-pup litter consisted of both full-siblings (11 pairs; TrioML range: 0.4012–0.5902) and half-siblings (10 pairs; TrioML range: 0.2277–0.3404), with Pups 1 and 4 being half-siblings to their litter mates (Fig. 2a). Pairwise estimates for the mother and her pups were consistent with a parent–offspring relationship (TrioML range: 0.4768–0.5378; Fig. 2a). Sibship analysis with COLONY indicated that the litter was sired by two males, thereby corresponding to two full-sibling families. Both families had high inclusion and exclusion probabilities (>0.99), with Father 1 siring two pups and Father 2 siring five pups (paternal skew = 2:5) (Table 2). The pups identified as half-siblings to the rest of their litter mates (Pups 1 and 4) were consistent with the result obtained with the pairwise-relatedness estimates. The calculated PISibs for the dataset of 1720 loci was 1.4 × 10⁻¹²³.

Fig. 2.  Pairwise relatedness estimates (TrioML) between pups and their mothers for (a) common sawshark (Pristiophorus cirratus) and (b) southern sawshark (Prisitophorus nudipinnis). Dark grey cells indicate non-applicable values.

Table 2.  Results from the sibship analysis performed in COLONY for the single litters of Pristiophorus cirratus and Pristiophorus nudipinnis.

For P. nudipinnis, pairwise relatedness in COANCESTRY showed that the two pups in the litter were half-siblings (TrioML estimate: 0.3028; Fig. 2b). A parent–offspring relationship between the mother and the two pups was also confirmed (Fig. 2b). The sibship analysis in COLONY supported a half-sibling relationship between the pups with two full-sibling families being constructed (Table 2). Thus, the litter was sired by two males (paternal skew = 1:1). The calculated PISibs for the dataset of 493 loci was 1.3 × 10⁻³³.

Principal-component analysis suggested the presence of a hybrid. The PCA plot (Fig. 3) showed a clear separation of the samples into two distinct genetic groups by the first axis (PCA Axis 1), which explained 90.2% of the variation. These two groups corresponded to the samples identified as P. cirratus and P. nudipinnis. The suspected misidentified sample was situated between these two groups (Fig. 3), which is consistent with what would be expected for a hybrid individual.

Fig. 3.  Principal-component analysis based on 4277 SNPs to identify a possible hybrid between two sawshark species. Each dot represents an individual sawshark.

The admixture and NewHybrids analyses further supported the hypothesis of a hybrid. The putative hybrid had Q-values of 0.57 and 0.43 for the ancestral populations representing P. nudipinnis and P. cirratus, respectively, while all other individuals had Q-values of 0.99 for one of these populations (Fig. 4). Such admixture proportions are indicative of a first-generation (F₁) hybrid. All five runs of NewHybrids provided identical results, with the putative hybrid being assigned as an F₁ hybrid with a posterior probability of 1. All samples of P. cirratus and P. nudipinnis were also assigned as separate parental species with posterior probabilities of 1. Thus, the three analyses performed in this study indicate that the sample previously identified as P. nudipinnis is an F₁ hybrid.

Fig. 4.  Plot of the admixture coefficients for each individual calculated in the LEA R package based on 759 SNPs. Each bar represents an individual sawshark, with the colours (brown and grey) corresponding to the respective species. The putative hybrid is indicated by the red arrow.

The results of this study deliver new insights into the reproductive biology of Australian sawsharks. Through the application of nuclear-derived SNPs, we provide the first evidence of multiple paternity occurring in both P. cirratus and P. nudipinnis and show that the two species are capable of hybridisation.

Detection of multiple paternity in this study indicates that the mating strategies of P. cirratus and P. nudipinnis involve polyandry. However, given that only a single litter was analysed for each species, it is not possible to determine the frequency of polyandry. For a more complete assessment of the mating system of P. cirratus and P. nudipinnis, additional sampling of pregnant females and their pups is required. Given that rates of multiple paternity can differ among populations (Chabot and Haggin 2014; Barker et al. 2019b), examination of this in genetically differentiated sawshark populations would also be prudent.

As stated previously, the one limitation of this study is the analysis of only single litters for each species. This is particularly an issue for P. nudipinnis, which had a litter size (n = 2) much lower than the reported litter size for this species (7–14 pups; Walker and Hudson 2005; Last and Stevens 2009). This reduced sample size could be due to the female aborting most of the litter following capture, with only two pups remaining within the uterus when samples were collected. Capture-induced parturition is common in elasmobranchs (Adams et al. 2018) and has previously been reported for sawsharks (Bass et al. 1975; Hudson et al. 2005; Walker and Hudson 2005).

Although further research is required to understand multiple paternity in sawsharks, some possible reasons for its presence are discussed. As a reproductive strategy, multiple paternity has been proposed to offer several benefits, including the generation of genetically diverse offspring with greater fitness, avoidance of inbreeding or genetically incompatible mates, and increased fecundity (Jennions and Petrie 2000; Neff and Pitcher 2005; Slatyer et al. 2012). However, studies examining multiple paternity in sharks have not found any evidence of these benefits (Feldheim et al. 2004; DiBattista et al. 2008; Daly-Engel et al. 2010; Boomer et al. 2013). For example, lemon sharks (Negaprion brevirostris) from singly-sired and multiply-sired litters show no difference in survival rates (DiBattista et al. 2008) and fecundity does not increase with multiple matings in smooth-hound sharks (Mustelus spp.) (Boomer et al. 2013; Marino et al. 2015).

Instead, the prevailing hypothesis for the occurrence of multiple paternity is convenience polyandry, whereby females will accept additional or superfluous matings from males if the costs of doing so are lower than the costs of resisting mating attempts (Portnoy et al. 2007; DiBattista et al. 2008; Lyons et al. 2021). Mating in elasmobranchs can often cause harm to the female [e.g. bite wounds from males (Pratt and Carrier 2001; Ritter and Amin 2019)] and this may also be true for sawsharks. While their mating behaviour has not been observed, the toothed rostrum used by sawsharks for sensing, capturing, and manipulating prey (Nevatte et al. 2017b; Burke and Williamson 2021; Wueringer et al. 2021) could inflict injuries when copulating. Although obvious mating wounds have not been observed on the two sawshark species assessed here, female and male sawfishes [a group of rays that also possess a toothed rostrum of similar function (Wueringer et al. 2012)] have been found to bear wounds attributed to the rostrum as a result of their mating habits (Papastamatiou et al. 2015; Brame et al. 2019). Thus, accepting superfluous matings from amorous males may be the least costly course of action for sawsharks.

Although convenience polyandry may contribute to multiple paternity, it should not be considered the sole reason for its occurrence, as other biological factors may also play a role (Lyons et al. 2021). For example, sperm competition, whereby sperm from different males compete for the fertilisation of the ova, is also likely to influence polyandry (Schlegel et al. 2012), and is a known issue in elasmobranchs (Fitzpatrick et al. 2012). Indeed, in male sharks, the morphology of the sperm (flagella length) and testes size have been shown to increase with an increasing incidence of multiple paternity (Rowley et al. 2019a, 2019b). Sperm storage could also play a role as females may exert post-copulatory control on the siring of their litters through the storage and selection of sperm from favoured males (Fitzpatrick et al. 2012). In species that store sperm, this may also increase the amount of sperm competition (Lyons et al. 2021). Although the ability to store sperm has been recorded for many chondrichthyans (Dutilloy and Dunn 2020), this has not been investigated in sawsharks.

Another factor may be a combination of ovulation patterns and embryonic diapause. In round stingrays (Urobatis halleri), differences in the sizes of the embryos and their positioning within the uterus suggest that females release their eggs in a staggered manner (Lyons et al. 2017). This is hypothesised to allow females to control the paternity of their litters (Lyons et al. 2017). While it is unknown whether female sawsharks release all their eggs at once or a few at a time, a staggered ovulation pattern over the breeding season could allow different males the opportunity to sire pups within a single litter. In the sawshark litters examined here, there does not appear to be any difference in the developmental stage of the embryos (e.g. Fig. 1a), which could indicate that females ovulate all their eggs at once. However, it has been suggested that female P. cirratus and P. nudipinnis are capable of embryonic diapause (Hudson et al. 2005; Walker and Hudson 2005; Waltrick et al. 2012), whereby the development of the embryo is suspended temporarily after fertilisation. Embryonic diapause would allow females to retain fertile eggs from multiple mating events and could also mask the effects of a staggered ovulation pattern. Embryos sired by different males or during different mating events would develop simultaneously and therefore have no discernible differences in size. Thus, the combination of staggered ovulation and embryonic diapause could explain the detection of multiple paternity in this study.

The level of multiple paternity may also be influenced by the rate at which potential mates are encountered (the encounter-rate hypothesis) (Daly-Engel et al. 2010; Nash et al. 2021). Breeding aggregations are likely to increase the chances of encounters between mature males and females and, therefore, lead to higher frequencies of multiple paternity. The aggregative behaviour displayed by the common smooth-hound (Mustelus mustelus) has been attributed to its relatively high frequency of polyandry (Rossouw et al. 2016), whereas the low frequency of multiple paternity observed in the Hawaiian spurdog (Squalus hawaiiensis), formerly the shortspine spurdog [Squalus mitsukurii (Daly-Engel et al. 2018)], has been attributed to its asynchronous reproductive cycle and absence of regular mating aggregations (Daly-Engel et al. 2010). However, in mating aggregations of finetooth sharks (Carcharhinus isodon), which display both annual and biennial reproduction, no significant difference in the rate of multiple paternity was detected between sharks using either reproductive mode (Nash et al. 2021). Sharks reproducing biennially were expected to have a higher incidence of multiple paternity because of the potential for more copulations over a single reproductive event.

The numbers of males and females within these breeding aggregations may also affect the frequency of polyandry. For example, the overall frequency of multiple paternity is low in leopard sharks (Triakis semifasciata) from La Jolla, California, which form aggregations consisting predominantly of females (Nosal et al. 2013). However, the frequency of multiple paternity increased in the year in which there was an increase in the number of males (Nosal et al. 2013). Male-biased aggregations are thought to increase multiple-paternity rates owing to the potential for mobbing behaviour around females, although this may be species-specific. Rates of multiple paternity in rig (Mustelus lenticulatus), a species that displays heavily male-biased aggregations, are only slightly higher than those in the closely related gummy shark (Mustelus antarcticus), which is not known to form breeding aggregations (42% M. lenticulatus; 31% M. antarcticus) (Boomer et al. 2013).

Regarding sawsharks, it is unknown whether P. cirratus or P. nudipinnis have mating grounds, but seasonal changes in catch rates suggest that these sharks may migrate into shallow waters for breeding in the austral autumn (Raoult et al. 2020). It is also unknown whether there is a bias towards one sex in these hypothesised breeding aggregations. Identifying whether breeding aggregations occur, and their composition, will be important for furthering our understanding of multiple paternity in these species.

Further information about the hybrid sawshark can be gleaned from molecular and morphometric data obtained during previous studies. For example, the mitochondrial DNA data reveal the maternal species involved because of the maternal inheritance of this marker. Because the hybrid possessed P. cirratus mitochondrial DNA [haplotypes PC Cyt-b 11 (GenBank Accession: MT376141) and PC ND5 7 (GenBank Accession: MT376162) identified in Nevatte et al. (2021)], it is clear that it was the result of a mating between a female P. cirratus and a male P. nudipinnis. The morphometrics showed that the hybrid female had grown to a total length (TL) of 1092 mm, a size consistent with that of adult sawsharks of either species [P. cirratus max. TL = 1490 mm; P. nudipinnis max. TL = 1240 mm (Last and Stevens 2009)]. This indicates that hybrid offspring are capable of surviving to adulthood in the wild. Furthermore, examination of the reproductive organs suggested that the hybrid had reached sexual maturity (max. ovarian follicle diameter = 5 mm; females are considered mature by Walker and Hudson (2005) if max. follicle diameter >3 mm). However, it was not possible to determine whether this female was fertile. The discovery of backcrossed individuals or second-generation hybrids would be required to resolve this question.

Detection of hybridisation in these sawsharks raises questions regarding the rarity of such events and how/why it may have occurred. Previous studies documenting hybridisation in elasmobranchs have generally reported a small number of hybrid individuals (1–4; Table 1), suggesting that the phenomenon might be rare. However, in some species, such as blacktip (Morgan et al. 2012) and hammerhead (Barker et al. 2019a; Barker et al. 2021) sharks, hybridisation appears to be quite extensive (Table 1). Given the number of samples examined in this study, it is suggested that hybridisation between P. cirratus and P. nudipinnis might be rare (only one hybrid was detected in the sample of 82 individuals: ~1%). However, because this study is among the first to examine the genetics of sawsharks in the region, further sampling may reveal the presence of additional hybrids. A possible hybrid sawshark may have been detected in the DNA barcoding studies conducted by Ward et al. (2005, 2008), where one P. nudipinnis sample possessed the mitochondrial haplotype of P. cirratus. Misidentification or mislabelling was suggested as the reason for this because hybridisation had yet to be documented in sharks, but no further investigation is possible because the sample was not retained as a voucher specimen (Ward et al. 2008).

Several factors have been proposed as drivers of hybridisation. The rarity of conspecifics is one such factor (Wirtz 1999; Montanari et al. 2016) and is based on the premise that matings between heterospecifics are more likely to occur when the number of conspecifics in an area is reduced. In particular, females of the rarer species are thought to engage in interspecific matings more than males (Wirtz 1999). Such unidirectional hybridisation has been documented in hammerhead sharks, with the rarer Carolina hammerhead (Sphyrna gilberti) often being the maternal species when hybridising with the more common scalloped hammerhead [Sphyrna lewini; (Barker et al. 2019a)]. However, in this study, the opposite was found. Fisheries data indicate that P. cirratus is more abundant than P. nudipinnis across their overlapping depth range (Raoult et al. 2020), and P. cirratus was detected as the maternal species for the hybrid. Thus, it is more likely that a general lack of conspecifics at the time of breeding may have contributed to the hybridisation. However, further research is needed to determine whether hybridisation is unidirectional or bidirectional between the two sawshark species.

Another factor proposed as a driver of hybridisation is niche overlap, which includes an overlap in habitat use and diet (Montanari et al. 2016). Although both species have overlapping depth ranges and occupy the same types of habitat, stable isotope analysis has shown that P. cirratus and P. nudipinnis display resource partitioning in regions where they co-occur (Raoult et al. 2015). Thus, overlap in diet is unlikely to be a major driver of hybridisation between these species. The breeding behaviour of either species is unknown. As such, it is uncertain whether P. cirratus and P. nudipinnis breed in the same areas or at similar times, although changes in catch rates suggest seasonal movements of these sharks (Raoult et al. 2020). The sharing of breeding habitat could facilitate hybridisation events.

This study has shown that P. cirratus and P. nudipinnis exhibit multiple paternity and that hybridisation between the two species can result in viable hybrids that can survive to reproductive maturity. While preliminary, the results provide important new information on the reproductive biology of the two sawshark species caught in Australian commercial fisheries. Further research in this area, including the collection of additional mother and pup samples and identifying potential mating grounds, is required to better understand the reproductive biology of P. cirratus and P. nudipinnis. Also, the combined use of mitochondrial and nuclear DNA markers in future molecular studies on these sawsharks will be necessary to ensure that specimens are correctly identified to species and any potential hybrids are detected.

Supplementary material is available online.

The data that support this study will be shared upon reasonable request to the corresponding author.

The authors declare no conflicts of interest.

Funding for the SNP genotyping was provided by a Save Our Seas Foundation Keystone Grant (#2019/441) awarded to J. E. Williamson, and the School of Natural Sciences at Macquarie University. R. J. Nevatte was supported by a Macquarie University Research Training Program (MQRTP) scholarship while conducting this research.