Macroevolutionary patterns of sexual size dimorphism among African tree frogs (Family: Hyperoliidae)
Portik, Daniel; Blackburn, David; McGuire, Jimmy (2020), Macroevolutionary patterns of sexual size dimorphism among African tree frogs (Family: Hyperoliidae), Dryad, Dataset, https://doi.org/10.5061/dryad.t76hdr7z6
Sexual size dimorphism (SSD) is shaped by multiple selective forces that drive the evolution of sex-specific body size, resulting in male or female-biased SSD. Stronger selection on one sex can result in an allometric body-size scaling relationship consistent with Rensch’s rule or its converse. Anurans (frogs and toads) generally display female-biased SSD, but there is variation across clades and the mechanisms driving the evolution of SSD remain poorly understood. We investigated these topics in a diverse family of African treefrogs (Hyperoliidae). Hyperoliids display traits considered rare among amphibians, including sexual dichromatism and protogynous sex change. Using phylogenetic comparative methods, we tested if adult ecology, sexual dichromatism, and sex change were predictors of body size or SSD. We also tested whether hyperoliids displayed allometric interspecific body-size-scaling relationships. We found a majority of hyperoliid taxa display female-biased SSD, but that adult ecology and sexual dichromatism are poor predictors of sex-specific body size and SSD. Regardless of the groupings analyzed (partitioned by clades or traits), we found support for isometric body-size scaling. However, we found that sex change is a significant predictor of SSD variation. Species in the Hyperolius viridiflavus complex, which putatively display this trait, show a significant reduction in SSD and tend to be sexually monomorphic in size. Although protogynous sex change needs to be validated for several of these species, we tentatively propose this trait is a novel mechanism influencing anuran body size evolution. Beyond this association, additional factors that shape the evolution of anuran body size and SSD remain elusive.
We use a newly available time-calibrated phylogeny of Afrobatrachia (including Arthroleptidae, Brevicipitidae, Hemisotidae, and Hyperoliidae) produced by Portik et al. (2019) as a basis for our phylogenetic comparative methods. This phylogeny was created through a multi-step approach, which we briefly summarize here. A hyperoliid species tree of 140 lineages was constructed using 1,047 exons obtained from transcriptome-based exon capture experiment (Portik et al. 2016b), and this species tree was used to constrain hyperoliid relationships in an expanded divergence dating analysis of Afrobatrachia based on a supermatrix of one mitochondrial (16S) and five nuclear loci (FICD, KIAA2013, POMC, TYR, and RAG-1) (Portik et al. 2019). The time-calibrated phylogeny contains 173 hyperoliid lineages, which were pruned to match the species representation in our morphological data set. We note that Hyperoliidae is in a state of taxonomic flux with recent studies recommending both the synonymy of species and splitting of species complexes (Rödel et al. 2002; Rödel et al. 2003; Wollenberg et al. 2007; Rödel et al. 2009; Rödel et al. 2010; Schick et al. 2010; Conradie et al. 2012; Dehling 2012; Channing et al. 2013; Conradie et al. 2013; Greenbaum et al. 2013; Liedtke et al. 2014; Loader et al. 2015; Portik et al. 2016a; Barratt et al. 2017; Bell et al. 2017a; Conradie et al. 2018). Therefore, we adopt the naming convention of Portik et al. (2019) to label taxa included in this study, where genetically and/or geographically distinct units within large species complexes are distinguished using integers (e.g, Afrixalus dorsalis 1, Afrixalus dorsalis 2).
We included six additional lineages of the H. viridiflavus complex represented in our morphological dataset that were not present in the phylogeny of Portik et al. (2019). These include H. marmoratus verrucosus, H. m. argentovittis, H. viridiflavus ferniquei, H. v. goetzi, H. v. rubripes, and H. v. viridiflavus. We obtained relevant 16S sequences from GenBank for these six lineages, which were generated by Wieczorek et al. (2000) in their mtDNA study of the H. viridiflavus species complex. We aligned these sequences with other 16S sequences from the H. viridiflavus complex available from Portik et al. (2019) using MAFFT with the automatic alignment algorithm selection option (Katoh & Standley 2013). We reconstructed a phylogeny for the H. viridiflavus group using a ML approach in RAxML v8 (Stamatakis 2014). We subsequently added the six missing tips to the time-calibrated Afrobatrachia phylogeny based on these phylogenetic results, using the ‘bind.tip’ function in ‘phytools’ package (Revell 2012) in R (R Core Team 2015).
We obtained body size measurements in the form of snout-urostyle length (SUL) from 2,771 preserved specimens representing 123 hyperoliid lineages from the following natural history collections: Burke Museum of Natural History and Culture, California Academy of Sciences, Cornell University Museum of Vertebrates, and the Museum of Vertebrate Zoology (Supplemental Appendix 1). Measurements were made with a Mitutoyo Series 500 Digimatic Caliper (Mitutoyo U.S.A., Illinois) and recorded to the nearest 0.1 millimeters. We obtained additional body-size data from several published sources (Schiøtz 1967; Nussbaum & Wu 1995; Schiøtz 1999; Channing 2001; Rödel et al. 2002, 2003; Glaw & Vences 2007; du Preez & Carruthers 2009; Rödel et al. 2009; Harper et al. 2010; Rödel et al. 2010; Amiet 2012; Conradie et al. 2012; Channing et al. 2013; Conradie et al. 2013; Loader et al. 2015; Conradie et al. 2018) to either supplement the number of specimens included for a given species (n = 9 species) or provide data for species we could not access (n = 26 species) (Supplemental Appendix 1). The combination of empirical and published data resulted in a total of 142 lineages that have SUL data available for both sexes, of which 138 are represented in the revised phylogeny.
We classified species as sexually dichromatic or sexually monochromatic following Portik et al. (2019). The six additional species of the H. viridiflavus complex were all considered to be sexually dichromatic (Schiøtz 1999). Protogynous sex change was initially documented in two non-sister taxa in the H. viridiflavus superspecies complex (including H. viridiflavus ommatostictus and H. marmoratus taeniatus; Grafe & Linsemair 1989), and never studied again. For the purpose of our analyses, we assumed this trait occurs throughout lineages currently or formerly classified as part of the H. viridiflavus complex, including all subspecies of H. viridiflavus, H. parallelus, and H. marmoratus. For our study, this included 14 lineages. We note that the clade containing these three species groups also contains several additional species (Fig. 3), but given the uncertainty of this trait we provide a conservative estimate of its occurrence.
We calculated a sexual size dimorphism index (SSDi) to quantify SSD across species using the equation SSDi = [(larger sex / smaller sex) – 1], arbitrarily set negative if males are the larger sex and positive if females are the larger sex (Lovich & Gibbons 1992). This SSDi has been widely used, is properly scaled around zero, and has high intuitive value because positive SSDi values indicate female-biased SSD and negative values indicate male-biased SSD (Lovich & Gibbons 1992). To calculate the SSDi, we used the mean body size of each sex when possible, but for some literature records that only provided body size ranges we used the range midpoint as an alternative.
When SSDi is calculated, a species is automatically classified as displaying male or female-biased SSD, even if the difference between male and female body size is minimal. In other words, a statistical distinction between sexual size dimorphism and monomorphism is not made based on this measure. To address this, we developed a new approach to statistically classify species as dimorphic or monomorphic, which relies on intraspecific variation in body size. For a given species with body size data available for multiple males or females, we calculated an SSDi for all possible male-female pairwise comparisons, and obtained the mean value. To distinguish this metric, we refer to it as the pairwise SSDi (vs. standard SSDi). We then performed a permutation test with 10,000 bootstrap replicates to evaluate the null hypothesis that male and female body sizes come from the same population (e.g., SSDi = 0). For each replicate, we randomly shuffled the labels of males and females, calculated all possible pairwise SSDi values, and obtained the mean SSDi value. The 10,000 simulated SSDi values of the permuted data represent the estimate of the sampling distribution under the null hypothesis. We then assessed whether the empirical mean pairwise SSDi was outside of the critical values of the simulated distribution (2.5% and 97.5%, representing a 5% significance level), and calculated the P-value for the permutation test following Phipson & Smyth (2010). In this context, a P-value < 0.05 allows rejection of the null hypothesis and classifies a species as displaying sexual size dimorphism. Alternatively, a nonsignificant P-value indicates there is not sufficient evidence to support the alternative hypothesis (e.g., male and female body sizes come from different populations, SSDi ≠ 0). Here, we classify these species as displaying sexual size monomorphism, but acknowledge the failure to reject the null hypothesis could result from a true lack of body size difference between the sexes, or from artifacts such as insufficient sample sizes or high variation. We performed theses analyses for all species with available data, which included 120 hyperoliid species. Using this method, the mean pairwise SSDi calculated for a given species takes intraspecific variation into account. Consequently, we were interested in determining how similar this value was to the standard SSDi calculated from an average male and female body size. To determine this, we simply calculated the absolute difference between the two SSDi measures for each of the 120 species. To perform all calculations of SSDi, pairwise SSDi, and permutation tests, we created a Python program (SSDi-Calculator.py) which is freely available from: https://github.com/dportik/SSDi-Calculator.
Phylogenetic Comparative Analyses
To visualize the distribution of sex-specific body sizes and SSDi on the hyperoliid phylogeny, we performed maximum likelihood (ML) ancestral state reconstruction for continuous traits using the ‘contMap’ function of the ‘phytools’ package (Revell 2012) in R (R Core Team 2015). To visualize body sizes, we used log-transformed values of SUL.
We tested for phylogenetic signal in sex-specific body size and SSDi to assess if the evolution of these traits deviates from expectations based on Brownian motion. We quantified the amount of phylogenetic signal using Blomberg’s K-statistic (Blomberg et al. 2003). For this metric, a K value close to 1 indicates trait evolution occurs according to a Brownian motion model. A K value of less than 1 indicates relatives are less similar than expected from Brownian motion, and a K value greater than 1 indicate relatives are more similar than expected. We used 1000 permutations in a randomization test to examine if the variance in the empirical data was significantly different than variance in randomized data sets with no phylogenetic signal. As a complement we also examined phylogenetic signal using Pagel’s λ, which is a scaling parameter used to transform the phylogeny such that it ensures the best fit to a Brownian motion model (Pagel 1999; Freckleton et al. 2002). If Pagel’s λ is close to 1, the trait evolves according to Brownian motion, whereas values approaching 0 indicate the trait has evolved independently of the phylogeny (e.g., a lack of phylogenetic signal). We evaluated if the fitted value of λ was significantly different from λ = 0 using a likelihood-ratio test. We conducted both types of tests for phylogenetic signal using the phylosig function in the ‘phytools’ package in R, using SSDi and log transformed values of SUL.
We performed phylogenetic generalized least squares (PGLS) regression to examine potential correlations between several predictor variables (ecology, sexual dichromatism, and protogynous sex change) and several response variables (male SUL, female SUL, and SSDi). PGLS incorporates phylogenetic nonindependence into generalized linear models in the form of a phylogenetic variance-covariance matrix (Freckleton et al. 2002). We determined the most appropriate correlation structure for the residuals for each regression by comparing the fit of Brownian motion, Brownian motion plus the λ scaling parameter, and Ornstein-Uhlenbeck (OU) models using AIC scores.
Testing Rensch’s Rule
For systems with female-biased SSD, Rensch’s rule states that the magnitude of SSD decreases with increases in body size, and the converse of Rensch’s rule states that the degree of SSD should increase with mean body size (Fig. 1). These allometric scaling relationships are generally tested by determining if the regression slope between log-transformed male and female body size differs significantly from a value of one (Abouheif and Fairbairn 1997; Fairbairn 1997; Smith 1999). Model I regressions (including ordinary least squares, OLS) assume the values for the x-axis are known without error, and the OLS equation is not symmetrical, producing different predictions depending on which sex is defined as the independent variable (Warton et al. 2006; Smith 2009). Reduced major axis regression (RMA) is a Model II regression that is symmetric, meaning that a single line defines the bivariate relationship regardless of which variable is X and which is Y. Furthermore, it assumes that variables on both axes are measured with error (Warton et al. 2006; Smith 2009). We therefore chose to use RMA because our estimates of the body size for each sex contain error, and because the assignment of male or female size as X or Y is arbitrary. When female size is assigned to the x-axis and male size is assigned to the y-axis, allometry conforming to Rensch’s rule is indicated by a slope greater than one, whereas the converse of Rensch’s rule is supported when the slope is less than one (Fig. 1). When the slope is not significantly different from one, male and female body size scale isometrically (Fig. 1).
To account for shared evolutionary history, we performed a phylogenetic reduced major axis (pRMA) regression (Ives et al. 2007) using the phyl.RMA function in ‘phytools’. Hypothesis testing for a slope significantly different than 1 was calculated following Ives et al. (2007) and implemented using the phyl.RMA function. We performed pRMA regressions for multiple taxonomic levels (family, subfamily, clade), and for groups displaying particular traits (sexual dichromatism, protogynous sex change). Two subfamilies are currently recognized within Hyperoliidae: Kassininae (~26 species), which is composed of both arboreal and terrestrial species, and Hyperoliinae (~206 species), which contains exclusively arboreal species (Fig. 3). Variation in ecological and reproductive traits is greatest between these subfamilies, but important variation in arboreal species is also captured in two major clades within Hyperoliinae: one clade contains species of Afrixalus, Heterixalus, and Tachycnemis (~43 species, hereafter referred to as Clade A), and the other contains the hyperdiverse genus Hyperolius and two species-poor genera (Cryptothylax, Morerella) (~156 species, hereafter referred to as Clade B) (Fig. 3). To disentangle any potential clade-specific patterns, we performed pRMA regressions for Hyperoliidae, Kassininae, Hyperoliinae Clade A, and Hyperoliinae Clade B. To investigate body-size scaling as it relates to sexual dichromatism, we performed separate pRMA regressions for dichromatic and monochromatic species occurring in Hyperoliinae. To investigate scaling as it relates to protogynous sex change, we also performed a pRMA regression for species in the H. viridiflavus complex. To visualize these data, we created bivariate scatterplots from the natural-log transformed SUL of males and females for a given comparison (with female body sizes plotted on the x-axis) and plotted the relevant pRMA regression line.
National Science Foundation, Award: DEB: 1311006