DOI: 10.5061/dryad.t76hdr7z6
title: Macroevolutionary patterns of sexual size dimorphism among African tree frogs (Family: Hyperoliidae)

description: 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.
description: <i>Molecular Phylogeny</i> 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
 (<i>16S</i>) and five nuclear loci (<i>FICD, KIAA2013,
 POMC, TYR, </i>and<i> RAG-1</i>) (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, <i>Afrixalus dorsalis</i> 1<i>, Afrixalus dorsalis
 </i>2).             We included six additional
 lineages of the <i>H. viridiflavus</i> complex represented in
 our morphological dataset that were not present in the phylogeny of Portik
 et al. (2019). These include <i>H. marmoratus verrucosus, H. m.
 argentovittis, H. viridiflavus ferniquei, H. v. goetzi, H. v. rubripes,
 </i>and <i>H. v. viridiflavus</i><i>. </i>We
 obtained relevant <i>16S</i> sequences from GenBank for these
 six lineages, which were generated by Wieczorek et al. (2000) in their
 mtDNA study of the <i>H. viridiflavus</i> species complex. We
 aligned these sequences with other <i>16S</i> sequences from
 the <i>H. viridiflavus</i> 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 <i>H. viridiflavus </i>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). <i>Trait
 Data</i> 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 <i>H. viridiflavus</i> complex
 were all considered to be sexually dichromatic (Schiøtz 1999). Protogynous
 sex change was initially documented in two non-sister taxa in the
 <i>H. viridiflavus</i> superspecies complex (including
 <i>H. viridiflavus ommatostictus </i>and <i>H.
 marmoratus taeniatus</i>; 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
 <i>H. viridiflavus</i> complex, including all subspecies of
 <i>H. viridiflavus</i>, <i>H. parallelus</i>, and
 <i>H. marmoratus</i>. 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
 <i>P</i>-value for the permutation test following Phipson
 & Smyth (2010). In this context, a <i>P-</i>value
 &lt; 0.05 allows rejection of the null hypothesis and classifies a
 species as displaying sexual size dimorphism. Alternatively, a
 nonsignificant <i>P</i>-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 (<i>SSDi-Calculator.py</i>) which is freely
 available from: https://github.com/dportik/SSDi-Calculator.
 <i>Phylogenetic Comparative Analyses</i>
 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.   <i>Testing
 Rensch’s Rule</i> 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
 <i>Afrixalus, Heterixalus, </i>and <i>Tachycnemis
 </i>(~43 species, hereafter referred to as Clade A), and the other
 contains the hyperdiverse genus <i>Hyperolius</i> and two
 species-poor genera (<i>Cryptothylax, Morerella</i>) (~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 <i>H. viridiflavus</i> 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.