In accordance with predictions of the size-advantage model, comparative evidence confirms that protogynous sex change is lost when mating behavior is characterized by weak size advantage. However, we lack comparative evidence supporting the adaptive significance of sex change. Specifically, it remains unclear whether increasing male size advantage induces transitions to protogynous sex change across species, as it can within species. We show that in wrasses and parrotfishes (Labridae), the evolution of protogynous sex change is correlated with polygynous mating, and that the degree of male size advantage expressed by polygynous species influences transitions between different types of protogynous sex change. Phylogenetic reconstructions reveal strikingly similar patterns of sex allocation and mating system evolution with comparable lability. Despite the plasticity of sex determination mechanisms in labrids, transitions trend towards monandry (all males derived from sex-changed females), with all observed losses of protogyny accounted for by shifts in the timing of sex change to prematuration. Likewise, transitions in mating system trend from the ancestral condition of lek-like polygyny toward greater male size advantage, characteristic of haremic polygyny. The results of our comparative analyses are among the first to confirm the adaptive significance of sex change as described by the size-advantage model.

Trait data compilation

Data on species-specific mating systems and sex allocation pathways were compilied from the primary literature (file: TableS5.xlsx, references are listed within). The types of sex allocation and mating systems expressed by labrid fishes, their definitions and predicted associations with the degree of male size advantage are detailed in figure 1. Mating system classifications focused only on terminal phase males and did not consider the mating strategies of initial phase males – although it is known that the reproductive output of initial phase males can outweigh that of terminal phase males for some species dependent on location-specific population dynamics (Warner and Hoffman 1980a,b; Warner 1982). We applied the consensus classification of the predominant mating system (i.e. supported by multiple authors) whenever possible, and otherwise relied on the most recent observations. We restricted sexual ontogeny data to accounts of protogyny that were distinguishable as either monandric or diandric based on gonad histology, population demographics or both. Cases where males are derived from females that have not passed through a functional stage were categorized as functionally gonochoristic following previous work (Sadovy and Shapiro 1987; Sadovy de Mitcheson and Liu 2008; Kazancıoğlu and Alonzo 2010; Erisman et al. 2013). Mating system and sex-change data were available for 89 labrid species. 

Phylogenetic inference and divergence time estimation

Gene sequences were aligned separately in Geneious Pro R8.1.9 (http://www.geneious.com/) using default settings and the alignments were manually adjusted through the insertion or deletion of gaps and trimmed to minimize the amount of missing data. Alignments were concatenated and partitioned by gene region, with separate partitions for the 3^rd codon position of protein coding genes.

Bayesian inference (BI) using partitioned mixed models in MrBayes v3.2.3 (Ronquist and Huelsenbeck 2003) was implemented on the CIPRES Science Gateway portal (Miller et al. 2010) to estimate the tree topology and branch lengths. We sampled across the entire general time reversible (GTR) model space using reversible jump Markov chain Monte Carlo (rjMCMC) methods to integrate out uncertainty about the correct substitution model for each partition (nst = mixed). This allowed us to quantify posterior probabilities of the substitution models sampled (i.e. the probability of each substitution model conditional on the data; table S3). Parameters were unlinked across partitions. Substitution rates and stationary nucleotide frequencies were allowed to evolve under different rates using a flat Dirichlet prior. The shape of the gamma distribution of rate variation evolved under an exponential prior with a mean of one. Branch lengths were unconstrained under an exponential prior with a rate of 200 (mean = 0.005). Posterior probabilities of clades were calculated following two 40 million generation Markov chain Monte Carlo (MCMC) analyses, each with eight chains (temp = 0.02) and two swaps, sampling every 2,000 generations. Convergence was assessed in Tracer v1.5 (Rambaut et al. 2009). Upon examination of the trace files, a conservative burn-in of 45% was discarded from each run and a 50% majority-rule consensus tree was computed using the remaining sampled trees (file: MrBayes_ConsensusTree.nex).

To estimate the temporal component of evolution, we first converted the majority-rule consensus tree to a rooted, ultrametric tree using the chronos function in the R package ape v3.2 (Paradis et al. 2004). We set the smoothing parameter lambda to 0.9 and specified a relaxed model of substitution rate variation among branches and seven age constraints (table S4). We then performed a divergence time analysis in BEAST v1.8.1 (Drummond et al. 2012) using the resultant ultrametric tree as the starting tree. Partitioning followed the scheme above and models of molecular evolution were specified using the parameters and priors corresponding to the model with the highest posterior probability from the MrBayes analysis and empirical base frequencies (see table S3 for substitution model details). Divergence times were estimated under a relaxed uncorrelated lognormal clock model (Drummond et al. 2006) and the birth-death process (Gernhard 2008). Evidence from six fossils informed exponential priors on corresponding nodes (table S4) to time-calibrate the trees. Exponential priors were used because they include the minimum age of the fossil and accommodate uncertainty in the age of the clade relative to the age of the fossil with soft upper-bounds. Monophyly of the Labridae was enforced based on previous phylogenetic reconstructions (Westneat and Alfaro 2005; Alfaro et al. 2009; Cowman and Bellwood 2011; Baliga and Law 2016). Posterior samples from three independent MCMC analyses, each with 80 million generations, sampling every 4,000 generations, were assessed for convergence and appropriate burn-in in Tracer v1.5 (Rambaut et al. 2009). Tree files were combined using LogCombiner v1.8.1 (Drummond et al. 2012) following the removal of 27.5–40% burn-in, and resampling every 16,000 states (file: PosteriorDistribution_9625TimeTrees.trees). A maximum clade credibility tree was constructed using TreeAnnotator v1.8.1 (Drummond et al. 2012) to display median ages and 95% highest posterior density (HPD) intervals (fig. S1; file: MaximumCladeCredibility_TimeTree.nex). 

We randomly sampled 1,000 time-calibrated phylogenies from the posterior distribution and pruned them to match our trait dataset (n = 89 species; file: RandomSample_1000TimeTrees_89spp.nex). Further pruning was required for the second trait correlation analysis we ran in Discrete to assess predicted correlations between polygyny and protogyny with different degrees of male size advantage (species coded as either lek-like of haremic and diandric or monandric; n = 70 species; file: RandomSample_1000TimeTrees_70spp.nex).

TableS5.xlsx

Mating system or sex allocation data are available for other labrid species. We considered only those species with both types of data avialable that met our source requirements.