Truncated jarid2 and kdm6b transcripts are associated with temperature-induced sex reversal during development in a dragon lizard.

Sex determination and differentiation in reptiles are complex. In the model species, Pogona vitticeps, high incubation temperature can cause male to female sex reversal. To elucidate the epigenetic mechanisms of thermolabile sex, we used an unbiased genome-wide assessment of intron retention during sex reversal. The previously implicated chromatin modifiers (jarid2 and kdm6b) were two of three genes to display sex reversal-specific intron retention. In these species, embryonic intron retention resulting in C-terminally truncated jarid2 and kdm6b isoforms consistently occurs at low temperatures. High-temperature sex reversal is uniquely characterized by a high prevalence of N-terminally truncated isoforms of jarid2 and kdm6b, which are not present at low temperatures, or in two other reptiles with temperature-dependent sex determination. This work verifies that chromatin-modifying genes are involved in highly conserved temperature responses and can also be transcribed into isoforms with new sex-determining roles.

INTRODUCTION

Sex in vertebrates is determined by various factors, ranging from genes on sex chromosomes [genetic sex determination (GSD)], environmental factors such as temperature [temperature-dependent sex determination (TSD)], and interactions between genes and the environment (–). While the GSD pathways of mammals are well understood, the more complex systems involving both genes and environmental cues of other vertebrate lineages are still poorly characterized. Since the discovery of TSD in the lizard, Agama agama, over 50 years ago (), 15% of vertebrate species, particularly fish and reptiles, have been shown to have TSD or other forms of environmental sex determination (). Although many hypotheses have been proposed, precisely how an environmental cue is received and transmitted to the regulatory and epigenetic processes that determine sex remains a mystery. It is not also clear how, once the environmental signal is captured by the cell, the signal is transduced to influence the regulatory processes that determine sexual fate.

The dragon lizard Pogona vitticeps is a model organism well placed to explore these questions. It has a ZZ/ZW system of chromosomal sex determination (), but this system is subject to sex reversal of the ZZ genotype to a female phenotype by high temperatures (above 32°C), both in the laboratory () and in the wild (, ). Sex-reversed ZZf females are viable and fertile. Lysine demethylase 6B (KDM6B) acts on the histone H3K27me3 (tri-methylation of lysine 27 on histone H3) (). A second gene product JARID2 (encoded by Jumonji and AT-rich interaction domain containing 2, jarid2) is a component of PRC2 (polycomb-repressive complex 2), a complex responsible for the methylation of histone H3 to yield H3K27me3 marks (, ) and so is thought to act in opposition to KDM6B, although this action has not yet been demonstrated in reptiles. The chromatin dynamics that genes kdm6b and jarid2 regulate are essential for proper embryonic development through the establishment of cell fate and ensuring the transmission of the chromatin state through cell divisions (–).

Evidence is gathering of a role for alternative splicing as an important regulatory process during embryonic development, producing different transcript isoforms that can encode divergent protein variants in embryos and adults, and at different stages of embryonic development (–). Temperature-dependent alternative splicing is commonly observed in TSD reptiles and has been suggested to be an essential process in both TSD (–) and sex reversal (, ). For example, genes kdm6b and jarid2 and their differential splicing patterns have been implicated in the regulation of TSD in two species, the American alligator (Alligator mississippiensis) and the red-eared slider turtle (Trachemys scripta) (–). Functional work in T. scripta has shown that KDM6B demethylates the promotor of a key sex gene dmrt1 () to determine the trajectory of the sex differentiation pathway (, ). Intron retention in kdm6b and jarid2 also occurs in gonads of adult central bearded dragons (P. vitticeps) that were sex-reversed by high incubation temperature (). Homology of the kdm6b and jarid2 elements, which display intron retention, indicates an attribute that has remained conserved over approximately 200 million years of divergence ().

Our previous research suggested that intron retention in kdm6b and jarid2 during sex reversal alters histone demethylation actions and modifies the dosage of unknown sex-determining genes, thus redirecting developmental canalization that otherwise would have followed the homogametic sex (ZZ, male) (). Analysis of embryonic transcriptomes in P. vitticeps showed significant up-regulation of kdm6b and jarid2 in sex-reversing embryos when compared with embryos at a normal incubation temperature (). The next step is to examine splicing patterns during embryonic development and specifically during sex reversal in P. vitticeps.

