Multi-year pair-bonding in Murray cod (Maccullochella peelii)

Identifying pair-bonds in an Australian Freshwater Fish

The iconic Australian Murray cod (Maccullochella peelii) (Mitchell, 1938) is one of four morphologically cryptic but genetically distinct species within Maccullochella  (Rowland, 1993; Nock et al., 2010). The Murray cod is Australia’s largest freshwater fish and can grow to as large as 180 cm in length with a maximum recorded weight of 113.6 kg. Its fecundity ranges from 9,000–120,000 eggs annually. The species is highly sought after by anglers. It is limited to parts of the Murray-Darling Basin (MDB) and is listed as vulnerable under the Environment Protection and Biodiversity Conservation Act (1999) (Australian Department of Agriculture, Water and the Environment, 2020). Murray cod still form an important recreational fishery (Lintermans & Phillips, 2005; Koehn & Todd, 2012; Ye et al., 2016).

Murray cod is a species worth studying for potential evidence of long-term pair-bonding and alternative mating strategies. Male Murray cod are known to provide parental care to both eggs and larvae for up to 20 days after preparing a nest area. They are long-lived (up to 48 years), slow-maturing (Lintermans, 2007) and show high site fidelity (Koehn & Nicol, 2016). These are life history factors associated with long-term pair-bonding in other species Kleiman, 1981; Barlow, 1988; Arnold & Owens, 1998; Hatchwell & Komdeur, 2000; Chapple, 2003. Murray cod have been shown to exhibit monogamous mating, as well as polygyny and polyandry over three breeding seasons when held captive in ponds (Rourke et al., 2009).

In this study single nucleotide polymorphisms (SNPS) were used to look for sibling relationships of Murray cod larvae and these data, in combination with relevant metadata, were used to infer the existence of multi-year pair-bonding as a mating strategy. The probability of multi-year pair-bonding occurring is considered. Although SNPS are increasingly being used to identify parentage (Vandeputte & Haffray, 2014; Huisman, 2017) confidence in the methods is still developing. For this reason stable isotope evidence is provided to support the genetic findings. This reduces the possibility that any relatedness found was not due to contamination during sample processing, or is merely an artefact of the probabilistic approach to inferring sibship relationships in the absence of known parentals.

Study site

We sampled Murray cod larvae from six sites along a 50 km upland reach of the in Murrumbidgee River in the Australian Capital Territory (ACT), Australia (Fig. 1).

Larval fish collection, identification, preparation and aging

Data were collected as previously described in Couch et al. (2016). Specifically we examined 261 of 2,607 Maccullochella larvae which were collected in 2011, 2012, and 2013 from the six sites (Fig. 1) using larval driftnets. The sympatric Maccullochella species Trout cod (M. macquariensis) and hybrid larvae were identified using a combination of SNPs and mitochondrial DNA sequencing (Couch et al., 2016). These were excluded from the dataset leaving 251 larval Murray cod used for this analysis. Larvae and tissue samples were preserved in 95% ethanol at room temperature until 2014 after which samples were stored at −20 °C. Fish were collected under ACT Government licences LT2011516, LT2012590 and LT20133653. The research was conducted under approvals CEAE 11-15 and CEAE 13-17 from the University of Canberra Committee for Ethics in Animal Experimentation.

Putative nest location

Individual nest sites are unknown so an estimate of the putative nest site was made based on a mean larval dispersal velocity of 700 m per day for the duration of the time available for dispersal based on the larva age (Couch, 2018).

Genomic DNA extraction and sequencing

