The spectacular diversity of insect male genitalia, and their relative insensitivity to the environment, have long puzzled evolutionary biologists and taxonomists. We asked whether the unusual evolvability of male genitalia could be associated with low morphological integration of genitalic traits, by comparison with male somatic traits and female traits. We also asked whether this pattern was robust to variation in resource availability during development, which affects adult condition. To address these questions, we manipulated larval diet quality in a split-brood design and compared levels of integration of male and female genitalic and somatic traits in the neriid fly, Telostylinus angusticollis. We found that male genitalic traits were substantially less integrated than male somatic traits, and less integrated than female genitalic traits. Female genitalic traits were also less integrated than female somatic traits, but the difference was less pronounced than in males. However, integration of male genitalic traits was negatively condition-dependent, with high-condition males exhibiting lower trait integration than low-condition males. Finally, genitalic traits exhibited lower larval diet family interactions than somatic traits. These results could help explain the unusually high evolvability of male genitalic traits in insects.

We utilised a full-sib, split-family breeding design where randomly chosen individuals from these populations were paired to create 17 mating pairs at 15 ± 2 days old. Each pair was allowed 48 hours to mate within 120 mL cages provided with a nutrient-rich oviposition medium and access to sugar, yeast and water ad libitum. Following the 48-hour period, from each mating pair we transferred 20 eggs to a nutrient-poor larval diet and 20 eggs to a nutrient-rich larval diet, also based on [33]. Upon emergence, virgin adult offspring were allowed 24 hours for their exoskeletons to sclerotize fully and then frozen at -80°C for dissection and morphological measurements. We quantified six genitalic and 12 somatic traits on each of 93 males, and four genitalic and 11 somatic traits on each of 96 females. All traits were quantified by measuring the lengths of the structures (see electronic supplementary material), except for testes (TE), for which we measured area in mm². For both sexes we used thorax length as an index of body size [34]. To minimise the loss of samples for multivariate analyses, missing trait values were replaced with the mean value for the family × larval diet × sex combination (where > 3 individuals were available for that treatment combination).

All analyses were carried out using R 3.5.3 [35]. For each set of traits (genitalic and somatic) and group combinations (sex and larval diet), we estimated morphological integration as the relative standard deviation of eigenvalues, SD_rel (λ) [36]. The higher the value of SD_rel (λ), the more variance is explained by the first few principal components of the trait matrix, and therefore the higher the integration. Integration was estimated from principal components analysis (PCA) performed separately on the correlation matrix for each sex × larval diet × trait type combination (Figures S3, S4). Standard errors for integration values were obtained by a Jackknife (n-1 traits) procedure. As a measure of environmental effects (i.e., larval diet quality), we computed marginal effect sizes, which indicate the variance explained by fixed effects [37]. To estimate the maximum genotype × environment interaction (G × E), we estimated the family × larval diet interaction. We computed conditional effect sizes, which reflect the variance explained by both fixed and random effects [37], and estimated the family x larval diet interaction by comparing the magnitude of the marginal and conditional effect sizes (see electronic supplementary material for details). Median parameter values were compared between treatment groups using nonparametric Kruskal-Wallis or Wilcoxon tests, with trait as the unit of replication.

Male and female trait values are transformed using log10. Repeatabilities were obtained by re-positioning each sample/trait and taking multiple images of the same sample. There are missing values in datasets because some samples were either damaged or major outliers (>3 sd) in the data and removed.