Sexual size dimorphism and male reproductive traits vary across populations of a tropical rainforest dung beetle species (Onthophagus babirussa)

Abstract Sexual size dimorphism (SSD) arises when natural selection and sexual selection act differently on males and females. Male‐biased SSD is rarer in insects and usually indicates strong sexual selection pressure on male body size in a species. Patterns of SSD can also vary between populations of species that are exposed to different environmental conditions, such as differing resource availability and diversity. Here, we investigate intraspecific variation in SSD as well as relative investment in precopulatory (horn length) and postcopulatory traits (sperm length and testes weight) in a tropical rainforest dung beetle Onthophagus babirussa across Singapore and Peninsular Malaysia. Overall, three out of four populations displayed significant male‐biased SSD, and SSD was greater in populations with smaller overall body size. Average male body size was similar across all populations while female body size was significantly smaller in Singapore, suggesting that the pronounced SSD may also be due to stronger sexual selection on male body size in Singapore populations. All populations showed significant investment in horns as a weapon likely used in male‐male competition, while postcopulatory traits showed no clear scaling relationship with body size, suggesting a higher priority on precopulatory sexual traits in the mating system of this species.


| INTRODUC TI ON
Sexual selection is defined as selection on heritable traits that vary between individuals within a population that influence reproductive success and fitness (Andersson, 1994). When individuals within a population have differential reproductive success (Panhuis et al., 2001), this can occur prior to copulation (precopulation), when males compete for access to females, leading to evolution of sexual dimorphism in size and secondary sexual traits such as ornaments and weapons (Simmons & García-González, 2008). Sexual selection can also occur postcopulation, for example, in the form of cryptic female choice, where females can influence the success rate of insemination by males and/or via sperm competition, where sperm from different males compete to fertilize the ova (Birkhead & Pizzari, 2002).
One of the most common traits that is subject to sexual selection is body size. Sexual size dimorphism (SSD) arises when the effects of natural selection and sexual selection act differently on males and females (Blanckenhorn, 2005). In most invertebrates, such as insects, species often display female-biased SSD, where females are larger due to strong fecundity selection (Blanckenhorn, 2005;Esperk et al., 2007;Rudoy & Ribera, 2017;Stillwell et al., 2010).
Larger male body size is usually a derived trait in most insect lineages and an evolutionary reversal of the ancestral state of femalebiased SSD (Blanckenhorn et al., 2004;Blanckenhorn et al., 2007).
Most studies on the evolution of male-biased SSD in insects focus on the effects of intraspecific factors on SSD, such as male-male competition and runaway selection of female-preferred traits associated with body size (Burkhardt & de la Motte, 1988;Fairbairn & Preziosi, 1994;Pomfret & Knell, 2006;Simmons & Tomkins, 1996;Wilkinson & Reillo, 1994). Fewer studies consider sexual selection in relation to broader ecology, such as external biotic factors. One example would be Beckers et al. (2015) that explored the effect of differential resource competition on divergence in life history traits in separate populations of the dung beetle Onthophagus taurus, finding effects of developmental plasticity, parental effects, and genetic background on different traits. In this research, we investigate the differences in SSD and pre-and postcopulatory traits in in situ populations of a dung beetle species that differ in resource diversity and availability.