Genomic DNA Extraction and Sequencing was performed as previously described in Couch et al. (2016). Total DNA was isolated from whole larval heads. The DNA extraction protocol is detailed in Couch & Young (2016) and is based on a turtle DNA extraction protocol (FitzSimmons, Moritz & Moore, 1995). Sequencing was done using Diversity Arrays Technology’s (DArT) DArTseq™ which represents a combination of DArT complexity reduction methods and next-generation sequencing platforms (Kilian et al., 2012; Courtois et al., 2013; Cruz, Kilian & Dierig, 2013; Raman et al., 2014). Sequences generated were processed using proprietary DArT analytical pipelines (http://www.diversityarrays.com).

Marker scoring and statistical analysis

Marker Scoring and Statistical Analysis was performed as previously described in Couch et al. (2016). Specifically, DArTsoft (Diversity Arrays Technology, Building 3, University of Canberra, Australia), a software package developed by DArT PL (http://www.diversityarrays.com/software.html), was used to both identify and score the markers that were polymorphic.

Carbon and nitrogen stable isotope analysis

Dried muscle material from each fish larva (0.86 mg ± 0.17 SD), the bulk of the posterior portion of its body without head and gut of the fish, were encapsulated in tin. Samples were combusted in an elemental analyser mass spectrometer (Sercon, Crewe, United Kingdom) at the Australian National University Research School of Earth Sciences Radiocarbon Laboratory, on a fee-for-service basis and assayed for δ¹⁵N and δ¹³C stable isotope ratios and C:N ratio. Isotopic signatures were determined based on Australian National University isotopic standards (USGU41, USGU40, Caffeine and Gelatine). Measurement precision was approximately 0.08 ‰for 13C and 0.15 ‰for 15N. Isotope values are expressed as the relative parts per thousand (‰).

Inferring existence of family groups using spatio-temporal data

Plots were made using hatch date and putative nest location for larvae from each annual cohort. Individual clusters, corresponding to putative nests, were identified and nominated.

Inferring sibship using related

The probability of relatedness (r) was calculated for individual larval dyads based on the trioml algorithm using ‘related’ (Pew et al., 2015).

A simulation was run to identify a probability density function for full siblings, parent–offspring, half-siblings, and unrelated individuals using the allele frequencies from the DARTseqs for each larva. A probability estimate that would best estimate the cut-off probability between full siblings and half-siblings was identified and used to subset likely full-sibling dyads from the larval dyads. Full siblings were assigned a name representing their putative ‘mother’. This was arbitrary and could have been assigned a name representing putative ‘father’.

The relationships between individuals based on the dyads and the probabilities was visualised by making network graphs and plotting them using Gephi (Bastian, Heymann & Jacomy, 2009) and ‘r’ package ‘iGraph’(Csardi & Nepusz, 2006). The set of dyads containing full siblings was used to prepare network graphs to facilitate visualisation of the family groups and the assignation of a name to putative female parent.

Inbreeding coefficient

The two likelihood algorithms - dyadML and trioml - as well as the lynchrd, and ritland algorithms within the ‘r’ package ‘related’ can account for inbreeding in their estimates of relatedness. The command “allow.inbreeding=TRUE” was set to output an inbreeding coefficient for each individual under each of the three algorithms above.

Comparison of inferred family groups and sibships

The plot illustrating family groups previously identified and named using spatio-temporal data alone was then coloured by the name of the putative mother. This comparison increased precision in identification of ‘family’ groups.

Accounting for outbred population

Using the measures of relatedness (r) calculated using the trioml algorithm it became apparent that there were two distinct groups within the larvae sampled; those that were highly outbred (r<−0.4), and the rest which were not strongly inbred or outbred (r>−0.4 and <0.3). Relatedness amongst the non-outbred fish was used to determine common parentals. Principal component analysis of those larvae considered to be outbred suggested a difference between the two populations. These differences were not correlated with location or year. Fish with a coefficient of inbreeding below −0.4 were considered outbreds and separate from the ‘river’ fish (Figure in supplementary material). These fish were excluded from subsequent related analysis as they were considered likely to be Murray cod introduced from a recent re-stocking program.

Inferring existence of family groups using spatio-temporal data alone

Initially, the existence of family groups was inferred from a scatterplot of spatial and temporal separation of larvae. Figure 2 shows clusters of larvae distributed over space and time.

Nomination of putative female parent

The set of dyads containing full siblings was used to prepare network graphs to facilitate visualisation of the family groups and the assignation of a name to putative female parent. Four female Murray cod mated with the same male for more than one year of the three years sampled (Fig. 3). One pair mated for three years sequentially and three other pairs mated for two of the three years studied, one pair detected mating sequentially and two pairs detected mating in non-sequential years.

Comparison of inferred family groups and sibships

A spatio-temporal plot was coloured by the putative female parent of the larvae (Fig. 4) and it is a clear correlation between identification of groups by spatio-temporal factors alone and those identified by genetic relatedness of the nest groups.

When the Carbon and Nitrogen Isotope ratios of the larvae are considered the strong clustering of CN isotope ratios and female parent illustrates a high correlation between full sibling status and body isotope signature provided to the larva from its female parent (Fig. 5). This also supports the validity of parentage assignment based on related analysis.

The identification of sibling relationship is also supported by the limited number of observations (just 3 individuals of 35 full sibling pairs) that were found at separate sites, and those three were only at detected at sites immediately adjacent to the other sibling. These larvae sibling pairs are 140,179; 133,105 and 170,191 can be seen in Fig. 4.

Inferring sibship using related

Simulation using the allele frequencies present in the population produced the probability density function in Fig. 6. In this way, the overlap between relatedness values can be assessed. In this case, a ‘cutoff’ value of any r value above 0.4 was selected to identify most full sibling dyads while minimising the possibility of inadvertently misclassifying half-siblings as full siblings (Fig. 6). The possibility of parent–offspring relationships is obviated because all larvae in the analysis were collected within three years and the sexual maturity of Murray cod is greater than 4–5 years (Lintermans, 2007).

After identifying the optimal value of r to eliminate dyads least likely to represent full siblings, the set of dyads were filtered to include full siblings only. This resulted in 35 dyads (pairs of full siblings) that were assigned to a family group.

Probability of observed multi-year pair-bonding

The multi-year bonding identified in this case is unlikely to be due to chance alone. The probability of multi-year pair-bonding occurring within the cohorts by chance can be estimated if:

y = number years (2 or 3) with the same mate; and

n = number of pairs (138) including singletons (larvae with no siblings found).

and we assume 50/50 sex ratio in accordance with the findings of previous research on the species (Cadwallader, 1977; Koehn & O’Connor, 1990) then the probability of pairing is: (5)${\displaystyle p={\left(1∕n-1\right)}^y-1.}$

Applying Eq. (5) for an individual female, the probability of mating with the same male for two years is <0.008 and for each individual female the probability of mating with the same male for three years is <0.00005. Clearly, the probability of four individuals each choosing the same mate for multiple years reduces this probability even further.

To put it another way; for p to even approach non-significance at the 0.05 level, the number of available males would need to be as low as five. Thus it seems unlikely to be random mate selection.