Here, we present the first global analysis of RNA splicing during sex reversal in gonadal tissues of embryonic P. vitticeps. We show that previously unknown splicing events occur exclusively during sex reversal, suggesting a role for alternative splicing in responding to high temperatures by initiating and maintaining the sex reversal cascade in P. vitticeps. Contrary to the expression patterns we previously observed in adult P. vitticeps (), intron retention in jarid2 and kdm6b occurs at normal incubation temperatures during embryogenesis, whereas perfect splicing of these genes occurs only at sex-reversing temperatures. We also show that two N-terminally truncated isoform variants of jarid2 and kdm6b transcripts are present uniquely in sex-reversed embryos. The splicing profile of sex-reversing P. vitticeps embryos also differs from that of two TSD species, A. mississippiensis and T. scripta, suggesting that splicing regulation of jarid2 and kdm6b during sex reversal is distinct from typical TSD pathways. This research highlights an important role for specific splicing events during sex reversal that may occur in species with thermolabile sex more broadly and suggests that whereas ancient epigenetic modifier genes may be commonly implicated in thermosensitive sex determination pathways, the means by which they are regulated may not be conserved.

RESULTS

Global analysis of intron retention during embryonic development

An unbiased analysis of intron retention events across the whole gonadal transcriptome using IRFinder () identified only a few genes with differentially retained introns (data S1). Aside from kdm6b and jarid2 (discussed in detail below), only one gene, fbrs (fibrosin), exhibited differentially retained introns between sex-reversed ZZf females and embryos of both sexes (ZWf and ZZm) incubated at 28°C at every developmental stage. Transcripts fbrs retain introns 16 and 17 more frequently in ZZf females, while intron 13 was retained more frequently in both ZWf and ZZm embryos at 28°C compared with 36°C sex-reversed ZZf (fig. S1). All retained introns contain numerous stop codons in every phase, so intron-retaining transcripts are expected to be nonfunctional. Fbrs was up-regulated during sex reversal (), and its association with the noncanonical PRC1 complex may indicate that it has a role in chromatin remodeling in the sex reversal cascade (). Intron retention may therefore fine-tune the availability of this regulator. Beyond this link, fbrs has no known function in sex determination or differentiation, so the role of these isoforms in P. vitticeps development remains unknown and is not explored further here.

The generality of intron retention in reptile sex determination was examined by comparing transcripts in alligator and turtle to those in P. vitticeps, and particular note was taken of genes that showed differential intron retention across reptile clades. IRFinder analysis of a published RNA sequencing (RNA-seq) dataset of an embryonic time series of A. mississippiensis () and T. scripta () showed that jarid2 and kdm6b were the only genes with differentially retained introns that were shared between all three species, consistent with our previous results ().

One other gene, mov10, showed an infrared event in both P. vitticeps and T. scripta but not alligator. This gene [Mov10 RISC complex RNA helicase] encodes an RNA helicase with diverse roles in mRNA export and nonsense-mediated decay and has no known roles in sex determination (). Two other genes show intron retention events in both A. mississippiensis and T. scripta, but not P. vitticeps: sra1 (steroid receptor RNA activator 1) and rsrp1 (arginine- and serine-rich protein 1). Sra1 functions as an RNA binding protein and transcriptional regulator, particularly of hormone-related genes, but has not previously been implicated in sex determination (). Rsrp1 (paralogue of srsf6) encodes a splicing factor, the phosphorylation of which is controlled by CDC-like kinase 4 (CLK4) (). In T. scripta, six other genes displayed intron retention events that were not observed in the other species (data S1). We then focused on the two chromatin remodeling genes jarid2 and kdm6b, as they were the only genes that showed differential intron retention in all three species.

Splicing of chromatin remodeling gene jarid2

In P. vitticeps, retention of intron 15 in jarid2 (Fig. 1A) and intron 18 in kdm6b (Fig. 2A) occurred only in normal ZWf and ZZm embryos incubated at 28°C. This is the opposite pattern to the one previously described for adult P. vitticeps, in which intron retention occurred in sex-reversed ZZf females, while retention was rare in ZWf females and ZZm males ().