Species belonging to the dung beetle genus, Onthophagus
Latreille, 1802, (i.e. the most species rich genus in the animal kingdom), have been gaining increased interest as models in evolutionary research. Recent studies show that their morphology and genetic variation can be influenced by sexual selection, parental investment, and environmental variation via a multitude of complex mechanisms (Dury et al., 2020;Hu et al., 2020;Schwab et al., 2019;Snell-Rood et al., 2016). They are particularly popular in sexual selection research because many species display strong sexual dimorphisms (Parzer & Moczek, 2008). Males often possess horns, a precopulatory sexual trait, on the head and/or thorax, which are used in defending breeding tunnels occupied by females (Garcia-Gonzalez & Simmons, 2011;Kijimoto et al., 2009;Simmons & García-González, 2008). Some species exhibit tradeoffs between male horn length and investment in postcopulatory traits such as testes and sperm (Moczek & Nijhout, 2004;Reynolds & Byrne, 2013). Studies in Onthophagus have shown alternative mating strategies where smaller males prioritize investing more in testes size and sperm production over horn investment (Simmons & Emlen, 2006;Simmons & García-González, 2008;Simmons et al., 2007). Sperm length has been shown to be under extreme selection in other insect groups such as in Drosophila flies where long sperm are better able to displace sperm from competing males (Lüpold et al., 2016;Snook & Karr, 1998), while shorter sperm has been found to confer higher fertilization success in dung beetles (García-González & Simmons, 2007). These pre-and postcopulatory phenotypes are determined during larval development and affected by the environment and maternal investment such as food provisioning (Emlen, 1994;Emlen, 1997a;Moczek, 1998;Silva et al., 2016).
Body size of specimens from Peninsular Malaysia appeared similar between the sexes, contrary to specimens collected from Singapore, despite relatively close proximity (~316 km). Intraspecific differences in SSD between separate populations have been observed in other species (Cox & Calsbeek, 2010;Liao et al., 2015;Piross et al., 2019;Rossi & Haga, 2019;Teder & Tammaru, 2005), including a complete reversal of SSD in the dung fly Sepsis punctum (Puniamoorthy et al., 2012), and these are usually due to differences in sexual selection pressures acting on each population. Differences in sexual selection pressure can in turn be influenced by external factors such as resource availability (Forsgren et al., 1996;Ghislandi et al., 2018).
In this study, we investigate the variation in SSD and relative investments in pre-and postcopulatory traits within and between four separate populations of Onthophagus babirussa from Singapore and Peninsular Malaysia (henceforth, SG and MY, respectively). The precopulatory trait examined in this study was male horn length, while testes weight and sperm length were measured as postcopulatory traits. Static allometries were calculated to estimate relative investment in the traits as a function of body size, following standard protocol (Eberhard et al., 2018;Knell, 2009). Resource availability differs between SG and MY since mammal diversity is much higher in the latter. We hypothesize that since dung resources are scarcer and less diverse in SG, competition between males over monopoly of access to dung and females would be higher, and thus male-biased SSD would be more pronounced in populations from SG than from MY. In line with this, we predict that pre-copulatory selection acting on SG population is likely stronger than post-copulatory selection; since male horns are important for male-male combat and mate acquisition (Beckers et al., 2017;Moczek & Emlen, 2000;Simmons & Ridsdill-Smith, 2011), we hypothesize a greater relative investment in horn length rather than in testes size and/or sperm length. 2.1.2 | Sampling and sorting protocol Dung beetle sampling was conducted using baited pitfall traps and baited funnel pitfall traps with human dung as the bait because it is widely accepted to be the best bait to attract a wide variety dung beetles (Howden & Nealis, 1975;Kudavidanage et al., 2012;Larsen & Forsyth, 2005). Exact details of trap materials and construction are appended (Appendix 1: Figure A1). Traps were retrieved after 24-48 h, and captured beetles were brought back to the laboratory for morphological identification and sorting using an Olympus . For these specimens, the right mid femur was dissected into 7 μl of QuickExtract solution, and the DNA was extracted by following the manufacturer's protocol (Lucigen, 2018). Then, 313 bp fragments of the COI gene were amplified via PCR (see Appendix 2 for detailed protocol), sent for next-generation sequencing (NGS) and used for DNA barcoding.
Sequence analysis was then conducted with reference to the analysis pipeline detailed by (Meier et al., 2016), and a well-established 3% threshold for uncorrected pairwise distances was used to delimit different species (Hebert et al., 2003;Meiklejohn et al., 2011;Srivathsan & Meier, 2012). All specimens examined in this study fell F I G U R E 1 Map of sampling sites located in Singapore and Malaysia. Colours represent the different sites that were treated as separate populations for analyses. within the same molecular cluster under this 3% threshold, and a cluster fusion diagram with representatives from each population is appended in Appendix 2 ( Figure A4), along with the full protocol for morphological and molecular sorting. The molecular barcodes were congruent with our morphological sorting and general consensus with the geographical sampling.

