The macroevolution of size and complexity in insect male genitalia

Morphological measurements

Limnebius includes ca. 150 species with an almost cosmopolitan distribution, all of them aquatic, living in all types of continental waters with the only exception of saline habitats (Perkins, 1980; Jäch, 1993; Hansen, 1998; Rudoy, Beutel & Ribera, 2016). There is very few information on the ecology and behaviour of the species of Limnebius, including their sexual habits. In a recent work, Limnebius has been shown to be divided in two sister lineages with an estimated Oligocene origin, the subgenera Bilimneus and Limnebius s.str., with ca. 60 and 90 described species respectively (Rudoy, Beutel & Ribera, 2016).

We characterised the size and complexity of the male genitalia of all species included in the molecular phylogeny (see below and Table S1). The structure of the female genitalia of Limnebius is poorly known, mostly due to the lack of sclerotized parts (Rudoy, Beutel & Ribera, 2016), and was not studied here. We measured body length of adults (males, lm and females, lf) as the sum of the individual maximum lengths of pronotum and elytra, as the different position of the articulation between the two could alter the total length when measured together. Similarly, the head was not measured, as in many specimens it was partly concealed below the pronotum. Measures were obtained with stereoscope microscopes equipped with an ocular micrometer.

Male genitalia (aedeagi) were dissected and mounted on transparent labels with dimethyl hydantoin formaldehyde (DMHF). For size measurements we used as a single value the average of each measure in all studied specimens of the same species (Table S1). For the measures of complexity a single specimen was used, as species show in general a very constant shape of the aedeagus (Jäch, 1993; Rudoy, Beutel & Ribera, 2016), and there is no evidence for the widespread existence of cryptic species within the genus (unpublished data). Male genitalia were always photographed in the same standard positions, orientated according to the foramen in ventral and lateral views. We measured the maximum length of the male genitalia (lg) on digital images obtained in the same standard position, orientated in ventral view as determined by the foramen (Rudoy, Beutel & Ribera, 2016). We did not include setae or apical membranous structures, but included appendages when they were longer than the median lobe (as in e.g., some species of the L. nitidus group, Rudoy, Beutel & Ribera, 2016). Measurements were directly obtained from the digital images using ImageJ v.1.49 (US National Institutes of Health: http://imagej.nih.gov/ij/) (Fig. S1). We estimated experimental error by measuring the same specimen of three species on three different sessions, using two sets of images.

We used two different measures to characterise the complexity of the aedeagus:

(i) Perimeter (per) of the aedeagus in ventral view, including the median lobe and the main appendages. We obtained an outline of the genitalia from digital images using ImageJ, with the “free hand outline” option. The total perimeter was the sum of the values of the different parts of the genitalia (median lobe and left parameter, plus main appendages if present, see Rudoy, Beutel & Ribera, 2016). We standardised the values by dividing the perimeter by the length of the aedeagus, to obtain a measure of complexity by unit of length (Fig. S1 and Table S1).

(ii) Fractal dimension (fd). We estimated the fractal dimension of the outline of the aedeagus in ventral view on images of standard size (2100x2100 pixels, 2000 pixels from base to apex of the aedeagus) with the software Fractal Dimension Estimator (http://www.fractal-lab.org/index.html). This software estimates the Minkowski fractal dimension of bidimensional images using the box-counting method (Falconer, 1990). The software converts the image to binary data, selects the scaling window of the box, and counts how many boxes are necessary to cover the image. The absolute value of the slope of a log–log graph of the scale with the number of boxes is the fractal dimension of the image (Fig. S1, Table S1).

Phylogenetic analyses

We obtained a phylogeny of Limnebius using molecular data obtained from Rudoy, Beutel & Ribera (2016), which included 71 of the ca. 150 known species of the genus. The taxon sampling was denser for the Palaearctic lineages in subgenus Limnebius s.str., including the full range of body sizes and aedeagus structural variation (Rudoy, Beutel & Ribera, 2016). Groups with the lowest number of sampled species show in general a more homogeneous and simpler aedeagus (L. mundus and L. piceus groups and subgenus Bilimneus). We used as outgroup and to root the tree the genus Laeliaena, considered to be sister to Limnebius based on multiple morphological synapomorphies (Hansen, 1991; Jäch, 1995; Perkins, 1997; Beutel, Anton & Jäch, 2003).

We reconstructed a phylogeny of the genus with Bayesian analyses in BEAST 1.8 (Drummond et al., 2012) using a combined data matrix with three partitions, (i) mitochondrial protein coding genes (two cox1 fragments plus nad1); (ii) mitochondrial ribosomal genes (rrnL plus trnL) and (iii) nuclear ribosomal genes (SSU plus LSU) (Table S1). We used a Yule speciation process as the tree prior and an uncorrelated relaxed clock. Analyses were run for 100 MY generations, ensuring that the number of generations after convergence were sufficient as assessed with Tracer v1.6 (Drummond et al., 2012) and after removal of a conservative 10% burn-in fraction.

To calibrate the phylogenetic tree we used the rates estimated in Cieslak, Fresneda & Ribera (2014) for a related group (familiy Leiodidae, within the same superfamily Staphylinoidea, Beutel & Leschen, 2005) and the same gene combination based on the tectonic separation of the Sardinian plate 33 Ma. We set as prior average mean rate a normal distribution with average 0.015 substitutions/site/Myr for the mitochondrial protein genes, 0.006 for the mitochondrial ribosomal genes, and 0.004 for the nuclear ribosomal genes, all with a standard deviation of 0.001. It must be noted that for our objectives only relative rates are needed. An absolute calibration would only be necessary to obtain absolute estimates of character change, which is not our main objective and does not affect our conclusions.

We reconstructed the ancestral values of the morphological variables using the values of the terminals (extant species) in BEAST 1.8. The reconstruction of ancestral characters should ideally include a combination of extant and fossil data (e.g., Webster & Purvis, 2002; Slater, Harmon & Alfaro, 2012; Bokma et al., 2015). However, as there are no known fossils of species of Limnebius we can only use extant species. The use of extant species without fossil data has resulted in valuable contributions to the study of character evolution (e.g., Cooper & Purvis, 2010; Baker et al., 2015), and in some cases it has been reported that the introduction of fossil data did not alter substantially the conclusions obtained with extant species only (e.g., Puttick & Thomas, 2015).

We implemented a Brownian movement model of evolution (BM), a null model of homogeneous evolution in which variation accumulates proportionally with time, with incremental changes drawn from a random distribution with zero mean and finite constant variance (Hunt & Rabosky, 2014; Adams, 2014). The reconstruction of ancestral values using a BM model of evolution can be biased toward average or intermediate values (Pagel, 1999; Finarelli & Goswami, 2013), which together with the lack of fossil data can result in an underestimation of the rates of evolution of some characters. Due to these limitations our reconstruction needs to be understood as the simplest null model explaining the evolutionary change in the studied characters. It must be noted though that some of our conclusions do not depend on these reconstructions, and for others we compare two groups (Bilimneus and Limnebius) likely affected by the same biases. In any case, all analyses were repeated using only the data of the terminal branches, with a decreased power due to their lower number but less prone to reconstruction biases. Terminal branches are truncated with respect to internal branches, as they do not reflect a full interval between speciation events but the time elapsed since the last cladogenetic split. This does not introduce any bias in our measures when change is anagenetic (i.e., proportional to time), but will overestimate evolutionary rates if change occurs mostly at speciation (i.e., according to punctuated models of evolution; Gingerich, 2009; Hunt & Rabosky, 2014).

Statistical analyses

We studied the evolution of the morphological characters trough the full evolutionary path of species (i.e., from root to tips) and in the individual branches, using phylogenetic ancestor-descendant comparisons (PAD; Baker et al., 2015). We first deleted the terminals with species with unknown males, as well as the duplicated specimen of L. cordobanus and the outgroup (Laeliaena) (Table S1). We then estimated the phylogenetic signal of the morphological variables in the whole tree using the K metric proposed by Blomberg, Garland & Ives (2003), which tests whether the topology and branch lengths of a given tree better fits a set of tip data compared with the fit obtained when the data have been randomly permuted. The higher the K statistic, the more phylogenetic signal in a trait. K values of 1 correspond to a BM model, which implies some degree of phylogenetic signal. K values closer to zero correspond to a random or convergent pattern of evolution, while K values greater than 1 indicate strong phylogenetic signal. We used the R package ‘Picante’ (Kembel et al., 2010) to compute K and the significance test. To test the possible decrease of power or the randomization test due to the low number of species in Bilimneus (15, vs. 50 in Limnebius s.str.) we randomly pruned the tree of Limnebius s.str. to 15 species in 1,000 replicas and computed the K values for all of them. We also measured the correlation between the variables across the whole tree with a regression of phylogenetic independent contrasts with the PDAP package in MESQUITE (Midford, Garland & Maddison, 2011).

For the analyses of the individual branches of the tree (PADs) we compiled the reconstructed values in BEAST of all variables for the initial and final node of each individual branch, together with their length. We measured three values for each of the individual branches (including terminals): (i) amount of phenotypic change, equal to the arithmetic difference between the final and initial values of the branch; (ii) absolute amount of phenotypic change, equal to the absolute value of the amount of phenotypic change; (iii) phenotypic change measured in darwins (Haldane, 1949), computed as the absolute value of the natural logarithm of the ratio between the final and initial values divided by the length of the branch in million years (Myr) (Table S1). The use of the natural logarithm standardises the change so it is proportional and directly comparable among species with different sizes (Haldane, 1949; Gingerich, 2009). Darwins are known to be strongly dependent on the length over which are measured (Gingerich, 2009; Hunt, 2012; Hunt & Rabosky, 2014), so to obtain the correlation with time of the different variables we used instead the absolute amount of phenotypic change, which showed a more linear relationship with branch length for most of the variables (see Results). As individual branches are in principle independent from each other we analysed these variables with standard statistical procedures (see e.g., Baker et al., 2015). We used a Bonferroni correction for multiple tests and considered as significant a p < 0.01 level (corresponding to five morphological variables), although we also report some results at the standard p < 0.05 level in some of the analyses. In the comparisons between the two subgenera of Limnebius (Bilimneus and Limnebius s.str.) we did not include the respective stem branches.

Phylogenetic signal and correlation between complexity and size

The reconstructed ancestral values of Bilimneus were in the upper range of the values of the extant species, especially for measures of complexity (fd and per). The reconstructed ancestral values of the whole genus were, on the contrary, larger than any extant species of Bilimneus, with the only exception of the length of the genitalia (lg) (Table 1). For the subgenus Limnebius s.str. all reconstructed values, both of the genus and the subgenus, were within the range of the extant species, and close to their centre of distribution (Table 1; Tables S1 and S2).

The phylogenetic signal (K) was very significant for all variables both in the genus Limnebius as a whole and the subgenus Limnebius s.str., with values close to 1 (i.e., a BM model) for the measures of complexity (per and fd), clearly over 1 (i.e., with a stronger phylogenetic signal) for the size of males and females (lf and lm), but with values lower than 1 (i.e., closer to a global random distribution, although still significant) for the length of the genitalia (lg) (Table 2).In contrast, none of the variables had a significant K in subgenus Bilimneus (Table 2). More than 99% of the 1,000 pruned replicas of Limnebius s.str. maintained significant values of K for all variables at a p < 0.05 level, and more than 90% at p < 0.01, with the only exception of lg, which become not significant (Table 2).

All pairwise correlations between the morphological variables were significant at a p < 0.01 level as measured with PDAP both in Limnebius as a whole and in subgenus Limnebius s.str., with the exception of the correlation between lg and the measures of complexity (per and fd) (Table 3). For Bilimneus, with a lower number of contrasts (14 vs 49 in Limnebius s.str.), results were similar except that the length of the female (lf) was not significantly correlated with any variable, and the length of the aedeagus (lg) was significantly correlated with its fractal dimension (Table 3).

Rates of evolution

There were no significant differences between terminal and all branches of the tree for any of the measured variables (including branch length) at a p < 0.01 level, both in average value (2 tailed t-Student with unequal variances when differences were significant) or variance (Fisher’s F-test) (Table 4). At a p < 0.05 only for the measure of darwins the variances of per and lf were larger in the terminals than in all branches, and the variance of lg smaller (Table 4). Some of the largest increases or decreases of absolute change were in terminal branches (e.g., L. nitiduloides and L. cordobanus respectively, branches 94 and 72 in Fig. S2).

The variable with the highest values of phenotypic change and the highest variance was lg, both measured in absolute phenotypic change or in darwins, and both in Bilimneus and Limnebius s.str. (Table 4). The lowest values of change were in fd, also in both subgenera. Males were in general more labile than females, with larger relative (but not absolute) differences in Bilimneus than in Limnebius s.str. (Table 4).

All comparisons of absolute change and darwins between Bilimneus and Limnebius s.str. were significantly different (always lower in Bilimneus) both for the variance and the average, with the only exception of the average lg, never significantly different between them regardless of how it was measured (Table 4; Fig. S3). When the sign of the change was considered variances were still significantly different, but only the average fd showed differences (positive in Limnebius s.str. and negative in Bilimneus), although only at a p < 0.05 level (Table 4; Fig. S3).

Branches of Limnebius s.str. were on average significantly shorter and with a lower variance than that of Bilimneus (Table 4; Fig. S3). In the whole Limnebius and the subgenus Limnebius s.str. all the measures of absolute phenotypic change were significantly correlated with the length of the branch, with lg having the lowest correlation (Fig. 2A; Table 5). On the contrary, for Bilimneus only lg and per had significant correlations. When only the terminal branches were used the correlations between the measures of absolute phenotypic change and branch lengths were still significant at a p < 0.01 level for Limnebius s.str. with the only exception of lg, but none of them was significant for Bilimneus (Fig. 2B and Table 5).

Distribution of the reconstructed phenotypic change

The amount of change in the individual branches for all variables was in general not significantly different from a normal distribution, as measured with a non-parametric Kolmogorov–Smirnov test. Only the variables of complexity when all the branches in the whole genus were considered could be said not to follow a normal distribution at a p < 0.01 level (Table 6). When measured in darwins (and taking into account the sign of the branch), the amount of change of the individual variables was in general less normally distributed than when the change was measured with independence of the time, although at a p < 0.01 level only per and lf in the whole genus, and per in Limnebius s.str., were significantly non-normal (Table 6).

Despite the general normal distribution of the amount of change, in Bilimneus there was an excess of negative branches (i.e., a change to lower values) for all variables except for lg. On the contrary, in Limnebius s.str. there were significantly more positive than negative changes for the variables of complexity (Table 7). In both cases differences were mostly significant only at a p < 0.05 level, and in Limnebius s.str. only when all branches of the tree were considered. The average amount of phenotypic change per branch was negative for all measures in Bilimneus and positive in Limnebius s.str., but in all cases the 95% confidence intervals included the zero (data not shown).

For all measurements the general tendency was for a branch to have the same type of change than that of the preceding. In Bilimneus this meant that the most frequent combination was two branches with negative change, while in Limnebius s.str. two with positive change, although lg had the lowest differences between the frequency of two consecutive positive or negative branches in both subgenera (Table 7). The terminal branches showed the same tendency except for the variables measuring length in Limnebius s.str., for which the most frequent combination was two negative branches (Table 7).