Fig. 1.

Isoforms of jarid2 present during embryonic development in P. vitticeps at normal (28°C) and sex reversal–inducing (36°C) temperatures.

(A) Read depth of stage 15 sex-reversed female (red), concordant female (light blue), and male (dark blue) depicting the intron retention event in intron 15 in jarid2. (B) The rate of intron 15 retention across the three developmental stages assessed for the three sex classes (note that it is zero for all ZZf at 36°C). (C) Normalized read depth for jarid2 depicting the marked increase in read depth in sex-reversed females (red) likely indicating a 5′ truncated isoform. (D) Average read depth of upstream (colored fill) and downstream (no fill) regions of jarid2 showing a significant increase in read depth in sex-reversed females (red).

Fig. 2.

Isoforms of jarid2 present during embryonic development in P. vitticeps.

(A) Transcript structure of the three observed jarid2 isoforms and the potential resulting protein variants with important functional domains indicated. (B) Hypothesized function of JARID2 isoforms during sex determination in P. vitticeps. Because of truncation of a PRC2-binding and -stimulating domain, the ability of ΔN-JARID2 to bind PRC2 and stimulate might be negatively affected. It is likely still able to bind typical target regions and blocks PRC2 from methylating these sites. The C-terminally truncated isoform produced from the retention of intron 15 may be detained in the nucleus (), or it may still be translated into a protein (ΔC-JARID2) but with reduced DNA binding affinity or altered target specificity, owing to the missing C-terminal zinc finger domain. (C) Relative abundance of the three jarid2 isoforms observed at stage 6 in the three sex classes. UIM, ubiquitin interaction motif; ZF, zinc finger; lncRNA, long noncoding RNA.

For jarid2, premature stop codons in intron 15 would lead to premature termination of translation, producing truncated protein isoforms lacking the C-terminal zinc finger domain (ΔC-JARID2). Thus intron retention occurs in normal ZZm males and ZWf females in the embryos but in sex-reversed ZZf female adults in P. vitticeps.

We investigated this unexpected finding by studying jarid2 intron retention in detail during embryogenesis in P. vitticeps. At 28°C, approximately 40% of jarid2 transcripts retained intron 15 in both ZWf females and ZZm males (Fig. 1B), compared with nearly zero intron retention at 36°C across all developmental stages assessed (Fig. 1B).

We also found another jarid2 isoform by examining read depth throughout the gene. We found that read depth of jarid2 in sex-reversed ZZf embryos markedly increased from slightly downstream from the start of exon 7 to the end of the transcript (Fig. 1C). This strongly indicated the presence of a jarid2 transcript isoform lacking the 5′ region in ZZf embryos, which was absent in both normal ZWf females and ZZm males incubated at 28°C. IRFinder analysis detected partial retention of intron 6, the intron preceding exon 7, in ZZf embryos (data S1). More than 90% of isoforms lacked the 5′ part of the transcript in ZZf embryos (Fig. 1D). The read coverage over intron 6 in ZZf embryos is unlikely to result from an intron retention event, so an alternative transcription start site within intron 6 or at the beginning of exon 7 is a more likely explanation (Fig. 3A). This interpretation is supported by uneven read depth over intron 6, which increases from 5′ toward 3′, and the higher number of reads spanning the intron 6–exon 7 boundary compared to the exon 6–intron 6 boundary. The most significant read depth increase occurred at the beginning of exon 7, suggesting that most transcripts start at the beginning of exon 7. The first two potential AUG start codons after the strongest incline in read depth are in frame with the canonical protein and display strong protein translation initiation sites (Kozak sequence AAGAUGAGA and AAAAUGGA).

Fig. 3.

Isoforms of kdm6b present during embryonic development in P. vitticeps at normal (28°C) and sex reversal–inducing (36°C) temperatures.

(A) Read depth of stage 15 sex-reversed female (red), concordant female (light blue), and male (dark blue) depicting the intron retention event in intron 18 in kdm6b. (B) The rate of intron 18 retention across the three developmental stages assessed for the three sex classes (note that it is zero for all ZZf at 36°C). (C) Normalized read depth for kdm6b depicting the marked increase in read depth in sex-reversed females (red) that corresponds with the ΔN-KDM6B isoform. (D) Average read depth of upstream (colored fill) and downstream (no fill) regions of kdm6b showing a significant increase in read depth in sex-reversed females (red).

