DOI: 10.5061/dryad.dbrv15f32
title: Epigenetics underpins phenotypic plasticity of protandrous sex change in fish

description: 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.
description: 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.
description: Please refer to the readme ('README.txt') and metadata
 ('00.0_c3_raw_sequence_metadata.csv') files for
 details.