Collection of eggs and adults

We collected eggs from snapping turtles nests in June 2008, 2009, and 2010. Adult turtles were collected in June, July, and August of the same years. Collection sites extended from the Iowa border to the Canadian border within the state of Minnesota. Adult snapping turtles were also collected on the Lake Lewisville U.S. Army Corps of Engineers property, Lewisville, Texas, in May of 2013 and 2014. Eggs and adults were transported to the animal quarters in the Biology Department at the University of North Dakota. Experiments were carried out according to a protocol approved by the Institutional Animal Care and Use Committee at the University of North Dakota (Protocol #0905-1).

Egg incubation and tissue collection

Eggs were washed in tepid water and candled to separate fertile and infertile eggs. Eggs from each clutch were randomly assigned to various experimental treatments, placed in plastic containers, and completely covered in moist vermiculite (1 part vermiculite:1 part water by mass). Egg containers were placed in foam box incubators set to a nominal temperature of 26.5°, which produces exclusively male hatchlings. Depending upon the particular study, eggs were either maintained at 26.5° throughout incubation or shifted to 31° for varying periods of time during the TSP. Embryos from each clutch were periodically sampled to monitor developmental rate and ensure embryos were at the desired stages for treatments [staging according to Yntema (1968)]. Equal (or roughly equal) numbers of eggs from each clutch were randomly assigned to each experimental group within each study. Thus, we were able to use clutch identity as a blocking factor/independent variable in all studies.

This experimental design is based on an extensive set of studies where eggs were shifted between different female- and male-producing temperatures at various stages of development to define the TSP for sex determination (Rhen et al. 2015). The entire TSP broadly encompasses Yntema stage 14 to stage 21 (Supplemental Material, Figure S1), but varies with the order of temperature shifts (male to female vs. female to male), the precise temperatures used (some temperatures are more potent at masculinizing or feminizing embryos than other temperatures), as well as the duration and timing of exposure to different temperatures. Those studies uncovered a brief sex-determining period in snapping turtles exposed to a specific thermal regime. Moving embryos from a male temperature (26.5°) to a female temperature (31°) for 5 days of a 65-day incubation period induces ovarian development in all stage 17 embryos. The sex-determining period for this thermal regime in snapping turtles (∼8% of embryogenesis) is comparable to the sex-determining period in mice (∼10% of embryogenesis) (Figure S1).

We used this temperature shift paradigm to produce mixed sex ratios instead of a constant pivotal temperature for two major reasons. First, we can focus on a brief period when high temperature (31°) first specifies (days 1–3 of the shift) and then determines (days 4–5 of the shift) ovarian fate (exposure starting at stage 17–17.5). Gonads in embryos kept at 26.5° during this same 5-day period are fated to become testes. This facilitates identification of causal TSD genes by allowing analysis of molecular changes precisely when temperature specifies gonad fate. In contrast, we would have to sample embryos over a 2-week period from stage 14 to stage 21 (the entire TSP) if eggs were incubated at a constant pivotal temperature. Second, incubation of eggs at a constant pivotal temperature throughout development is unnatural. Eggs in nature normally experience temperature fluctuations on various temporal scales (Wilhoft et al. 1983; Packard et al. 1985; Congdon et al. 1987; Kolbe and Janzen 2002). Thus, a temperature-shift paradigm is arguably a more natural way to produce mixed sex ratios and reveal phenotypic and genetic variation in temperature sensitivity.

Euthanasia of embryos and hatchlings was by rapid decapitation. Adrenal–kidney–gonad complexes were quickly dissected from embryos and hatchlings and either placed in RNAlater or fixed in 4% paraformaldehyde. Tissues in RNAlater were kept at −20° for long-term storage. Tissues in 4% paraformaldehyde were kept overnight at 4°, then washed in 1× PBS, and transferred to 70% ethanol for long-term storage.

RNA extraction from gonads and Complementary DNA synthesis

Gonads were microdissected from the underlying adrenal–kidney complex as previously described (Rhen et al. 2007). We then extracted total RNA from gonad pairs of individual embryos using the PicoPure RNA Isolation Kit (Life Technologies). A subset of samples from one study was extracted with the RNeasy Micro Kit (Qiagen) for comparison to the PicoPure Kit. Purified RNA (150 ng) from each gonad pair was reverse-transcribed in 20-µl reactions using the iScript cDNA Synthesis Kit (Bio-Rad, Hercules, CA). Gonad complementary DNA (cDNA) was diluted to the equivalent of 1.25 ng input RNA/µl for use in real-time PCR reactions. Thus, we measured gene expression in pure gonad preparations that were completely separated from adrenal and kidney tissue.

DNA extraction from kidneys

Genomic DNA was isolated from kidneys using TRI Reagent according to the manufacturer’s instructions (Molecular Research Center, Cincinnati). Genomic DNA was isolated from blood using a phenol/chloroform/isoamyl alcohol protocol available online at http://www.genome.ou.edu/protocol_book/protocol_partIII.html#III.H. The concentration and purity of DNA was measured with a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE). The integrity of DNA was determined by agarose gel electrophoresis. Purified genomic DNA was stored at −20° until used for genotyping.

PCR methods for measuring gene expression and for genotyping

We used dye-based and probe-based PCR to measure CIRBP expression and to determine CIRBP genotypes. Real-time PCR was performed using the Bio-Rad CFX 384 Real-Time PCR System with SsoFast EvaGreen supermix according to the manufacturer’s instructions (Bio-Rad). Each 10-μl PCR reaction contained 5 μl of the 2× EvaGreen supermix, 0.3 μl of the forward primer (0.3 μM final concentration), 0.3 μl of the reverse primer (0.3 μM final concentration), 2.4 μl of water, and 2 μl of cDNA (input = 2.5 ng total RNA). Standard curves were made via serial dilution of purified PCR products as described in Rhen et al. (2007).

We designed Taqman probes (Life Technologies) and primers to distinguish the SNP at position c63A > C in the open reading frame of snapping turtle CIRBP. SsoFast probes supermix (Bio-Rad) was used for allele-specific expression assays and genotyping. We isolated PCR products that contained only the A allele or only the C allele. We empirically determined the optimal annealing and extension temperature (64.3°) for distinguishing these alleles. We then made standard curves that were specific for each CIRBP allele: the A probe was hydrolyzed only when the A template was present (Figure S2A), while the C probe was hydrolyzed only when the C template was present (Figure S2B). Slopes and intercepts for standard curves were statistically homogeneous. We could therefore directly and quantitatively compare the abundance of transcripts for each allele and accurately genotype individuals.

We also designed primers for high-resolution melt temperature (HRM) analysis of PCR products to confirm CIRBP genotypes for the SNP at position c63A > C. We determined optimal conditions for distinguishing AA homozygotes, AC heterozygotes, and CC homozygotes using 10 ng of genomic DNA as template. Amplification was via three-step PCR. Samples were heated to 98° for 2 min and then subjected to 40 PCR cycles: 98° for 5 sec (denature), 53° for 1 sec (anneal), and 60° for 5 sec (extend). After 40 cycles, we collected HRM data from 60° to 85° at 0.2° increments and 10-sec read times. HRM data were imported into Bio-Rad Precision Melt Analysis Software. Optimal settings for CIRBP genotyping were determined using pure A standard as template, pure C standard as template, and an equal mix of A and C standards as template (premelt temperature = 70.4–70.5°; postmelt temperature = 75.7–75.8°; temperature shift bar height = 0.20; melt shape = 35; Tm = melt temperature difference = 0.25°).

In addition, we designed primers for HRM genotyping of loci presumably linked to CIRBP. The snapping turtle genome has not yet been sequenced so we used comparative genomics to identify seven protein-coding genes that display conserved gene order in human, chicken, painted turtle, and Chinese softshell turtle genomes: STK11, ATP5D, MIDN, CIRBP, C19orf24, EFNA2, and CAMK4L. Given the synteny across mammals, birds, and divergent turtle species, these genes are most likely linked in the snapping turtle. We used Illumina sequencing data from the snapping turtle to identify SNPs within these genes (more details are provided in Sequencing and Bioinformatics, SNP identification, and RNA-Seq below). We optimized HRM genotyping for these SNPs as described above. Sequences for PCR primers and probes are provided in Table S1.

Temperature shifts to characterize CIRBP mRNA and protein expression patterns

We conducted an experiment to define the time course for CIRBP mRNA induction. We collected gonads from embryos incubated at 26.5° and from clutch mates shifted from 26.5° to 31° at stage 16. Approximately equal numbers of embryos from each temperature were sampled at 6, 12, 24, and 48 hr of the temperature shift. Tissue was collected and RNA extracted from embryonic gonads as described above.

A similar experiment was conducted to examine CIRBP protein expression during and after the sex-determining period. Eggs were incubated at 26.5° until embryos reached stage 16. At that time, half of the eggs were kept at 26.5° to produce males, while half were shifted to 31° for 6 days to produce females. Embryos were sampled from both temperatures on days 1, 2, 3, 4, 5, and 6 of the shift. We also sampled embryos after eggs at 31° were returned to 26.5° (days 11, 21, 31, and at hatch). Adrenal–kidney–gonad complexes were fixed in 4% paraformaldehyde for 24 hr at 4° and then transferred to 70% ethanol for long-term storage.

CIRBP expression–sex ratio correlation study

We conducted an experiment to examine the relationship between CIRBP mRNA expression in gonads during the TSP and sex ratio at hatching. Snapping turtle eggs from nine randomly collected clutches were incubated at 26.5° until embryos reached stage 17.5 of development. At stage 17.5, eggs were shifted to 31° for 2.5 days and then shifted back to 26.5° for the rest of development. This is the developmental stage when snapping turtle embryos are most sensitive to the feminizing effect of high temperatures (Rhen et al. 2015). Roughly equal numbers of embryos from each of the nine clutches were sampled at 24 and 48 hr of the shift to 31°. Tissue was collected and RNA was extracted from embryonic gonads as described above. A subset of eggs from each clutch was allowed to hatch to determine sex ratios for each clutch. The sex of hatchlings was diagnosed based on gross morphology of gonads and reproductive tracts as previously described (Rhen and Lang 1994).

Sequencing

We used cDNA from the gene expression–sex ratio correlation study described above as template for PCR amplification of the entire CIRBP-coding sequence in multiple individuals. Purified PCR products were Sanger-sequenced using BigDye as previously described (Rhen et al. 2007). Sequencher software was used to align sequences and search for SNPs.

We also collected 12 clutches from southern Minnesota and 13 clutches from northern Minnesota for RNA-Seq (range in latitude = 43.511° North to 48.720° North). Approximately equal numbers of eggs from these clutches were separated into two experimental groups. One group of eggs was incubated at 26.5° throughout development to produce males. The other group of eggs was incubated at 26.5° until stage 17.5, and eggs were then shifted to 31° for 6 days to produce females and returned to 26.5° for the rest of development. Embryos from both groups were sampled on days 1, 2, 3, 4, and 5 of the temperature shift. Tissue was collected and RNA was extracted from embryonic gonads as described above.

We combined an equivalent amount of total RNA from gonad pairs from individual embryos to make 20 pools for RNA-Seq: 2 temperatures × 5 days × 2 biological replicates. Total RNA from 25 embryos was pooled for each biological replicate. All together, we sampled 500 embryos from 25 clutches collected across a wide geographic range to capture as much genetic variation as possible (5° range in latitude from the southern-most clutch to the northern-most clutch within Minnesota). Ten pools (2 temperatures × 5 days) were from eggs collected in southern Minnesota (12 clutches). The other 10 pools (2 temperatures × 5 days) were from eggs collected in northern Minnesota (13 clutches). Clutches of eggs from northern and southern Minnesota were separated by latitude of 46° 55′ N. Total RNA was sent to the W. M. Keck Center for Comparative and Functional Genomics at the University of Illinois at Urbana-Champaign for sequencing. Twenty individual cDNA libraries were prepared with unique barcodes and sequenced on Illumina HiSeq 2000 (100 bp, single-end reads; 156.4 million reads).

Additional sequence data using tissues derived from all three germ layers enabled de novo assembly of a more complete snapping turtle transcriptome. The data described above were supplemented with sequence data from eight hypothalamus/pituitary libraries (50 bp, single-end Illumina reads; 172.2 million reads), two intestine libraries (50 bp, single-end Illumina reads; 31.8 million reads), as well as normalized embryonic gonad libraries (454 Sequencing, average read length = 350 bp; 2.8 million reads).

Individual and population-based genetic association studies

Snapping turtle eggs from another 15 clutches from northern Minnesota were incubated at 26.5° until stage 17.5 of development. Eggs were shifted to 31° for 2.5 days and then shifted back to 26.5°. Similar numbers of embryos were sampled from each clutch at 24 and 48 hr of the exposure to 31°. Tissue was collected and RNA was extracted from embryonic gonads as described above. A subset of eggs from each clutch was allowed to hatch to determine sex ratios and to collect tissue for genotyping. Hatchling sex was diagnosed, adrenal–kidney–gonad complexes were dissected, snap-frozen, and stored at −80° until genomic DNA was isolated from kidneys. Hatchlings were genotyped using TaqMan MGB probes. Genotypes were verified using a different set of PCR primers and high-resolution melt temperature analysis.

Equal numbers of adult turtles were collected from the southern half and the northern half of Minnesota to compare allele frequencies in meta-populations experiencing gene flow, but differing slightly in TSD pattern. A small number of turtles from Texas were available for comparison to a geographically distinct population unlikely to experience gene flow with turtles in Minnesota. We drew blood from the subcarapacial vein of adult turtles as described by Moon and Hernandez Foerster (2001). Whole blood was transferred to a microfuge tube and kept on ice until genomic DNA was extracted (within a few hours).

We genotyped adult snapping turtles captured in northern Minnesota (n = 28), southern Minnesota (n = 28), or northern Texas (n = 9). Genotypes for 481 individual embryos and hatchlings were also available from 15 clutches collected in northern Minnesota (i.e., from the genetic association study). To avoid pseudoreplication (i.e., siblings within a clutch are not independent samples), we used genotype frequencies within families to infer parental genotypes. Most clutches (n = 14) exhibited genotype frequencies consistent with Mendelian inheritance from a single mother and a single father. One clutch displayed genotype frequencies, suggesting multiple paternity. Thus, we inferred 31 genotypes for parents of embryos and hatchlings from the genetic association study.

Bioinformatics, SNP identification, and RNA-Seq

We used CLC Genomics Workbench (CLC bio, Cambridge, MA) for read mapping, SNP verification, and RNA-Seq analyses. We mapped reads to full-length CIRBP cDNA and determined the number of reads that were mapped and not mapped from each cDNA library (20 total gonad libraries). We used the SNP detection tool to identify polymorphisms in CIRBP cDNA. Parameters for SNP detection required a minimum coverage of 20 reads for each SNP and a minimum variant frequency of 10%. Illumina reads uniquely mapped to CIRBP were then used to estimate allele frequencies in populations from northern and southern Minnesota. A recent article shows that pooled RNA-Seq data can be used to estimate allele frequencies with accuracy similar to pooled genome sequencing (Konczal et al. 2013). As indicated in the previous section, we also conducted a completely independent study in which we extracted DNA and genotyped adult and hatchling turtles using traditional methods.

The same procedure was used to identify SNPs in a set of genes flanking CIRBP. The snapping turtle genome has not yet been sequenced so we used comparative genomics to identify genes likely to be linked to CIRBP. Six protein-coding genes are syntenic with CIRBP in human, chicken, painted turtle, and Chinese softshell turtle genomes in the following chromosomal arrangement: STK11, ATP5D, MIDN, CIRBP, C19orf24, EFNA2, and CAMK4L. These genes span ∼650 kb in a painted turtle scaffold (gi|371559100|gb|JH584733.1) and exhibit highly correlated physical distances in the homologous Chinese softshell turtle scaffold (gi|351669362|gb|JH209284.1) (r = 0.9988, n = 7). We used sequence data from snapping turtle gonads, hypothalamus/pituitary, and intestine to identify 12 common SNPs in these genes.

We used Illumina single-end reads to measure allelic-specific CIRBP expression. We normalized gene expression using reads per kilobase of exon per million mapped reads (RPKM) (Mortazavi et al. 2008). We explicitly estimate allele-specific expression by counting only those reads that contain the SNP (Jiang and Wong 2009). We calculated expression values for alternative alleles using uniquely mapped reads (n = 4,744) that included position 63 within the CIRBP-coding sequence. We used 199 bp as the length of the transcript for allele-specific RPKM calculations because Illumina reads from gonads were 100 bp long and had to overlap the SNP to be informative. Other methods, such as RNA-Seq by Expectation-Maximization, use uninformative reads to calculate expression values for splice isoforms (or alleles) by assigning fractions of uninformative reads in proportion to counts from informative reads (Li et al. 2010). The assumption that uninformative reads are from a particular isoform (or allele) is unfounded and could lead to biased expression estimates for splice variants and for alternative alleles. Thus, we have chosen to be conservative and only use informative, uniquely mapped reads for our analysis of allelic-specific expression. Two other SNPs were detected in the 3′ UTR of CIRBP. Although we can unambiguously map reads containing these SNPs, we cannot determine linkage relationships among SNPs separated by >100 bp from single-end RNA-Seq data. We therefore analyzed each SNP separately.

Linkage analysis

We used PCR and HRM analyses to determine genotypes for markers in STK11, ATP5D, MIDN, CIRBP, C19orf24, EFNA2, and CAMK4L in 56 adult snapping turtles from Minnesota. We assigned numerical values of 0, 1, or 2 based to the genotype at each locus: homozygotes for rare allele = 0, heterozygotes = 1, and homozygotes for common allele = 2. We excluded ATP5D from further analysis because the frequency of the rare allele was too low (0.027) to be useful as a marker. We calculated pairwise correlations among genotypes for all other loci and squared this value (r²) to get a standardized measure of LD as outlined by Rogers and Huff (2009).

Immunohistochemistry

Paraformaldehyde-fixed tissues were embedded in paraffin, sectioned at 6 µm, and placed on microscope slides. Sections were de-paraffinized and rehydrated using standard immunohistochemistry protocols. The primary CIRBP antibody was diluted 1:500 (caalog no. ARP40348_P050, Aviva Systems Biology, San Diego). Negative controls were incubated in blocking solution without the primary antibody. The VectaStain ABC Elite kit (Vector Laboratories) was used to detect the primary/secondary antibody complex. Staining was performed using fresh 3, 3′-diaminobenzidine (Vector Laboratories) and counterstained with hematoxylin. Images were captured with an Olympus BX-51 microscope equipped with an Infinity 2 digital camera (Lumenera, Ottawa) using Rincon HD software (Imaging Planet, Goleta, CA).

Statistical analyses

We used three-way analysis of variance (ANOVA) with clutch, temperature, and day as main effects to analyze CIRBP mRNA expression in the quantitative PCR (qPCR) study. Significant main or interaction effects were followed by comparisons among specific groups of interest. The Dunn-Sidák method was used to correct for multiple comparisons. The nominal significance level was calculated as α’ = 1 − (1 − α)^1/k, where k is the number of comparisons for an experiment-wise α = 0.05. Spearman’s rank correlation was used to test for a correlation between CIRBP mRNA expression and sex ratio in clutch mates. This nonparametric test is appropriate when data do not meet assumptions for the Pearson product-moment correlation (i.e., a linear relationship between variables). We used likelihood-ratio (LR) chi-square tests or Fisher’s Exact Tests for variation in sex ratios among clutches and for associations between genotype and sex. In another test for association between CIRBP genotype and sex, we used a matched-pairs t-test to control for genetic background: allele frequencies in males and females were compared after pooling clutches with very similar numbers of males and females (Risch and Teng 1998; Teng and Risch 1999). We used the Cochran–Mantel–Haenszel test to compare allele frequencies between northern and southern Minnesota populations while controlling for temperature effects on transcript levels in the RNA-Seq study. In other words, data were stratified by temperature, thus allowing comparison of northern and southern Minnesota populations. Statistics for gene expression and genetic associations were performed using JMP 5.0.1.2 software (SAS Institute, Cary, NC). Effects were considered significant at P ≤ 0.05. Throughout the text, means are given ± SEM (± SE).

Data availability

Data files are uploaded as supplementary files: Data for Figure 1 contain raw and processed expression data for Figure 1. Data for Figure 3 contain clutch sex ratios and mean expression data for Figure 3. Data for Figure 3 contain raw expression data for Figure 3. Data for Figure 4A contain read counts for alternative CIRBP alleles from RNA-Seq study and summary data plotted in Figure 4A. Data for Figure 4B contains raw expression data and summary data plotted in Figure 4B. Data for Figure 4C, Table 2, and Figure S2 contain individual genotype and phenotype data for those figures and table. Data for Figure 5 contain raw data presented in Figure 5. Data for Figure 6 and Figure 7 contain individual genotype data for those figures.