Genotypic Sex Determination and Analysis

Genotypic Sex Determination

Genotypic sex (ZZ or ZW) was determined using a PCR test on DNA from tail tips or blood.
The PCR test followed the method described by Holleley et al. (2015) using sex-linked markers identified by Quinn et al. (2010).
PCR products were visualized on a 1.5% agarose gel.
ZW females had two bands, while ZZ individuals had only one band.
Internal ZZ and ZW control samples were run on every gel.
Animals with genotype-phenotype discordance were classified as sex-reversed, according to Holleley et al. (2015).
Sex-reversal was reported in two ways:
- Sex-reversal rate: proportion of all ZZ individuals that are female ( $pZZ$ ) (Grossen et al., 2011; Schwanz et al., 2020).
- Proportion of all females that are sex-reversed ( $p_f$ ) (Holleley et al., 2015; Castelli et al., 2021).
Chi-square tests were used to determine the association between sex-reversal ( $pZZ$ ) and sampling year, aggregated into 8-year periods (2003-2011 and 2012-2019).

Genotype by Sequencing and Population Genetic Analysis

Reduced representation sequencing (Kilian et al., 2012) was used to genotype individuals.
This was done to determine levels of population genetic structure (as an indicator of gene flow) across the study area.
Genomic DNA was extracted from tail tissue and blood samples using:
- Qiagen's Gentra Puregene DNA purification kit for tail tips.
- Whatman Elute quick extraction protocol for blood on Whatman FTA Elute Cards.
DNA was extracted, sequenced, and informative SNP markers identified using Diversity Arrays Technology (DART Pty Ltd.).
SNP genotyping was performed using complexity reduction with restriction enzymes, implicit fragment-size selection, and next-generation sequencing (Kilian et al., 2012).
The restriction enzyme combination of PstI (5'-CTGCA | G-3') and SphI (5'-GCATG | C-3') was used for complexity reduction by double digestion of genomic DNA.
Sequences were processed using proprietary DART analytical pipeline (Kilian et al., 2012) to yield polymorphic SNPs.
Calling quality was assured by high average read depth per locus (medium coverage, average ca 10×).
Approximately one-third of samples were processed twice from DNA to allelic calls as technical replicates.
Scoring consistency (repeatability) was used as the main selection criteria for high-quality and low error-rate markers.
For further detail refer to Georges et al (2018).
Sequences were then filtered using the R package dartR (v.1.1.11; Gruber et al., 2018).
Loci were removed with a call rate of less than 95%.
Loci were removed with repeatability less than 99%.
Loci were removed if they were sex-linked across all individuals.
Secondary SNPs residing on a single sequence tag were filtered, retaining only one at random.
After filtering loci, individuals with a call rate of less than 95% were removed.

Population Genomic Analyses

The R package dartR (Gruber et al., 2018) was used for population genomic analyses of the DARTseq SNP data.
Principal coordinates analysis (PCoA) and F statistics ( $F<em>{IS}$ , $F</em>{IT}$ , $F_{ST}$ ) were used to identify population structure across the study area.
Isolation by distance was used to assess the relationship between genetic and geographic distance.
Genetic distance between demes was estimated using $F_{ST}$ among the sexes and using the proportion of shared alleles between pairs of individuals.
A Mantel test in the R package adegenet (Jombart, 2008) was used to test for an association between genetic distance and geographic distances.
Pairwise analysis of $F_{ST}$ was used to examine populations for sex-biased dispersal.
A redundancy analysis (RDA) was applied to decompose the contribution of geographical coordinates (latitude and longitude) and sex on allele frequencies (Meirmans, 2015).

Study Site, Frequency of Sex-Reversal, and Morphometrics

An exhaustive survey of I. vitticeps was conducted at one site within the area described in Section 2.1.
A population of P. vitticeps was studied continuously from September 2018 to March 2020 at Bowra Wildlife Sanctuary (a 14,000 ha property).
Operational sex ratio and rates of sex-reversal were quantified.
Vegetation type at Bowra was determined using a stereoscopic photo interpretation of 2015 aerial photography (1:30,000, resolution 0.21 m), with ground truthing by vegetation surveys.
Vegetation was classified into six types (woodland, grassland, mixed forest woodland, claypan/herbland, shrubland, and other) based on dominant vegetation (Table S1).
Surveys for lizards were conducted by foot and by vehicle within the property, timed to match the diurnal activity patterns of the lizards during their most active months (October-December).
Foot and road surveys were standardized by the number of hours and observers.
Foot surveys included all vegetation types at least once each week and all roads driven at least once each week (Figure S1).
Capture effort for each lizard was estimated as the sum of total search hours (or driving hours) divided by the total number of lizards captured.
Capture effort was then compared by survey method (search/vehicle count) and compared across sex classes.
Chi-square tests were used to estimate if sex-ratio was biased in the population between males (ZZm) and females (ZZf +ZWf).
Upon capture, each lizard was sexed phenotypically and a blood sample taken for genotyping (following procedures outlined in Section 2.1).
Snout-vent length (SVL) was measured using a ruler (±1 mm) and mass was measured using an Optek digital scale (±0.5 g).
Body condition index (BCI) was calculated for each lizard as the residual from an ordinary least squares (OLS) regression of SVL and mass (after Schulte-Hostedde et al., 2005).
ANOVA was used to compare how morphometric data (SVL, mass, and body condition index) varied among sex classes.