Translation of the 5′ truncated transcript is likely to be initiated at one of those two initiation codons. The resulting N-terminal truncated protein (ΔN-JARID2) lacks the PRC2-stimulating domain and the domain responsible for binding to SUZ12, a core subunit of the PRC2 complex (Fig. 3C). A similar variant has been recently described in human epidermal keratinocytes ().

Splicing of chromatin remodeling gene kdm6b

In sex-reversed ZWf embryos, approximately half of the kdm6b transcripts retain intron 18 (Fig. 2B). This would be translated into a truncated protein that retains the demethylating functions of the JmjC domain (ΔC-KDM6B) but lacks the C terminus containing zinc finger and C-helix domains (Fig. 4A). The role of the ΔC-KDM6B variant and its potentially altered target specificity during sex reversal requires further investigation.

Fig. 4.

Isoforms of kdm6b present during embryonic development in P. vitticeps.

(A) Transcript structure of the three observed kdm6b isoforms and the potential corresponding protein variants with functional domains indicated. (B) Hypothesized function of KDM6B isoforms during sex determination in P. vitticeps. The N-terminal region is poorly described for KDM6B, so we hypothesize that although the annotated cNLS are truncated, given the high expression of kdm6b, it is likely that ΔN-KDM6B is still translated. We also suggest that it retains its demethylation function owing to retention of the JmjC domain, allowing it to demethylate target genes required for the initiation of sex reversal. The C-terminally truncated isoform produced by retention of intron 18 may be detained in the nucleus and not translated (), or it may still be translated (ΔC-KDM6B) but have reduced DNA binding affinity or altered target specificity owing to the missing zinc finger and C-helix C-terminal domains. (C) Relative abundance of the three kdm6b isoforms observed at stage 6 in the three sex classes.

We found that, as for jarid2, read depth for ZZf embryos in kdm6b increased considerably from exon 10 to the end of the transcript, indicating the presence of an isoform lacking the 5′ end (Fig. 2D). Intron 9, preceding exon 10, was also detected in the IRFinder analysis as partially retained in ZZf embryos (data S1). An alternative transcription start site within intron 9 of kdm6b seems a likely explanation. In such a 5′ truncated transcript, the first AUG triplet occurs in exon 14, which contains four AUG triplets in frame with the canonical KDM6B protein. The N-terminal truncated protein translated from any of the four AUG codons would start just before the demethylating JmjC domain. Previous research characterizing KDM6B protein function in mammalian cell lines revealed the presence of classical nuclear localization signals (cNLSs) in the N-terminal region that are required for nuclear localization and subsequent demethylation actions of KDM6B (). Most kdm6b transcripts in sex-reversed ZZf embryos lack the 5′ end, so the putative ΔN-KDM6B protein is likely to have altered localization dynamics, but this remains to be tested (Fig. 4A).

ΔN isoforms unique to sex reversal in P. vitticeps

Investigation of jarid2 and kdm6b read depth in T. scripta and A. mississippiensis showed no strong indication of 5′ truncated transcripts, as were seen in P. vitticeps (figs. S2 and S3). The slight increase in read depth of downstream exons compared to upstream exons is likely due to a typical 3′ bias in polyadenylate-selected RNA-seq data. In both T. scripta and A. mississippiensis, there are alternative translation start sites (AUG) in jarid2 and kdm6b; therefore, the production of N-terminally truncated proteins may be possible, but not prominent, at the incubation temperatures examined. Thus, both chromatin remodeling genes, jarid2 and kdm6b, show differential intron retention and alternative transcription start sites at different temperatures in P. vitticeps. Truncation of the 5′ region of the transcripts leads to proteins lacking the functional domains in these regions. In ΔN-JARID, loss of the PRC2-interacting domains could lead to competitive inhibition of the canonical JARID, whereas in ΔN-KDM6B, the JmjC domain is retained, so it probably still exhibits lysine demethylase activity.

DISCUSSION

