Phenotypic plasticity is an important driver of species resilience. Often mediated by epigenetic changes, phenotypic plasticity enables individual genotypes to express variable phenotypes in response to environmental change. Barramundi (Lates calcarifer) is a protandrous (male-first) sequential hermaphrodite that exhibit plasticity in length-at-sex change between geographic regions. This plasticity is likely to be mediated by changes in DNA methylation (DNAm), a well-studied epigenetic modification. However, region-specific relationships between length, sex and DNAm in sequential hermaphrodites were previously unreported. To investigate these relationships, here we compare DNAm in four conserved vertebrate sex-determining genes in male and female barramundi of differing lengths from three regions of northern Australia. Despite a strong association between increasing length and male-to-female sex change, the data reveal that DNAm becomes more sex-specific (rather than more female-specific) with length. Significant differences in DNAm between males and females of similar lengths suggest that female-specific DNAm arises rapidly during sex change, rather than gradually with growth. The findings also reveal that region-specific differences in length-at-sex change are accompanied by differences in DNAm, and were concurrent with variability in remotely sensed sea temperature and salinity. Together, these findings provide the first in situ evidence for epigenetically and environmentally mediated sex change in a protandrous hermaphrodite, and offer significant insight into the molecular and ecological processes governing the marked and unique plasticity of sex in fish.

Animal collection and sampling design

Gonadal tissue samples were collected as part of the Fisheries Queensland, Department of Agriculture and Fisheries (DAF) barramundi biological monitoring program throughout the 2015 and 2016 barramundi wild harvest fishery open season (February through October) and stored in RNAlater (Thermo Fisher Scientific). Sex was determined macroscopically, whereby male and female gonads were identified by their distinct colour and texture (Fisheries Queensland, 2020a). No fish undergoing testes-to-ovary transition were recognised in this sample set, as their identification requires histological examination. Total length (to the nearest 10 mm) was measured from the anterior-most part of the snout to the posterior tip of the caudal fin. Fish age was determined by sectioned otolith examination including increment counts, edge classification and periodicity, as well as timing of opaque zone transformation. Preliminary analysis revealed differences in the length at which sex change occurs (length-at-sex change, herein) within GoC samples, thus, the southern GoC sub-population was split into ‘southern GoC’ (at ~ 16°S to land, near Karumba) and ‘mid-northern GoC’ (at ~ 13°S to 14°30"S near Aurukun). Only the northern Qld east coast region (at ~ 16°S to 18°S, near Cairns and Mission Beach), referred to as ‘NQ east coast’, was included due to insufficient sample size south of this latitude. Stratified random sampling was used to select gonad samples for downstream DNAm analysis. For each of the three geographic regions, 1 - 10 gonad samples (based on sample frequency and therefore availability) from fish from each 10 cm length category were analysed. Samples did not require collection permits or animal ethics approval, as the fish were taken by recreational, charter and commercial operators as part of usual fishing practices. Samples collected were fisheries regulation capture-dependent and restricted by minimum and maximum length limits of 58 and 120 cm, respectively.

Genomic DNA extraction

Gonadal tissue was removed from storage in RNAlater (Thermo Fisher Scientific), washed once in PBS and dried with a KimWipe (Kimtech Science) before immediate immersion into DNA extraction buffer [100 nM Tris-HCl, 1.4 M NaCl, 20 nM EDTA, 2 % Cetyl trimethylammonium bromide (CTAB), 2 % polyvinylpyrrolidone (PVP)]. Genomic DNA (gDNA) was extracted following the CTAB protocol, including overnight digestion with Proteinase-K followed by the addition of a phenol:chloroform:isoamyl alcohol (25:24:1) step to assist with the removal of proteins. Quantity and purity of gDNA was assessed using an ND-1000 spectrophotometer (Nanodrop technology) based on absorbance at 260 nm and 260/280 nm ratio, and integrity was assessed by visualisation on a 0.8 % agarose gel with lambda DNA standards at 50, 20, 10 and 5 ng/µl and a 1 kb Plus DNA ladder (Thermo Fisher Scientific).

Bisulphite conversion of gDNA and amplicon specific PCR

To analyse DNAm levels of the region surrounding the start codon of cyp19a1a, esr1, dmrt1 and nr5a2, a targeted bisulphite amplicon sequencing (BSAS) approach, adapted from Masser, Berg, and Freeman (2013), was applied. Following the manufacturer’s instructions, 500 ng of extracted gDNA was subject to bisulphite treatment using EZ DNA Methylation‐Gold™ (Zymo Research). Gene-specific primers spanning the start codon of each gene were designed using MethPrimer with the Illumina adapter overhang nucleotide sequences added. PCR amplification of these regions was achieved using Platinum® Taq DNA Polymerase (Thermo Fisher Scientific) following the manufacturer’s instructions. Reaction conditions were as follows: 95 °C for 2 min followed by 40 rounds of 95 °C for 30 s, 57.5 °C for 35 s, 72 °C for 40 s, with a final extension of 72 °C for 10 min. PCR products were purified using Sera-Mag SpeedBeads (GE Healthcare) and quantified using QuantiFluor (Promega) fluorometric nucleic acid quantitation and measured on an EnSpire Multimode plate reader (PerkinElmer).

NGS library preparation and DNA methylation quantification

Dual indexed libraries were generated using a Nextera XT Index Kit following the manufacturer’s protocol (Illumina). Purified PCR products were indexed using limited cycle number PCR under the following reaction conditions: 95 °C for 5 min followed by 12 rounds of 95 °C for 30 s, 58 °C for 35 s, 72 °C for 40 s, with a final extension of 72 °C for 10 min. Indexed amplicons were purified and quantified as described above, normalised by molarity to 4 nM and pooled into a final library. The final library was and quality checked on an Aligent TapeStation (Aligent Technologies) and quantified on a Qubit 3.0 Fluorometer (Invitrogen) for subsequent molarity determination. Libraries were then diluted to 8 pM, spiked with 20 % PhiX to compensate for low base diversity in the amplicon sequencing library, and loaded onto a 600 cycle V3 reagent cartridge for sequencing on an Illumina MiSeq (Illumina). FASTQ files were imported into Geneious Version 10.2.6. Paired reads were merged using BB merge paired Read Merger Version 37.25, sequences were then trimmed based on an error probability limit of 0.05 with a maximum of one ambiguity and resultant reads were aligned to in-silico bisulphite converted reference sequences using the Geneious in-built read mapper end to end read mapping and a minimum mapping quality of 30. Variants were detected using a minimum coverage of 500 reads to identify C-T conversions in CpG positions. Because only unmethylated cytosines are converted to thymines through bisulphite treatment, at a given CpG site, the proportion of reads with C-T variants indicates the proportion of unmethylated gene copies in the original gDNA sample. The inverse of this value gives the proportion of cytosine methylation (DNAm) for each CpG site for each gene amplicon.

Please refer to the readme ('README.txt') and metadata ('00.0_c3_raw_sequence_metadata.csv') files for details.