Electronic Supplemental Online Material

ELECTRONIC SUPPLEMENTAL ONLINE MATERIAL

TITLE: Body size evolution simultaneously creates and collapses species boundaries in a clade of scincid lizards

AUTHORS: Jonathan Q. Richmond and Elizabeth L. Jockusch

(a.) Models and priors for inferring tree topology and divergence times

The phylogeny in Figure 1 is based on a dataset expanded from Richmond and Reeder (2002) and was inferred using the following procedures. Nucleotide substitution models were determined using the Akaike Information Criterion in MrModelTest v2.2 (Nylander 2004). The tRNAs were treated as a single data partition because of their small size (206 nt combined) and similar functional constraints, and the ND4 region was treated as a second partition. The models selected were HKY + I + Γ and GTR + I + Γ for the tRNA and ND4 partitions, respectively. The priors used in Mr. Bayes 3 (Ronquist & Huelsenbeck 2003) were as follows: rate parameters uniform [0,1000]; alpha exponential (1.0); proportion of invariant sites uniform [0,1]; branch lengths exponential (10.0). Two analyses were run for 10,000,000 generations each, with samples taken every 1000th generation. Trace plots were generated in Tracer v1.3 (Rambaut & Drummond 2005) to assess stationarity and effective sample sizes. Because the model parameters converged to similar stable values, trees were combined from each run to create a consensus phylogram after eliminating burn-in samples.

Divergence times for haplotypes of each pair of populations used in the mate compatibility trials and for larger clades of interest were estimated using the following models and priors in Beast v1.3 (Drummond & Rambaut 2003): GTR + I + Γ nucleotide substitution model and a UPGMA starting tree; rate parameters uniform [0,500]; alpha exponential (1); proportion of invariant sites uniform [0,1]; discretized branch rates uncorrelated log normal prior (1, 0.1). The tree prior was derived from a non-linear skyline model with ten groups, and defaults were implemented for all scale operators. Nucleotide substitution rates were not correlated between parent and daughter branches on the phylogeny (8.25 x 10-3 analysis: covariance mean = -0.004, 95% CI = -0.11 – 0.10; 7.15 x 10-3 analysis: covariance mean = -0.003, 95% CI = -0.10 – 0.11), indicating that the use of an uncorrelated relaxed clock model was appropriate for estimating divergence times. Convergence and effective sample sizes were monitored in Tracer v1.3. Analyses were run for 10–18 million generations, with parameters sampled every 1000 generations. The results for both calibration rates are shown in Table ESM1.

(b.) Husbandry and experimental conditions

Skinks were housed individually in a climate-controlled room (10/14 hr dark/light cycle) and fed crickets dusted with a calcium-phosphorus supplement every 3 days. Testing chambers consisted of 40.0 X 76.0 X 40.0 cm glass terraria with open lids and full spectrum lighting. Three sides of the inner glass were lined with black shelving paper to eliminate reflections and to minimize outside visibility. Skinks were separated by a divider during a 10 minute acclimation period prior to each trial, and a pulley system was used to remove the divider at the onset of testing. All trials were conducted behind a blind and lasted for 45 minutes, or until copulation was complete.

(b.) Statistical analysis of mate compatibility and morphological data

The MCMC specifications for evaluating model fit and estimating regression coefficients adhered to the same general format in all WinBUGS v1.4.1 (Spiegelhalter et al. 2003) analyses. Four Markov chains were initiated from overdispersed starting values, and trace plots generated in WinBUGS were used to confirm that model parameters had reached stable values. Stationarity was also evaluated using the Gelman-Rubin statistic, which is close to 1.0 if the chains have converged properly (Gellman & Rubin 1992). Analyses were run for 100,000 iterations, and elimination of the first 10% of the samples provided conservative burn-in for all analyses. Samples were often thinned to retain every 10th to 50th cycle of the Markov chain, and thinning adjustments were made according to the autocorrelation monitor in WinBUGS. We specified vague normal priors for which all plausible values have nearly equal probability density.

By including intercept effects for each individual and trial series, the best-fit model for the mate compatibility experiments accounts for variation across individuals and trial series in the propensity to mate. This controls for the non-independence of trials resulting from the hierarchical trial series design and repeated use of individuals in estimating the intercept coefficient. A model allowing a separate slope coefficient for each trial day (i.e. controlling for non-independence in estimating the slope coefficient) failed to converge. However, we have two reasons to believe that the large effect of body size divergence is biologically real, and not a result of the repeated use of individuals. First, we note that this coefficient is estimated for pairs of individuals, and each pair occurs only once in the dataset. Second, we estimated the slope coefficient for three non-overlapping subsets of data in which each individual was used at most once. These subsets were produced by using data from only a single day (days 1-3) in each trial series. In all cases, the mean value of the body size coefficient was negative, and the 95% CIs excluded 0, supporting the conclusion that the probability of copulation decreases as body size divergence increases.

In cases where hormone supplements were needed, a pilot study indicated that two intraperitoneal injections of testosterone propanoate (~0.003µg/gram body weight for males) or estradiol propionate (~0.05µg/gram body weight for females) administered two days apart consistently induced receptivity. To test whether our results were consistent between hormone and no-hormone trials, we compared DIC scores for a model in which the intercept was allowed to vary between trials with and without hormones, and one in which it did not. Inclusion of this additional intercept parameter did not improve model fit (Table ESM2). Thus, there is no evidence that use of hormones induced results that differ from those in naturally receptive individuals.