In this study, we demonstrate temperature-specific alternative splicing events and alternative initiation sites in jarid2 and kdm6b, which all are likely to play key regulatory roles in sex reversal in P. vitticeps. We show that intron retention in kdm6b and jarid2 occurs only at low incubation temperatures during embryogenesis and changes to the opposite pattern in adults ().

As previously reported, the same intron retention event causing C-terminally truncated kdm6b and jarid2 isoforms occurs exclusively during low temperature incubation and not high temperature incubation of A. mississippiensis, T. scripta, and P. vitticeps (). Thus, our present observations on P. vitticeps embryos now demonstrate that the pattern of intron retention in embryos is the same in all three species and occurs at low temperatures in the three widely divergent reptiles. These low temperature intron retention events are distinct from the 5′ truncated kdm6b and jarid2 isoforms, which occur only during high-temperature sex reversal in P. vitticeps embryos, and are absent in T. scripta and A. mississippiensis embryos under the incubation temperatures investigated here.

These new findings have important implications for our understanding of the thermolabile sex determination systems in these three species. For P. vitticeps, we now know that the patterns of isoform expression are not consistent between embryos and adults, as was previously assumed (). It seems unlikely that the ΔC-JARID2 and ΔC-KDM6B variants have a role in sex determination, as they were observed in the GSD pathways of both ZZm males and ZWf females at 28°C (Fig. 5B). The dynamics of protein translation and function of both the ΔC and ΔN variants present exciting direction for ongoing research.

Fig. 5.

Hypothesized pathway for the role of the ΔN jarid2 and kdm6b transcript and protein isoforms during temperature-induced sex reversal in P. vitticeps embryos and the presence of isoforms at different incubation temperatures in P. vitticeps, A. mississippiensis, and T. scripta.

(A) We hypothesize that under high incubation temperatures, splicing of jarid2 and kdm6b is altered, potentially by CLK4-modulating SR proteins, to create the ΔN variants. If translated to a protein, then the ΔN-JARID2 isoform lacking the PRC2-interacting domains can compete with the canonical JARID2 protein and the PRC2 complex by binding at the same target sites. The ΔN-KDM6B isoform likely retains its demethylating functions at the JmjC that is intact. Together, the ΔN variants may decrease levels of H3K27 methylation (a repressive chromatin mark), causing an increase in female sex gene expression required for the initiation and maintenance of male to female sex reversal. This hypothesized pathway can be followed from left to right in the diagram. (B) Schematic representation of the presence of ΔN and ΔC jarid2 and kdm6b isoforms in the three species at the incubation temperatures that have been assessed. For A. mississippiensis and T. scripta, the dashed lines at high incubation temperatures denote the possibility that the ΔN variants may be present under these conditions but have not been assessed. The color gradients denote where the variants have been demonstrated at a given incubation temperature, and the question marks in A. mississippiensis and T. scripta denote high incubation temperatures where the ΔN variants may occur but have not been assessed.

It would be interesting to find whether more extreme incubation temperatures in T. scripta and A. mississippiensis induce expression of the ΔN-JARID2 and KDM6B isoforms (Fig. 5B). In T. scripta, it has been conclusively demonstrated that at the male-producing temperature (MPT), kdm6b is highly expressed and is required to demethylate the dmrt1 promoter and initiate male development (). We also know that the ΔC-KDM6B variant is expressed at the MPT in this species, but the functional implications of this for the induction of the male sex determination cascade remains unknown, as does the role of the ΔC-JARID2 isoform.

Much less is known about the sex determination pathways in A. mississippiensis. Of particular interest is that in this species, females can be produced at high and low incubation temperatures (approximately 35°C) (). However, no transcriptomic work has yet been conducted at such a high incubation temperature. It is possible that this more extreme female-producing temperature may generate a similar cellular effect as is seen during sex reversal in P. vitticeps, and the ΔN-JARID2 and KDM6B variants may be expressed (Fig. 5B). If this is the case, then it would have important implications for the evolution of the ΔN variants in thermosensitive sex determination systems. On the basis of previous research (), we propose a crucial role for the splicing regulator CLK4 in the P. vitticeps sex reversal cascade (Fig. 5A). This research assessed the activity profiles of CLK4 in different species and found that it functioned at temperatures relevant to each species (). The activity of CLK4, as well as the serine- and arginine-rich (SR) proteins that it regulates, in P. vitticeps during sex reversal remains unknown.