| Precopulatory trait measurements
To investigate the sexual size dimorphism in the four populations of O. babirussa, maximum pronotum width (Figure 2) of males and females was measured as a proxy for body size with the eyepiece reticle on the Olympus SZX10 microscope. This is widely used as a proxy for body size because the pronotum width does not change in adulthood and has been found to be the most appropriate measure for body size in dung beetles (Emlen, 1997a(Emlen, , 1997bKnapp & Knappová, 2013).
Horn lengths of male O. babirussa ( Figure 2) were measured to document variation in this precopulatory trait. Images were taken of the anterior habitus. Heads of the beetles were separated and suspended with Durex KY Jelly, with horns aligned parallel to the lens of the camera. Images were captured using the EOS 800D and 6D camera body with the Canon MP-E 65 mm f/2.8 1-5× lens at 5× optical zoom. The camera was suspended on the Dun, Inc. P-51, and the Camlift controller V2.9.3.0 software was used to take multiple images at different heights for focus stacking. EOS Utility Launcher software was used to access the images and stack them using the Zerene Stacker V. 1.04. software. Stacked images were imported to Adobe Photoshop CS5 V. 12.0 ×64, and a 1 mm scale bar was added to each image. Next, processed images were imported to ImageJ V.
1.51, and the horns were measured from the tip to the bottom of the outer edge of each horn, following previous studies (Moczek & Emlen, 1999  on a frosted slide. Then, sperms were teased out from the vesicles using an insect pin. Slides were dried in the oven, and sperms were fixed onto the slides with a solution of three parts methanol and one part acetic acid for 2 min. Next, the slides were washed in 1× PBS for 1 min, and the sperms were stained for 5 min in the dark with 4′,6-diamidino-2-phenylindolev (DAPI), which binds to DNA to form a fluorescent complex to allow for visualization of sperm heads under a fluorescent microscope. Following that, the slides were washed in 1× PBS and placed in the dark to dry. When the slides were dried completely, one to two drops of glycerol were added on the stained regions, coverslips were placed, and the edges were sealed with clear nail polish and left to dry in the dark. The sperms were visualized using an Olympus BX50 fluorescence microscope and measured using μManager and ImageJ V. 1.51 software.

| Statistical analyses
Box plots of average pronotum width were constructed with confidence intervals using the R packages ggplot2 (Wickham, 2016), dplyr (Wickham et al., 2020), and plotrix (Lemon, 2006) and tested for significance in body size difference between the sexes within each population using ANOVA, checking the residuals for normality after. To test if SSD varied between populations, we ran linear models testing for significant sex by location interaction. Post-hoc analyses using Dunn test were also conducted to determine which populations differ from the other for male and female body size. In addition, the sexual dimorphism index (SDI) was calculated for each population following the formulation by Lovich and Gibbons (1990), where the mean size of the larger sex is divided by the mean size of the smaller sex. A negative sign is arbitrarily added to the SDI as the males are larger (Lovich & Gibbons, 1990).
To determine whether populations differed with respect to relative investments in precopulatory and postcopulatory traits, the static allometries were calculated by first constructing log-log scatterplots of trait size against pronotum width. As the log-log scatter plot of horn length against pronotum width displayed a clear nonlinear relationship, we followed the recommendations by Knell (2009) and Parrett et al. (2021) and fitted (1) linear model, (2) quadratic model, (3) cubic model, and (4) breakpoint model using the R package segmented (Muggeo, 2008) to the pooled data with all four populations to characterize the trait size-body size relationship Figure 3.
Model selection was then conducted with the Akaike information criterion (AIC). The breakpoint model had the lowest AIC score for horn length (Table 1), indicating that this model is the best model for explaining the relationship between the variables (Knell, 2009).
Following this, allometries were also calculated for the overall data separated by (1) population and (2) minor or major morphs as determined by the breakpoint models applied to each population (see Appendix 3: Figure A5).

| Variation in sexual size dimorphism (SSD)
To test if SSD varied between populations, we ran linear models and found that the best fitted model with normally distributed errors included significant sex by location interaction, showing that SSD differed between populations (Table 2) Females from MY populations were significantly larger than females from SG populations while males from Langkawi Island MY were significantly larger than males from SG populations but did not differ significantly with Central Peninsular MY (Table 3). Body size also did not differ significantly between the SG populations and Central Peninsular MY (Table 3).

| Variation in male reproductive traits as a function of body size
Using a log-transformed data and the breakpoint model, a hyperallometric relationship (allometric coefficient, β > 1, Figure 5a, Table 4) was found between horn length and body size (pronotum width) for all four populations of O. babirussa. The adjusted R 2 values for equation one of the breakpoint models were high for all four populations, signaling a strong positive correlation. In addition, 95% confidence intervals (CIs) for equation 1 of all populations excluded zero, ruling out the likelihood of a zero slope, indicating a significant relationship between horn length and body size. These results suggest that body size is a significant factor in explaining horn length variation, where larger males have disproportionately longer horns. Interestingly, there is an overlap in CI values for all populations, which suggests that there were no significant population-level differences in allometric relationships (Figure 5a, Table 4). Overall analysis of horn length allometry separated by morphs found that both morphs showed hyperallometry, but minor morphs showed greater investment (β = 8) than major morphs (β = 2.1; Figure 5b).
On the contrary for postcopulatory traits, using log-transformed data, investments in both testes weight and sperm length increase somewhat, but the 95% confidence intervals overlap with both zero and unity (Figure 5c,e, Table 4). Thus, this increase in investment is not significant, and neither trait significantly deviates from

F I G U R E 3
Log-log scatterplot to determine the allometric relationship between horn length and body size (pronotum width) in male O. babirussa from Singapore. Following recommendations by Knell (2009), we fitted (a) linear model, (b) quadratic model, (c) cubic model, and (d) breakpoint model using the R package segmented (Muggeo, 2008) to the pooled data with all four populations to characterize the horn length-body size (n = 292).

| Sexual size dimorphism (SSD) varied among populations
Overall, our results showed that there is significant male-biased SSD in all populations except Central Peninsular MY, and that there is significant investment in precopulatory weapons, but no clear trend observed regarding investment in postcopulatory sexual traits.
Since mammals are of lower abundance and diversity in Singapore's forests, dung resources in Singapore are scarcer and less diverse, possibly leading to greater competition between males and higher sexual selection pressure. We thus hypothesized that malebiased SSD would be greater in Singapore than Peninsular Malaysia, and our results mostly agree with the hypothesis. Before discussing SSD, however, we first must address the finding that average body size and specifically female body size were found to be much smaller in both SG populations as compared with both populations from MY.
This disparity could be due to weaker fecundity selection on females in SG. We had planned to study this by examining female fecundity via measuring the spermathecae or rearing females and measuring clutch sizes, but we were unable to obtain enough data for either.  (Figure 4).
Theory suggests that strong precopulatory sexual selection drives male-biased SSD in insects as larger body size in males has been widely documented to increase mating success due to female choice or male-male competition (Blanckenhorn, 2005;Puniamoorthy et al., 2012;Stillwell et al., 2010). In many Onthophagus dung beetles and related taxa, males compete to gain access to females and body size is a predominant factor in determining fighting success (Emlen, 1997a(Emlen, , 1997bMoczek & Emlen, 1999 (Stillwell et al., 2010). As such, the male-biased SSD in O. babirussa is likely a derived trait that can be due to a relative increase in the intensity of sexual selection on male body size in this species. Our results also showed a strong investment in horns, a precopulatory weapon, further supporting that strong sexual selection is acting on males in this species via male-male competition.
One possible factor that could contribute to both lower female fecundity and stronger sexual selection on males in SG is resource availability, specifically dung resource. In SG, approximately 95% of forests were cleared over the last 200 years due to urbanization, causing high local extinctions of fauna such as birds and mammals in forest habitats (Bickford et al., 2010;Brook et al., 2003).
Singapore's remaining forests are mostly degraded, highly frag- presence of more and larger species provides female dung beetles more food and brood resources for oviposition opportunities (Qie et al., 2011;. Hence, lesser food resources in SG suggest that there could be a stronger viability selection on Singapore populations. On its own, this should lead to both males and females being smaller since viability selection acts on both sexes. However, fewer resources could also lead to greater intraspecific competition, especially between males competing over access to resources in order to gain access to potential mates.
The intensity of sexual selection on males could be strong enough to counteract the viability selection selecting for smaller body size, thus resulting in extreme male-biased SSD and males from SG reaching similar sizes to those from MY. Smaller females may produce fewer offspring but will still pass on their genes nonetheless, while smaller males may not even get an opportunity to mate. Body size could thus be such an important trait for males in SG that even under resource limitation, a minimum male body size must be achieved to even stand a chance in finding and securing a mate.
Alternative hypotheses to resource limitation that could affect body size and SSD differences between populations include Our results and the above discussion cover potential ultimate forces such as viability and sexual selection and how they could mediate differences in body size. Equally crucial factors to examine are potential proximate mechanisms driving these differences (Beckers et al., 2015). Based on our current findings, it is impossible to tell if the larger male body sizes in SG populations are due to genetic or environmental effects, such as differential gene expression or differential maternal resource partitioning to offspring of different  (Shafiei et al., 2001). If mothers from SG populations allocate more dung in the construction of male offspring brood balls than that of females, sex-biased differential maternal investment in offspring could be the driving proximate mechanism of male-biased SSD. If no significant differences are found in maternal investment, it is likelier that there is a genetic component such as differential gene expression between the sexes at play.
This study has shown that based on differing degrees of SSD  Despite the importance of possessing larger horns in gaining access to females, males with small body sizes and small horns were still regularly sampled and seem to persist in wild populations (Figure 5a).

| Investment in precopulatory and postcopulatory traits varied among populations
Besides common underlying causes for smaller body and horn size such as food limitation and larval competition, small-horned males of many Onthophagus species utilize alternative mating strategies in which they masquerade as females to sneak past guarding males with larger horns to gain access to breeding females (Beckers et al., 2017; Moczek & Emlen, 2000;Simmons & Ridsdill-Smith, 2011). Such an alternative mating strategy may exist in O. babirussa, which could explain the phenotypic variation in horn length observed in wildcaught populations (Moczek & Emlen, 2000).
Due to limited resources for growth and development, there may potentially be trade-offs in the investment of precopulatory and postcopulatory traits (Moczek & Nijhout, 2004). As there was a high relative investment in horn length, a precopulatory trait, we hypothesized that there would be a low relative investment in postcopulatory traits such as testes weight and sperm length. We would also then expect a lower allometric coefficient compared with horn length allometry. However, our results do not show a clear relationship between body size and both testes weight and sperm length across all populations. Looking at the data separated by minor and major morphs (Figure 5d,f), however, some trends can be observed. Testes weight for minor morphs showed a negative allometry, while major morphs showed slight hyperallometry. This could show morph-specific investment in postcopulatory traits, with minor morphs prioritizing investment in precopulatory traits, while major morphs can afford to invest in postcopulatory traits. This is supported by the much greater horn length allometric coefficient observed in minor males relative to major males ( Figure 5b).
Sperm length for both morphs was hypoallometric, but major males also showed a slightly steeper allometry and thus more relative investment in this postcopulatory trait. Overall, our findings suggest that investment in horns is more important, suggesting a lower relative investment in sperm length and testes weight than horns, which could be due to weaker postcopulatory selection in male O.
babirussa. Horns could be so important for mate acquisition that smaller, minor males prioritize investment in horns at the expense of postcopulatory investment, while major males could be at a com- It is also interesting to note that House and Simmons (2007) showed that in Onthophagus taurus, horn length allometry varied significantly with dung resource quality, while male genitalia exhib-

| CON CLUS IONS
This study reports population-level differences in SSD in the species We also thank members of the Reproductive Evolution Lab for assistance and support during this project.

CO N FLI C T O F I NTE R E S T
All authors certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers' bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The R scripts and data that support the findings of this study are openly available in Dryad at Yap, Sean; Toh, Kai Xin (2022)  thus able to differentiate between species that are closely related (Fraija-Fernández et al., 2018;Hebert et al., 2003;Waugh, 2007).

M O R PH O LO G I C A L SO RTI N G
Furthermore, the COI gene is an ideal species marker for insects due to the simple sequence alignment, primer sites that are robust and easily available and low likelihood of recombination and possession of introns (Foottit & Adler, 2009). A 3% threshold for uncorrected pairwise distances was used as it is commonly used in literature to F I G U R E A 1 (a) Diagram and (b) photograph showing structure and components of dung-baited pitfall traps. Human dung was wrapped in cloth to form a "dung ball" of 4 cm in diameter and suspended approximately 4 cm above a buried plastic cup using cotton twine and a shelter made of 15 cm by 15 cm corrugated board. Buried plastic cups were filled with 4 cm of water to wet the wings of dung beetles that have fallen in to prevent the beetles from escaping. differentiate species (Hebert et al., 2003;Meiklejohn et al., 2011;Srivathsan & Meier, 2012).
For all specimens, the right mid femur was dissected into 7 μl of QuickExtract solution and the DNA was extracted by following Lucigen's (the manufacturer's) protocol (Lucigen, 2018). After DNA extraction, COI amplification was conducted on extracted DNA. The and lastly, final extension (72°C, 5 min). In addition, primers were labelled with a sequence tag of 7-9 bp such that all specimens will have a unique combination of labelled forward and reverse primers.
The successfully amplified PCR products were pooled and sent for NGS sequencing, and sequences were analysed by constructing a cluster fusion diagram using uncorrected pairwise distances.

N G S A N D S EQ U E N CE A N A LYS I S
According to the manufacturer's instructions, the PCR products underwent bead cleanup to purify the products. Subsequently, the products were quantified, pooled in equimolar ratios and submitted for library preparation by Genome Institute of Singapore. Then, Next-Generation Sequencing (NGS) was conducted using MiSeq sequencing platform to obtain 313 bp fragment of COI gene.
Sequence analysis was then conducted with reference to the analysis pipeline detailed by (Meier et al., 2016). The reads from the pairedends were merged by the software PEAR 0.9.6 (Zhang et al., 2014).
The reads of each PCR products were then matched to their specific template specimen which was achieved due to primer pair combination that were uniquely labelled. A python script by (Srivathsan, unpublished) was used to (1) demultiplex data, (2) tally the reads for each sample, (3) identify and cluster identical reads into groups, (4) identify dominant groups of reads and combine with variants that were otherwise of identical length and lastly (5)

F I G U R E A 4
Cluster fusion diagram constructed based on uncorrected pairwise distances between COI barcode sequences from 26 representative specimens from across the main sampling sites.
next highest identity (Meier et al., 2016). Quality control was carried out by a set of criteria namely more than 50× read count, more than 10× barcode count and for the number of dominant reads to be five times or more than second most dominant reads (Meier et al., 2016).
This was to ensure that coverage attributed to each barcode was sufficient and not from confounding sequences such as contaminant DNA fragments. In addition, quality control rejects dominant sequences that may have arisen out of amplification error in the PCR step. Next, the sequences that passed the quality control were entered into the search query in Basic Local Alignment Search Tool (BLAST) to search for sequences that match >97% to non-Onthophagus taxa, which were contaminant sequences and thus eliminated from analysis. After quality control, MEGA7, an online software, was used to align the sequences to ensure that there were no stop codons. Then, a new Python script (Srivathsan, unpublished) was used to construct a cluster fusion diagram based on uncorrected pairwise distances, and a threshold of 3% was used to delimit species. This 3% pairwise distance threshold is widely used to distinguish between insect species in literature (Hebert et al., 2003;Srivathsan & Meier, 2012). A cluster fusion diagram of a subset of O. babirussa haplotypes from across the sampling sites is shown in Figure A4 below. The cluster fusion shows that there is little geographic pattern in the distribution of haplotypes, and that all haplotypes fall under the 3% pairwise distance threshold, supporting our morphological sorting of the specimens under the same species.
F I G U R E A 5 log-log scatterplot to determine the allometric relationship between horn length and body size (pronotum width) in male O. babirussa from (a) Central Catchment SG (n = 45), (b) Pulau Ubin Island SG (n = 61), (c) Central Peninsular MY (n = 46) and (d) Langkawi Island MY (n = 139), using the breakpoint model.