We also tested whether skinks were more or less likely to mate as a trial series progressed. To perform this analysis, a ‘trial order’ parameter was scored categorically (1-5), where numbers indicated the day of a trial over the five-day series, and allowed the intercept to vary between trial days. Our results show that inclusion of a trial order effect on the intercept did not improve model fit (Table ESM2), indicating that individuals did not suffer from trial fatigue because they were equally likely to mate on any given day of a trial series.

(c.) AFLP protocol

AFLPs were generated from genomic DNA using the restriction enzymes EcoR1 and Bfa1 and ligation of oligonucleotide adaptors. These reactions were followed by two rounds of selective PCR to reduce the total pool of AFLP fragments. The first round incorporated primers that matched the adaptor sequence plus one additional nucleotide (EcoR1A x Bfa1C – subscripts indicate the selective nucleotides), and the second incorporated primers that matched the adaptor sequence plus three additional nucleotides (EcoR1ACT x Bfa1CGC; EcoR1AGC x Bfa1CAG). The fluorescein-labeled dyes HEX and FAM were attached to the EcoR1 primers in the second round of PCR for genotyping.

Because AFLPs are dominant markers, calculating classical F-statistics (Wright 1951) using standard procedures is impossible without knowledge of the inbreeding coefficient within populations. However, the Bayesian approach used in this study allows for the direct application of Wright’s F-statistics by incorporating uncertainty into the measure of within-population inbreeding and relaxing the assumption of Hardy-Weinberg equilibrium (Holsinger et al. 2002; Holsinger & Wallace 2004). To estimate the posterior distribution of qB in Hickory v1.0.3 (Holsinger & Lewis 2004), four Markov chains were run for 500,000 generations. Chains were thinned by drawing samples every 50th iteration. Trace plots showed that exclusion of the first 50,000 samples provided ample burn-in.

Table ESM1: Divergence dates estimated in Beast v1.3 using two calibration rates. ESS is the effective sample size taking into account autocorrelation across generations. Mean divergence times, standard deviations and 95% CIs are reported in millions of years. Calibration rates are in substitutions/site/million years. Numbers for skiltonianus refer to the clades identified in Fig. 1.

Calibration rate → / 0.00825 / 0.00715
Lineage / mean ± s. d. / 95% CI / ESS / mean ± s. d. / 95% CI / ESS
All gilberti / 12.70 ± 0.14 / 8.79 – 17.00 / 240 / 14.70 ± 0.09 / 10.27 – 19.40 / 744
Southwestern-Inyo gilberti / 5.47 ± 0.05 / 4.04 – 7.12 / 227 / 6.31 ± 0.06 / 4.67 – 8.14 / 264
Southwestern gilberti-skiltonianus 8 / 4.55 ± 0.04 / 3.20 – 6.04 / 322 / 5.18 ± 0.05 / 3.64 – 6.87 / 380
Inyo gilberti- skiltonianus 6 / 1.72 ± 0.02 / 1.07– 2.47 / 376 / 2.01 ± 0.02 / 1.17 – 2.87 / 628

Table ESM2: Regression models used to explore the effects of hormone injections and order of testing within a trial series.

Table ESM3: Regression coefficients for size coefficients when allowed to vary across trial type (Table 1, model B; s = skiltonianus, g = gilberti)

Coefficient / Mean ± S.D. / 95% CI
β(1) s-s / -0.38 ± 0.07 / -0.52 – -0.25
β(2) s-g / -0.48 ± 0.09 / -0.69 – -0.33
β(3) g-g / -0.31 ± 0.07 / -0.47– -0.18

Figure ESM1. Experimental design for a five day trial series. (A.) Ovals represent mtDNA clades A-E, where A-C represents the three gilberti clades (subscript g) and D-E represent skiltonianus clades (subscript s). In this example, a gilberti male from clade A would be mated with a female from clade A on day one of the series. On day two, that same male would be mated with a female from clade B, on day three he would mate with a female from clade C, et cetera for a total of five days. The order of pairing of test individuals, as well as testing order for each day, was randomized. (B.) The 5 X 5 matrix showing each of the pairwise matings for a trial series (1 = successful copulation, 0 = no copulation).

REFERENCES

Drummond, A. J. & Rambaut, A. 2003 BEAST: Bayesian Evolutionary Analysis Sampling Trees v1.3: http://evolve.zoo.ox.ac.uk/beast/.

Gellman, A. & Rubin, D. B. 1992 Inference from iterative simulation using multiple sequences. Stat. Sci. 7, 457-511.

Holsinger, K. E. & Lewis, P. O. 2004 Hickory 1.0.3: http://darwin.eeb.uconn.edu/hickory/hickory.html.

Holsinger, K. E., Lewis, P. O. & Dey, D. K. 2002 A Bayesian approach to inferring population structure from dominant markers. Mol. Ecol. 11, 1157-1164.

Holsinger, K. E. & Wallace, L. E. 2004 Bayesian approaches for the analysis of population genetic structure: an example from Platanthera leucophaea (Orchidaceae). Mol. Ecol. 13, 887-894.

Nylander, J. 2004 MrModeltest 2.2: http://www.csit.fsu.edu/~nylander/.

Rambaut, A. & Drummond, A. J. 2005 Tracer - v.1.3: http://evolve.zoo.ox.ac.uk/software.html?id=tracer.

Ronquist, F. & Huelsenbeck, J. P. 2003 Mr. Bayes3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572-1574.

Spiegelhalter, D. J., Thomas, A., Best, N. & Lunn, D. 2003 WinBUGS (Bayesian Inference Using Gibbs Sampling) vers. 1.4.1: http://www.mrc-bsu.cam.ac.uk/bugs/welcome.shtml.

Wright, S. 1951 The genetical structure of populations. Ann. Eugen. 15, 323-354.