In P. vitticeps, the alternative jarid2 transcript is expected to translate into ΔN-JARID2 protein in ZZf females. This is a new variant of JARID2 (ΔN-JARID2) present in gonadal cells from stage 6. A similar truncated ΔN-JARID2 transcript was characterized recently in human epidermal keratinocytes () and was suggested to display competitive inhibition of PRC2 on target sites. Whereas full-length JARID2 recruits PRC2 to its target sites by binding both nucleosome and PRC2 domains and stimulates the H3K27-methylating activity of PRC2, ΔN-JARID2 was suggested to still bind to its target sites but to block PRC2 binding and/or activity at those sites because it lacks the PRC2-binding and -stimulating domains (Fig. 3B) (). Thus, in mammalian cells at least, the function of ΔN-JARID2 appears to be antagonistic to the canonical JARID2 isoform. We suggest that this mode of action may be critical during sex reversal to alter the expression of sex specific genes through modulation of H3K27 methylation, initiating male to female sex reversal from stage 6 of development (Fig. 5A).

Another interesting parallel between human and dragon ΔN-JARID2 variants is their association with cellular calcium concentration. In human cell lines, ΔN-JARID2 is observed when cellular calcium concentrations are increased (). The same may be true in dragons, as we know from previous work that activation of the calcium channel (trpv2) and a suite of genes associated with calcium signaling are up-regulated during high-temperature sex reversal (). Whether the ΔN-JARID2 isoform is present in other species, or whether it can be generated under heat stress conditions, and what its role may be, are currently unknown.

Human and dragon ΔN-JARID2 forms appear to be similar and their function may be analogous, but the variants are generated via different mechanisms, constituting an unusual form of convergent evolution. Human ΔN-JARID2 is generated via posttranslational cleavage of the protein, whereas dragon ΔN-JARID2 is generated via the activation of an alternative transcription start site.

We also identified a new 5′ truncated kdm6b transcript isoform, which could be translated into a truncated ΔN-KDM6B protein. Newly identified in this study, the function of the ΔN-KDM6B variant is unknown. We hypothesize that, as has been shown for ΔN-JARID2, it has an antagonistic function to the canonical JARID2, maintaining the active demethylated state of key genes necessary to initiate and promote the sex reversal pathway, diverting development away from the normal male development of the homogametic sex (Fig. 5). In mammalian cell lines, cNLSs in the N-terminal region regulate the nuclear localization of KDM6B (). As balanced KDM6B localization between the nucleus and cytoplasm is a mechanism by which H3K27me3 methylation levels can be regulated, we propose that modification of this histone is involved in P. vitticeps sex reversal, providing a compelling avenue for future research. KDM6B-mediated demethylation of male pathway gene dmrt1 has been demonstrated in T. scripta; however, note that it is not known if the canonical isoform, ΔC or both KDM6B isoforms, is responsible for dmrt1 demethylation (, ). We propose that, in the P. vitticeps male to female sex reversal system, ΔN-KDM6B demethylates female pathway genes, as this is what would be required to initiate the female pathway (Fig. 5A). Future experiments assessing the function of the ΔN-KDM6B variant in P. vitticeps are required and could include investigation of the target sites of this variant and changes in methylation dynamics of candidate genes such as cyp19a1 or foxl2.

Comparisons of transcript isoforms in P. vitticeps, A. mississippiensis, and T. scripta show that JARID2 and KDM6B are consistently implicated in thermosensitive sex determination cascades across evolutionarily disparate lineages. However, there are critical differences between the three species, highlighting that differing selective pressures have driven these genes to acquire differing functions via isoforms lacking DNA or protein-binding domains that change target specificity or ligand-binding efficiency. Chromatin remodeling genes such as KDM6B and JARID2, which can operate across many genomic sites, are ideal candidates for these types of dynamic evolutionary shifts. We show that the 5′ truncated transcript isoforms that encode the truncated ΔN-JARID2 and ΔN-KDM6B proteins appear only in sex-reversed embryonic P. vitticeps, suggesting a unique pathway for the evolution of these genes. Investigation of isoform variants in other species with thermosensitive sex determination is needed to determine whether these variants are specific to P. vitticeps or are associated with sex reversal or TSD systems more broadly.