Radiotelemetry and GPS Tracking

Adult lizards were tracked from September 2018 to March 2020 to determine how movement rates, space use, survival, and growth vary among sex classes.
Each lizard was fitted with a Pinpoint Beacon 250 (Lotek Ltd.) using a custom-fit backpack harness.
Each unit housed a GPS logger, a single stage VHF transmitter (150-151 Hz) and 3-axis accelerometer.
The unit and harness did not exceed 5% of the animal's body mass.
Dilution of precision (DOP; Langley, 1999) was recorded for every GPS fix as an indicator of quality of GPS position (≤5 m accuracy).
Estimates of movement and space use were calculated from GPS data only (DOP ≤3).
Lizards were located using a VHF receiver and Yagi antenna.
GPS points were acquired every day in spring (September-November), summer (December-February), and autumn (March-May), and once every three days in winter (June-August).
The status of animals (dead or alive) was determined by relocating all individuals once per week (September-December) and a minimum of once per month for the remainder of the year (January-August).
Females were palpated every two weeks during the breeding season (October-December) and monthly when lizards were handled to change telemetry equipment, to determine reproductive status.

Movement and Space Use Estimates

Movement rate (metres moved per day, m/d) was calculated as the straight-line distance between consecutive GPS locations divided by the number of days elapsed between relocations.
A linear mixed-effects model was used to test the effect of season, sex class, and their interaction on movement rate with individual lizard as a repeated (random) effect.
Within-individual variability of movement between gravid and non-gravid periods was compared by recording distance moved per day for gravid females (2 weeks pre/post lay date) with the distances moved while non-gravid (2-6 weeks post lay date).
The effect of reproductive status on distance moved per day (response variable) was tested, with reproductive status (gravid or not gravid) as a fixed effect, and female ID as a repeated effect.
Two home range estimates were used to compare space use:
- 100% minimum convex polygons (MCP; Burt, 1943; Row & Blouin-Demers, 2006).
- 95% fixed kernel density estimators (KDE; Worton, 1989; Young et al., 2018) using the adehabitatHR package (Calenge, 2006).
Least square cross validation was used to select a kernel-smoothing parameter for KDE.
Space use was estimated only for animals that had a minimum of 15 GPS locations, and estimates were log10-transformed for both 100% MCP and 95% KDE to fit normality assumptions.
A linear mixed-effects model was used to test the effect of season, sex class, body size (SVL) and their interaction on space use estimates with individual lizard as a repeated (random) effect.
Equations for each movement and space use model used can be found in Supplementary Table S2.

Survival and Growth Rate Estimates

In cases where a tracked lizard was found dead, the date of death was determined using the last movement recorded by the accelerometer housed on the GPS units or abnormal large-scale movements (>1 km) over a short period of time.
Maximum likelihood survival probabilities were estimated using known fate models in the program MARK (White & Burnham, 1999).
Akaike Information Criteria (AICc) was used to correct for small sample sizes and models with delta AICc of <2.0 from the best model were considered to have support.
Growth rates were calculated by dividing the change in SVL between the initial capture and subsequent re-captures by the total number of days elapsed (>20 days).
Differences in SVL growth rates across sex classes were determined using an analysis of covariance (ANCOVA) with sex class as a fixed (3-level) factor, growth rate as the response variable and initial log10-transformed SVL as a covariate.

Results

Frequency of sex-reversal across populations

A total of 207 lizards were captured during sampling trips across the study area (excluding individuals captured at Bowra).
Once individuals were genotyped for their chromosomal sex this amounted to:
- ZZf (n = 18)
- ZWf (n = 98)
- ZZm (n = 91)
Sex-reversal was documented in three out of four demes that were sampled.
Overall, the proportion of females ( $p_f$ ) that were sex-reversed was 16%.
The rate of sex-reversal ( $pZZ$ ) was 17% when considering all animals captured during the 16-year sampling period.
A decline in the rate of sex-reversal ( $pZZ$ ) from 24% to 6% ( $\chi^2$ = 5.36, df = 1, p < .05) was observed between the first half of that sampling period (2003-2011) to the second (2012-2019).