Our observations of these transcript variants in P. vitticeps raises important questions: What regulates splicing and transcription start site choice for these genes, and how are they related to thermosensitivity and sexual development? There are likely to be many other genes involved in thermal sensitivity of sex reversal in P. vitticeps. CLK4, a protein kinase that regulates the activation of SR proteins that control mRNA splicing (), is a very promising candidate. CLK4 activity controls jarid2 intron retention in two turtle species, which may play a role in regulating the splicing patterns we observed during sex reversal in P. vitticeps (Fig. 5A). The findings from this research open up many new avenues for future research, including determining the range of temperatures at which the ΔN variants may be functioning in P. vitticeps and demonstrating their function in regulating sex determination cascades. The effect temperature fluctuations during incubation in natural nests on splicing modulation also needed to be investigated ().

While much remains to be understood about thermosensitive sex determination systems, our demonstration of unique splicing events during reversal in P. vitticeps represents an important advance. We reveal complex and evolutionary dynamic roles for jarid2 and kdm6b across divergent lineages, bringing new insight into the elusive molecular mechanisms by which temperature determines sex. More broadly, our work adds to the body of evidence demonstrating that jarid2 exists in different forms that can have completely opposing functions. Fully understanding the complexities of splicing in these chromatin modifying genes may redefine our understanding of the dynamics of transcriptional activation and repression during vertebrate development.

MATERIALS AND METHODS

Egg incubations, sampling, and sequencing for P. vitticeps

The dataset used in this study has been described in full in (). Briefly, eggs from ZZf and ZWf mothers were collected at or soon after lay and incubated at 28°C. Following a period of developmental entrainment (10 days), eggs from ZZf females were switched to 36°C. Eggs at each temperature were sampled at stages 6, 12, and 15 according to the staging criteria established for P. vitticeps () to capture the bipotential gonads and early and late stages of gonad differentiation (). The isolated embryonic gonads were dissected, and total RNA was extracted using the QIAGEN RNeasy Micro Kit (catalog no. 74004) according to the manufacturer’s protocols. Sequencing libraries were prepared in randomized batches using the Roche NimbleGen KAPA Stranded mRNA-seq Kit (catalog no. KK8420), and samples were sequenced on the Illumina HiSeq 2500 system at the Kinghorn Center for Clinical Genomics (Garvan Institute of Medical Research, Sydney). On average, 25 million read pairs per sample were obtained.

Data analysis

The paired-end RNA-seq libraries for P. vitticeps were trimmed using cutadapt () with -q 20 -m 20 -max-n 4 -trim-n. Trimmed reads were aligned to the P. vitticeps genome assembly pvi1.1 (GCA_900067755.1; http://ftp.ensembl.org/pub/release-100/fasta/pogona_vitticeps/dna/Pogona_vitticeps.pvi1.1.dna.toplevel.fa.gz) () using STAR (v2.7.0f) () with splice-aware alignment guided by the accompanying Ensembl gene annotation (pvi1.1.100). Parameters were chosen to output only unique alignments and to assure compatibility with IRFinder (--outFilterMultimapNmax 1 --outSAMstrandField intronMotif --outFileNamePrefix --outSAMtype BAM Unsorted). To reduce the number of scaffolds in the genome assembly, those with no genes annotated were excluded. Paired-end sequencing libraries for A. mississippiensis were obtained from () data accessioned on the DNA Databank of Japan Sequence Read Archive accession number DRA004128-41. Paired-end sequencing libraries for T. scripta were obtained from () data accessioned on National Center for Biotechnology Information (NCBI)’s Sequence Read Archive database (Bioproject PRJNA331105). Because of insufficient replicates for robust analysis, A. mississippiensis samples DDR048634-DDR04636 and T. scripta stage 26 samples were excluded. The libraries were trimmed and aligned as described above for P. vitticeps. For A. mississippiensis, the genome assembly ASM28112v4 (RefSeq GCF_000281125.3) was used with accompanying NCBI gene annotation (GCF_000281125.3_ASM28112v4_genomic.gtf). For T. scripta, the genome assembly CAS_Tse_1.0 (RefSeq GCF_013100865.1) was used with accompanying NCBI gene annotation (GCF_013100865.1_CAS_Tse_1.0_genomic.gtf). Approximations of the developmental stages for A. mississippiensis and T. scripta with the staging system from P. vitticeps are provided in fig. S4.

To conduct an unbiased analysis of splicing patterns across whole transcriptomes obtained for the three species, IRFinder () was used to detect differentially retained introns between ZZf and ZWf embryos across development in P. vitticeps and between male and female producing temperatures across development in A. mississippiensis and T. scripta. For T. scripta, two jarid2 transcript isoforms were removed from the annotation gtf file (XM_034762362.1 and XM_034762363.1). They were not expressed in any sample but interfered with correct recognition of intron retention events in the major transcript isoform. Before differential intron retention analysis, introns were excluded with an intron depth (intronDepth) or a maximum of spliced reads (maxSplice) of ≥5 in less than two samples. For P. vitticeps and A. mississippiensis, the differential intron retention analysis was performed with DESeq2 (1.26.0) () in R version R 3.6.1 as outlined in the IRFinder manual. For T. scripta, an Audic and Claverie test with pooled samples was performed according to IRFinder guidelines for samples with two or less replicates. For differentially retained introns, the following cutoffs were applied: intron retention of ≥0.1 in either group, intron retention change of ≥0.25, adjusted P < 0.01 for P. vitticeps, and P < 0.01 for A. mississippiensis and T. scripta. Bedgraphs were created for each alignment file with bedtools genomecov (https://bedtools.readthedocs.io/en/latest/index.html) and used to extract read depth over jarid2 (P. vitticeps, ENSPVIT00000008211; T. scripta, XM_034762361.1; A. mississippiensis, XM_006270617.3) and kdm6b (P. vitticeps, ENSPVIT00000021647; T. scripta, XM_034757088.1; A. mississippiensis, XM_019493797.1) transcripts with costume scripts. To determine the relative abundance of the 5′ truncated jarid2 and kdm6b transcript isoforms in P. vitticeps, the read depth of a region upstream and downstream of the putative alternative transcription start site was compared. For P. vitticeps, the kdm6b upstream and downstream regions comprise 293 nucleotides (nt) (covering exons 8 and 9) and 898 nt (covering exons 14 to 16), respectively. The jarid2 upstream and downstream regions comprise 721 nt (covering exons 3 to 6) and 799 nt (covering exons 8 to 11), respectively. Read depth was normalized by library size and length of the region. The average ± SD of all samples per group is shown.

For similar analysis in T. scripta and A. mississippiensis, upstream and downstream regions were chosen according to the putative ∆N protein isoforms of P. vitticeps, that is, for jarid2 regions upstream and downstream of the JmjN domain and for kdm6b regions far upstream and downstream of the JmjC domain. For T. scripta, the kdm6b upstream and downstream regions comprise 65 nt (covering exon 2) and 149 nt (covering exon 4), respectively. There is an assembly gap in the T. scripta genome assembly GCF_013100865.1 within the beginning of kdm6b. Therefore, 1587 nt in the beginning of kdm6b gene annotation in GCF_013100865.1_CAS_Tse_1.0_genomic.gtf is missing. The jarid2 upstream and downstream regions comprise 407 nt (covering exons 1 and 2) and 587 nt (covering exons 4 and 5), respectively. For A. mississippiensis, the kdm6b upstream and downstream regions comprise 396 nt (covering exons 4 to 6) and 452 nt (covering exons 15 to 17), respectively. The jarid2 upstream and downstream regions comprise 716 nt (covering exons 3 to 6) and 798 nt (covering exons 8 to 11), respectively. Read depth was normalized by library size and length of the region. The mean ± SD of all samples per group is shown.

The relative abundance of all three jarid2 and kdm6b transcript isoforms (full-length, retained intron, and 5′ truncated) in ZZf, ZWf, and ZZm in P. vitticeps was inferred from integration of gene differential expression analysis (), differential intron retention analysis, and the abundance calculation for the 5′ truncated isoform as explained above.

We follow the standard set for the naming of genes (e.g., kdm6b for reptiles and KDM6B mammals) by Kusumi et al. (). Proteins are uppercase.