TY - JOUR AU - Dean, Matthew D AB - Abstract The rapid divergence of male genitalia is a preeminent evolutionary pattern. This rapid divergence is especially striking in the baculum, a bone that occurs in the penis of many mammalian species. Closely related species often display diverse baculum morphology where no other morphological differences can be discerned. While this fundamental pattern of evolution has been appreciated at the level of gross morphology, nearly nothing is known about the genetic basis of size and shape divergence. Quantifying the genetic basis of baculum size and shape variation has been difficult because these structures generally lack obvious landmarks, so comparing them in three dimensions is not straightforward. Here, we develop a novel morphometric approach to quantify size and shape variation from three-dimensional micro-CT scans taken from 369 bacula, representing 75 distinct strains of the BXD family of mice. We identify two quantitative trait loci (QTL) that explain ∼50% of the variance in baculum size, and a third QTL that explains more than 20% of the variance in shape. Together, our study demonstrates that baculum morphology may diverge relatively easily, with mutations at a few loci of large effect that independently modulate size and shape. Based on a combination of bioinformatic investigations and new data on RNA expression, we prioritized these QTL to 16 candidate genes, which have hypothesized roles in bone morphogenesis and may enable future genetic manipulation of baculum morphology. baculum, sexual selection, shape, size The rapid divergence of male genitalia is an almost universal evolutionary pattern among sexually reproducing organisms (Hosken and Stockley 2004; Eberhard 2009, 1985). Rapid evolution extends to the baculum, a bone in the penis of many mammalian species. Even a cursory examination of over 200 illustrations by Burt (1960) reveals a bewildering array of baculum diversity in size, shape, amorphousness, symmetry, and even the number of bones and patterns of branching. Mammalogists have long recognized the utility of the baculum in identifying species that are otherwise morphologically indistinguishable (i.e., Herdina et al. 2014; Thomas 1915), even positing that baculum morphology reinforces reproductive isolation between species (Patterson and Thaeler 1982, but see Good et al. 2003). The selective forces that drive the evolutionary novelty of the baculum are likely to include a combination of female choice, male-male competition, and conflict between male and female reproductive interests (House and Simmons 2003, 2005; Eberhard 2009, 1996, 2000; Breseño and Eberhard 2009; Arnqvist 1998; Arnqvist and Rowe 2005; Perry and Rowe 2012; Rönn et al. 2007; Rowe and Arnqvist 2012; Schulte-Hostedde et al. 2011; Hosken and Stockley 2004; Tasikas et al. 2009; Patterson 1983; Patterson and Thaeler 1982; Long and Frank 1968; Ramm 2007; Brennan et al. 2007, 2010; Arnqvist and Danielsson 1999; Higginson et al. 2012). The enormous interspecific divergence in baculum morphology, coupled with relatively low levels of intraspecific variation (e.g., Ramm et al. 2010; Klaczko et al. 2015), suggests that male genitalia are the targets of recurrent adaptive evolution. Several studies have linked baculum characteristics with male reproductive success. In seminatural enclosures of multiple male and female house mice, males with a wider baculum sired both more and larger litters (Stockley et al. 2013). Similarly, in wild populations and experimental populations subjected to higher levels of sperm competition, males developed bacula that were significantly wider than males from low sperm competition populations (Simmons and Firman 2014). Dominant males tend to develop wider bacula (Lemaître et al. 2012), suggesting additional links between baculum morphology and presumably adaptive behavioral phenotypes. Since sperm competition is a regular feature of mouse mating ecology (Dean et al. 2006; Firman and Simmons 2008), baculum size and shape are probably important components of fitness in natural populations. In mice, the baculum resides in the distal portion of the glans penis which enters the vagina during copulation (Rodriguez et al. 2011). Intromission may be an important aspect of “copulatory courtship,” which influences fertilization success beyond delivery of ejaculate (Adler and Toner 1986; Matthews and Adler 1978; Toner and Adler 1986; Toner et al. 1987; Diamond 1970), and it is possible that the baculum plays a role in these processes. The puzzling diversity of baculum morphology and its link to male reproductive fitness demand more attention. A genetic dissection of the ultimate sources of variation in baculum morphology has been hindered by two main obstacles. First, bacula lack true landmarks such as sutures, foramina, or processes. As a result, most studies summarize baculum complexity with simple length, width, and/or weight measurements (Long and Frank 1968). Some studies have captured baculum outlines using more modern computational techniques (Simmons and Firman 2014), but ideally, baculum size and shape would be captured in three dimensions. Here, we overcome this challenge using tools in computational geometry and geometric morphometrics to measure bacula without reliance on true landmarks, which should be applicable to many different biological structures. Second, although variation in baculum morphology is heritable (Simmons and Firman 2014), no study to date has mapped the genetic basis of variation in size or shape. Both genetic causes and molecular mechanisms underlying its rapid evolution are uncharted. The rapid evolutionary divergence discussed above predicts a relatively simple genetic architecture in which mutations in a few key loci lead to large changes in morphology. Furthermore, we might predict that size and shape would be modulated by different loci, so that divergence in one aspect is uncoupled from divergence in the other. Here, we unite our newly developed morphometric methods with the power of mouse genetics in a quantitative genetics framework, using two genetic models: a family of recombinant inbred lines (RILs) known as the BXDs (Taylor et al. 1973, 1999; Peirce et al. 2004), and a family of advanced intercross lines (AILs) known as the LGxSM (reviewed in Nikolskiy et al. 2015). Our study makes several advances toward understanding the genetic basis of baculum size and shape. We demonstrate that variation in baculum size and shape is heritable, and is controlled by a small number of QTL with comparatively large effects. We identified different QTL affecting size vs. shape that, when combined with new data on gene expression and bioinformatic filtering, highlight several compelling candidate genes. Future molecular analyses should eventually lead to a better understanding of genetic mechanisms of recurrent adaptive morphological evolution. Materials and Methods Specimens All protocols and personnel were approved by USC’s Institute for Animal Care and Use Committee (protocol #11394). Our study was based on two powerful mouse genetic models, the BXD RILs and the LGxSM AILs. Most mice were already euthanized and frozen as part of unrelated research programs; a few were raised in-house or ordered through common vendors, then euthanized via carbon dioxide exposure followed by cervical dislocation. The BXD RILs began in the 1970s with approximately 30 RILs (Taylor et al. 1973, 1999) and were expanded to ∼150 lines (Taylor et al. 1999; Peirce et al. 2004; R. W. Williams, unpublished results). All strains have been genotyped at over 3000 informative markers (www.genenetwork.org). Each BXD RIL began as a cross between two classical inbred strains, C57BL/6J (B6) and DBA/2J (D2). BXD1 through BXD42 were then maintained via brother-sister mating. BXD43 and higher were maintained through nonsib matings until generation 8–12 to accumulate larger numbers of recombinations and then inbred, making them an advanced intercross-rederived RIL family (Peirce et al. 2004), to create homozygous strains that on average have 50% B6 and 50% D2 genome, but with a different 50% segregating across different BXD RILs dependent upon the randomness of recombination during inbreeding. Individuals from the same BXD RIL are essentially genetically identical, and can be considered biological replicates of the same genotype. We sampled males that were 60–200 d old since the baculum is fully developed at this age (Glucksmann et al. 1976). The genomes of B6 and D2 have been nearly completely sequenced (Keane et al. 2011; Waterston et al. 2002; Wang et al. 2016), providing power to followup investigations of QTL. The LGxSM AILs began from crosses between the LG/J and SM/J strains, originally selected over many generations for large and small body size, respectively (Goodale 1938; MacArthur 1944). LGxSM AILs were not intentionally inbred, but rather maintained through random breeding of unrelated individuals. Approximately 100 families are created each generation by pairing males and females that are not full siblings. Each generation accumulates more recombination events that break up the genomes of the two parental strains while avoiding inbreeding. Published studies of the LGxSM AILs exist at the F2, F8, and F34 generations (Norgard et al. 2011; Parker et al. 2011, 2014, 2012; Cheng et al. 2010). As part of unrelated research, several hundred individuals of the F43 and F44 generations were genotyped at thousands of markers (J. M. Cheverud, unpublished results). These F43 and F44 mice were sampled here. Unfortunately, many of the frozen carcasses were missing the penis as a result of prior necropsy, and only 144 males could be sampled. The parental LG/J and SM/J genomes have been sequenced (Nikolskiy et al. 2015). For all specimens, the glans penis was cut proximal to the baculum, then placed in distilled water at room temperature for approximately 1 wk, at which point tissue was easily teased away using forceps and pressure of liquid flow from a squeeze bottle of 70% ethanol. Bacula were soaked in 200 μL of 1% ammonia solution for a maximum of 2 hr to clean any remaining tissue and remove oils, then dried and packed into thin layers of Aquafoam (Kokomo, IN) for micro-CT scanning. There is a distinct mass of fibrocartilage distal to the bony baculum; our study did not include this structure. For the remainder of the manuscript, “baculum” refers to the single bone in the baculum, excluding the distal fibrocartilage (Figure 1). Figure 1 Open in new tabDownload slide Visual representation of our morphometric pipeline for defining semilandmarks showing the two “parental extremes.” See text for details. Dorsal (A and B) and lateral (C and D) views of a typical D2 (A and C) or B6 (B, D) individual (specimens #DBA_1 and #C57_1, respectively). The gray background of each bone represents approximately 175,000 x-y-z points segmented from micro-CT scans. Moving from proximal (left sides of A–D) to distal (right sides of A–D), we sampled 50 slices. One slice is shown in more detail (E), which is looking down the center of the bone, displaying an empty internal medullary cavity. Exactly seven points were defined on the ventral and dorsal surfaces, as well as the leftmost and rightmost, totalling 16 semilandmarks per slice, indicated by spheres. The colors of the spheres indicate their contribution to shape differences (LD1) between D2 and B6 with red indicating regions that differ, and blue indicate more similar regions in comparisons between the two parental strains. Horizontal black axes indicate the z-axis (A–D) or the x-axis (E). Vertical black axes indicate the x-axis (A and B), or the y-axis (C–E). Micro-CT scanning of bacula Aquafoam slices containing bacula were stacked in a 35 mm cylindrical micro-CT sample holder. Micro-CT scans were acquired using a uCT50 scanner (Scanco Medical AG, Bruttisellen, Switzerland) at the USC imaging core facility, under the following settings: 90 kVp, 155 uA, 0.5 mm Al filter, 750 projections per 180° (360° coverage), exposure time of 500 msec, and voxel size of 15.5 μm. 3D transformation to “align” bacula Customized scripts were written in Python (www.python.org) and R (R Core Team 2014) to process microCT images, segment individual bacula, perform transformations, and call semilandmarks; all scripts and microCT images have been deposited in the Figshare Data repository (DOI:10.6084/m9.figshare.3080725.v1). Each baculum was segmented out of the z-stack of micro-CT images and converted to a 3D point cloud of x-y-z points. Each pixel in the point cloud represents bone, and the images include internal structures, not just the surface. Point clouds were transformed in three main steps borrowing methodology from Dines et al. (2014) and are graphically illustrated in Figure 1. First, the two points furthest apart in the point cloud were used to initially define a distal-proximal axis, which we consider a z-axis. Second, the 10% most proximal and the 10% most distal points were sliced out separately, and the centroids of their convex hulls calculated. The point cloud was then transformed such that the proximal centroid became 0, 0, 0 and the distal point cloud became 0, 0, z (where z was some positive value). This second transformation effectively controlled for slight variation associated with defining points in the first step. Third, we sliced out points that fell 15.00–15.25% along the length of the new z-axis and computed its minimum bounding rectangle (MBR) (Butterworth 2013). The entire point cloud was then retransformed so that the long end of this MBR ran parallel to a new x-axis, and the short end parallel to a new y-axis. Slight curvatures in the bone, especially toward the distal end, were exploited to ensure the dorsal-ventral axis was correct. All bone transformations were visually confirmed using the RGL library in R with customized R scripts. Defining semilandmarks The baculum lacks any true landmarks, so we mathematically defined 802 points along each bone (DOI:10.6084/m9.figshare.3080725.v1). These points can be thought of as “semilandmarks” (Mitteroecker and Gunz 2009; Bookstein 1997, 1991) since they correspond to regionally homologous regions across all specimens. Fifty point slices, each with thickness 0.25% the length of the z-axis, were evenly spaced along the anterior-posterior axis. One such slice is shown in detail in Figure 1E. Each slice was divided along the x-axis by seven lines running parallel to the y-axis. Points within 4% of the width of the slice to each line were projected onto that line, then the ventralmost and dorsalmost points defined as semilandmarks, specifically labeled according to slice and position. The leftmost and rightmost points of each slice were also defined as semilandmarks, for a total of 16 semilandmarks per slice (Figure 1E), and a total of 800 semilandmarks across the 50 slices. The single anterior-most and single posterior-most points of the baculum were added, bringing the total to 802 (Figure 1, A–D). Quantifying size variation Size was quantified as centroid size, the square root of the sum of squared distances of these 802 semilandmarks from their centroid. Counting the number of pixels per scan, or even the density of pixels (number of pixels divided by centroid size) yielded identical results, but we present those based on centroid size for simplicity. Quantifying shape variation We quantified shape difference between all possible pairs of bacula in a Generalized Procrustes framework, which standardized each set of 802 semilandmarks to a common size, translated them to a common origin, then optimally rotated their coordinates to minimize their Procrustes distance (Rohlf and Bookstein 1990; Slice 2007). During Procrustes superimposition, semilandmarks were allowed to “slide” along the bones’ surfaces using the function gpagen in the R package geomorph (Adams and Otárola-Castillo 2012), which improves alignment of corresponding anatomical regions lacking individual landmark homology (Bookstein 1997, 1996; Mitteroecker and Gunz 2009; Gunz et al. 2005). In short, sliding semilandmarks accompany uncertainty in specific placement of landmarks. The gpagen analysis resulted in a pairwise distance matrix (in the case of the BXD RILs, a 369 × 369 matrix; in the case of the LGxSM AILs, a 144 × 144 matrix). Our goal was to define a single metric that summarized each specimen’s shape in the context of the parental strains. Therefore, using only specimens of the parental strains (in the case of BXD RILs, 27 B6, and 18 D2 individuals; in the case of LGxSM AILs, 10 LG, and 9 SM individuals), we performed a linear discriminant analysis on the pairwise distance matrix, using the lda function in the R package mass (Ripley 2002) and parental strain as the separating factor. With two parental strains, a single linear discriminant function exists. We then projected the entire pairwise distance matrix into the space defined by this LD1, using the predict function. Repeatability To assess the repeatability of our micro-CT scanning and computational geometric manipulations, we scanned 34 bacula 81 times (32 bacula scanned twice, one baculum scanned six times, and one baculum scanned 11 times), each time removing the specimen and reloading it into Aquafoam and the micro-CT holder. For size, we calculated repeatability as the 1 – the median coefficient of variation of centroid size (unbiased standard deviation divided by the mean) across specimens. For shape, we could not calculate repeatability of LD1 because the measure spans 0, meaning we could be dividing by 0 to calculate coefficients of variation. Instead, we took the pairwise shape distance of each replicate bone to the other 368 bones in the dataset, then calculated the coefficients of variation per pairwise comparison, averaging across pairwise comparisons. Environmental input To assess the effect of environment, we focused on 28 B6 individuals (10 from the Jackson Laboratory, nine from the Levitt Lab at USC, nine from the Williams/Lu Lab) and 17 D2 (seven from the Jackson Laboratory, 10 from the Levitt Lab at USC). We performed an ANOVA with lab origin nested within strain for these 45 parental individuals. There were no other strains for which we sampled substantial number of individuals from different labs. Heritability Heritability, the proportion of phenotypic variance explained by genetic variance of size and shape measurements was estimated using a one-way ANOVA to summarize the amount of variance explained by strain. Among our 73 BXD RILs, we collected at least three males for 60. Heritabilities were estimated by focusing on these 60 (though the QTL mapping described below included all 73 RILs). The proportion of variance explained by strain identity was taken as the heritability. Mapping QTL We employed two main analyses to map loci affecting baculum size (centroid size) and shape (LD1). For the BXD RILs, phenotypic values were first averaged for any RILs for which we had sampled multiple individuals. Genotypes were downloaded from GeneNetwork, and included 3805 markers (2.2% mean missing data per RIL) spaced at a mean 0.54 cM throughout the genome. We kept the 2953 markers that were also part of a relatively recent effort to update the mouse genetic map (Cox et al. 2009). After removing parental strains (which are uninformative since they lack recombinant chromosomes), we employed the scanone function in the R package qtl, using Hayley-Knott regression to estimate the location and the effect of QTL (Broman and Sen 2009). To determine significance, we permuted phenotypes and genotypes 1000 times, and took the 95th quantile of the 1000 maximum LOD scores as our empirical significance threshold. We estimated the confidence intervals of any significant QTL using the lodint function in the R package qtl, by dropping 1.5 LOD units from the maximal LOD unit on the chromosome. Age and weight could not be included as covariates in analyses of the BXD RILs because not all labs that donated carcasses collected this information, and where they did occur they often varied dramatically given the individuals derived from very diverse research programs. Ignoring age and weight should only add noise to our analysis and inflate Type II errors. Unlike the BXD RILs, the LGxSM AILs were not maintained via brother-sister mating, so no two individuals were genetically identical but were also not equally related (Darvasi 1998; Darvasi and Soller 1995). The LGxSM AILs were analyzed using the R package qtlrel (Cheng et al. 2011), which accounts for background genetic relatedness prior to scanning for QTL. Genetic relatedness was estimated using all markers except the markers on the same chromosome as the one being analyzed (the “marker” option). Unlike the BXD RILs, age and weight could be included as covariates in our analyses of LGxSM AILs, because all individuals derived from the Cheverud Lab and were reared under identical conditions. Individuals were genotyped at 4716 diagnostic markers (0.1% mean missing data per specimen) spaced at a mean 0.36 F2 cM (J. M. Cheverud, unpublished data). Significance thresholds were estimated with 1000 permutations. Recombination distances were calculated from Cox et al. (2009). RNA sequencing To potentially narrow down QTL identified in the BXD RILs to candidate genes, we generated RNA-seq data from 10 5-wk-old B6 (N = 5) and D2 (N = 5) individuals. For logistical reasons, our B6 bacula were derived from C57BL/6N, which is a substrain that is over 99.9% genetically identical to C57BL/6J. The baculum is still not fully ossified at this point (Glucksmann et al. 1976), so genes expressed may ultimately affect its size and shape. Focusing on 5-wk-old males admittedly overlooks potentially important expression patterns that occur prior to or after this age. Therefore, we view these data as a means to “generate” hypotheses rather than definitively link genes under QTL to causality. Because we did not identify any significant QTL among the LGxSM AILs (see below), we only gathered RNA-seq data from the BXD parental strains. Bacula were homogenized in liquid nitrogen, then placed immediately into Trizol according to the Direct-zol RNA MiniPrep (Zymo Research R2052) protocol. RNA integrity was verified by Experion analysis (Bio-Rad). PolyA selection was carried out using Illumina Truseq V2 polyA beads. Libraries were prepared using Kapa Biosystems Stranded mRNA-Seq kit. We performed 12 PCR cycles to amplify the libraries, which were then visualized by Bioanalyzer analysis (Agilent) and quantified by qPCR (Kapa Biosystems Illumina Library Quantification Kit). Sequencing was performed on a NextSeq 500 using V2 chemistry. Paired-end 100 bp sequencing reads were mapped using Tophat (Trapnell et al. 2009) and aligned to the complete Mus musculus GRCm38.73 genome release available from Ensembl. Tophat v2.06 was run on a Linux x86 64-bit cluster with the following –read-edit-dist 8 –read-mismatches 8 –segment-mismatches 3 –min-anchor-length 12 –report-secondary-alignments. In addition, the GRCm38.73 General Transcript File (GTF file) was included with the “-G” and “–no-novel juncs” tags, ensuring that only known annotated exons were used. These mappings of full length and junction reads were subsequently used by Cufflinks (Trapnell et al. 2009) to generate gene counts. Cufflinks v2.1.0 was run on a Linux x86 64-bit cluster with the following parameters “–multi-read-correct–upper-quartile-norm–compatible-hits-norm–frag-bias-correct” and with the gene annotation included in the GTF file. Illumina sequencing generated 276 M paired end reads, of which 247 M (90%) mapped reads were used to access expression for the entire annotated transcriptome. A total of 13,926 genes had an FPKM value of at least 1 in at least one of the 10 specimens analyzed. Of these, 479 (610) were significantly differentially expressed at a corrected P = 0.05 (0.10) between the D2 and B6 parents. Because our main interest with the RNA-seq data was to generate hypotheses, we were willing to accept increased Type I errors (false rejection of a null hypothesis) in favor of reduced Type II errors (false acceptance of null), so we used a cutoff of P = 0.10 instead of the more traditional P = 0.05. From the RNA-seq data, we identified differentially expressed genes (P = 0.10 between parental strains after correction by Benjamini-Hochberg) as well as genes that were highly expressed in both parents (in the top 10% of genes expressed with FPKM at least 1). Even though this latter category may not show evidence of differential expression between parental strains, they could still lead to differential baculum development between parental strains, for example through posttranslational modification. Bioinformatics To further identify candidate genes under QTL, we characterized genetic variants that differed between the parental strains D2 and B6. Single nucleotide polymorphisms were downloaded using the Biomart tool from Ensembl version 81 (www.ensembl.org) and we identified genes with at least one nonsynonymous difference between B6 and D2 that also fell under our QTL. We also downloaded their estimates of dN/dS from “one-to-one” orthologs between mouse and rat to scan for genes with unusual rates of nonsynonymous evolution. Data availability Data and associated code from the entire pipeline described below are available on the Figshare Digital Repository (DOI:10.6084/m9.figshare.3080725.v1). Results We microCT-scanned 369 bacula representing 73 different BXD RILs and the two parental strains B6 and D2. Repeatability was 0.99 for centroid size, and 0.97 for shape divergence, respectively. Lab origin did not significantly influence shape (F3,1 = 1.3, P = 0.28), but did contribute to size (F3,1 = 12.68, P < 0.0001), explaining 18% of the variance in the latter case. In summary, our methods show good repeatability and any contribution arising from lab origin should not introduce systematic biases into our analyses. BXD RILs, size Strain was a significant predictor of size variation (F59,244 = 14.8, P < 0.01), accounting for 78% of the variance (heritability = 0.78), with BXD RILs largely falling between the parental B6 and D2 strains indicative of mostly additive genetic variance (Figure 2). Phenotype means that fell outside either parent suggested transgressive segregation or epistasis. Figure 2 Open in new tabDownload slide Baculum size and shape variation among the 73 BXD recombinant inbred lines (RILs) and the two parental strains D2 and B6. Each point represents the mean size and shape for each strain, bars indicate standard errors where possible. B and D indicate the two parental strains of the BXD RILs. B = B6, D = D2. From the 73 BXD RILs, we detected two QTL for size (Figure 3A). One occurred on chromosome 1 between 136.3–150.1 Mb, with a maximum LOD score of 7.45, which was higher than all 1000 permutations (95% quantile = 3.63). The other occurred on chromosome 12, 58.0–76.4 Mb, with a maximum LOD score of 4.21, which was greater than 988 permutations (P = 0.012). The two markers closest to the maximal QTL peaks on chromosomes 1 or 12 showed a clear difference in centroid size among BXD RILs carrying a B6 allele vs. D2 allele (Figure 4, A and B). These two size QTL were not in linkage disequilbrium (χ2 = 0.16, df = 1, P = 0.7), as expected since they reside on different chromosomes. Figure 3 Open in new tabDownload slide Results of scanone analyses for (A) baculum size (centroid size) and (B) shape (LD1) in the BXD recombinant inbred lines RILs. Dashed line indicates significance threshold determined from 1000 permutations of genotype and phenotype, lines indicate LOD (logarithm of the odds) scores testing the null hypothesis of no QTL (quantitative trait loci) along 2,953 markers in the genome. Figure 4 Open in new tabDownload slide Baculum size (centroid size) and shape (LD1) from the 73 BXD recombinant inbred lines (RILs). The 73 BXD RILs are separated according to which parental allele they carry at the three QTL (quantitative trait loci) identified in Figure 3. The black diamonds represent the phenotypic value at each marker for the parental strains, C57 or DBA respectively. (A) rs6202860, the marker closest to the size QTL on chromosome 1. (B) rs3716547, the marker closest to the size QTL on chromosome 12. (C) rs13476871, the marker closest to the shape QTL on chromosome 2. BXD RILs, shape Strain was a significant predictor of LD1 variation (F59,244 = 18.1, P < 0.01), accounting for 81% of the variance (heritability = 0.81), with BXD RILs largely falling between the parental B6 and D2 strains indicative of mostly additive genetic variance (Figure 2). From the 73 BXD RILs, we detected a single QTL on chromosome 2, 154.7–161.7 Mb, with a maximum LOD score of 4.32, which was higher than 988/1000 permutations (P = 0.012) (Figure 3B). At the marker closest to the maximum QTL peak of chromosome 2, there was a clear difference in LD1 between the 33 BXD RILs with a B6 allele compared to the 40 with a D2 allele (Figure 4C). This shape QTL was not in linkage disequilibrium with either of the two size QTL (χ2 = 0.01, 0.32, df = 1, P = 0.9, 0.6, respectively), as expected since they reside on different chromosomes. BXD RILs, size and shape correlated Interestingly, the mean size and shape per strain were correlated (Pearson’s correlation coefficient = –0.35, P = 0.002, Figure 2), even though major-effect QTL affecting size and shape were found on different chromosomes. This correlation suggests that small bacula tend to have the straight shape of D2, while large bacula tend to have the dorsoventral curvature of B6. BXD RILs, RNA-seq and bioinformatics There were 111 (162) protein-coding genes that occurred under the size QTL identified on chromosome 1 (chromosome 12), of which 93 (145) were expressed, 3 (5) were highly expressed, 1 (2) was differentially expressed, and 11 (14) had at least one nonsynonymous variant between B6 and D2 (Table S1 and Table S2). There were 160 protein-coding genes that occurred under the shape QTL on chromosome 2, of which 140 were expressed, 12 highly expressed, 2 differentially expressed, and 27 had at least one nonsynonymous variant (Table S1 and Table S2). A literature review of all protein-coding genes that fell under any one of our three QTL (Figure 3) and were also highly expressed, differentially expressed, or had at least one nonsynonymous variant, revealed 16 genes with potentially interesting effects on bone morphology (Table S3). Four of these – Ptprc, 2700049A03Rik, Aspm, and Kif14 – showed relatively high pairwise dN/dS estimates when compared to their one-to-one ortholog in rat (all above the 75% genome-wide quantile). LGxSM AILs, size and shape From the LGxSM AILs, we microCT-scanned 144 bacula, including 36 unrelated F43 and 89 F44 from 52 unique families, as well as 10 LG and 9 SM parental strains. Although both centroid size and LD1 variance showed evidence of largely additive genetic contribution (Figure S1), no significant QTL were identified for either (Figure S2). Furthermore, there was no correlation between the proportion of an individual’s genome that was from the LG parent and either baculum size (centroid size) or shape (LD1) (Pearson’s correlation coefficient, P = 0.98, 0.67, respectively). This is probably an indication that our sample size for LGxSM was underpowered, a hypothesis supported by the lack of significant QTL or correlation between proportion LG/J alleles and body size (Pearson’s correlation coefficient, P = 0.75), a trait that is highly heritable in larger studies of LGxSM intercrosses (Kramer et al. 1998; Cheverud et al. 1996). Discussion Baculum size and shape diverge rapidly across species, with known effects on male reproductive success. By combining a novel morphometric pipeline with the power of mouse quantitative genetics, we make three main discoveries. First, a few QTL of major effect influence baculum variation in the BXD RILs, with two QTL explaining 50.6% of the variance in size, and a third QTL explaining 23.4% of the variance in shape. Interestingly, no QTL were detected in a second genetic system, the LGxSM AILs. Second, these major-effect QTL affecting size were distinct from the major-effect QTL affecting shape, even though the two traits were correlated. Third, by combining our QTL with RNA-seq data as well as bioinformatic features and literature review, we identified 16 promising candidate genes that may enable a deeper understanding of the genetic basis of their rapid divergence. Testing predictions that stem from rapid divergence Our finding of a few QTL with large effects suggests that morphological divergence of the baculum may accumulate easily, acting via mutations in a few key targets. Several studies (all in insects, most in Drosophila) have found QTL of large effect (>10% of phenotypic variance) on aspects of the male genital apparatus. It should be emphasized that like the current study, none of these studies have estimated effect sizes from specific genes. Many of these studies were initiated by crossing two closely related species, then backcrossing the F1s to either parental species, including two species of carabid beetles (Sasabe et al. 2007), Drosophila simulans and D. mauritiana (True et al. 1997; Zeng et al. 2000; Liu et al. 1996; LeVasseur-Viens and Moehring 2014; Tanaka et al. 2015), D. simulans and D. sechellia (Macdonald and Goldstein 1999), and D. yakuba and D. santomea (Peluffo et al. 2015). Although large-effect QTL were detected in all these studies, the extent of linkage that will exist in these backcross designs, coupled with the limited number of markers used to interrogate the genome, may underestimate the true number of QTL and thus inflate the percent variance explained by any single QTL. Furthermore, it is possible that interspecific studies are biased toward finding large-effect QTL since the different species employed were already known to harbor highly divergent male genitalia. Alternatively, we might predict a bias toward finding small-effect QTL if genetic effects are masked by disruption of normal development in hybrid males. In contrast, studies within species tend to find QTL with smaller effects. Three QTL contributed 4.7–10.7% to variance in the shape of the posterior lobe of genital arch within an advanced intercross panel of D. melanogaster (McNeil et al. 2011). Similarly, a genome-wide association test among 155 inbred strains of D. melanogaster strains of the Drosophila Genetic Resource Panel (Mackay et al. 2012) found multiple SNPs throughout the genome that correlated with small to moderate differences in the size and shape of the posterior lobe (Takahara and Takahashi 2015). Another study in D. montana showed mostly small effect QTL explaining less than 10% of the variance (Schäfer et al. 2011). Taken together, QTL identified in these intraspecific studies are more numerous, with smaller effect sizes, than the interspecific studies discussed above. In this sense, our study is unusual in that it was an intraspecific study but identified a few QTL of large effect, although it should be noted that a small amount of interspecific introgression was detected in both B6 and D2 (Yang et al. 2011). In contrast, a recent morphometric study of mouse skulls found many QTL of very small effect (Maga et al. 2015). Perhaps bacula are unique in that a few major QTL enable their rapid divergence. In our study, QTL affecting size and shape were independent of each other, another genetic characteristic that might allow for more rapid morphological divergence, since mutations may be “less pleiotropic.” Size and shape QTL could not be separated in many of the studies mentioned above (Liu et al. 1996; Macdonald and Goldstein 1999), as the low number of markers employed precluded delineation of multiple QTL. Two studies within D. melanogaster found distinct QTL affecting size vs. shape of male posterior lobes (Schäfer et al. 2011; Takahara and Takahashi 2015). Interestingly, we found that the major QTL affecting size and shape of the baculum were independent of each other (Figure 3), even though size and shape were correlated (Figure 2). One explanation is that many QTL of small effect were not detected here but sufficient to drive an overall correlation between size and shape. Two additional predictions could not be properly evaluated here. First, we predicted that different genetic backgrounds of mice might yield distinct QTL, meaning there are multiple genetic pathways that can generate baculum variation. Some support for this hypothesis comes from a comparison of two quantitative genetic studies that mapped nonoverlapping QTL in two different populations of D. melanogaster (Takahara and Takahashi 2015; McNeil et al. 2011). However, one study of a D. mauritiana–D. simulans backcross identified identical QTL as found in a D. mauritiana–D. sechellia backcross (LeVasseur-Viens and Moehring 2014). In the current study, we found three large-effect QTL in the BXD RILs (Figure 3), none of which were found in the LGxSM AIL panel (Figure S2). However, the lack of QTL in the LGxSM AIL panel may simply be due to lack of power. Second, we predicted that baculum variation might arise via alterations in gene expression since both early genital development and bone morphogenesis proceed from highly conserved molecular pathways where protein-coding changes would be expected to have dire consequences (Carroll 2008). Our data do not yet speak to this question; across the 16 candidate genes identified under our three QTL (Table S3), seven were highly or differentially expressed among parental strains of the BXD RILs, and 10 had at least one nonsynonymous mutation. Further evaluation of whether baculum variation arises from expression or structural prediction requires finer scale mapping, more detailed expression studies, and/or stronger evidence that certain genes are good candidates for baculum variation. Rapid genital divergence vs. conserved genital development pathways One of the paradoxes of the rapid evolution of male genitalia, at least in vertebrates, is the highly conserved set of genes that appear to be expressed during early genital development, including Sonic Hedgehog (Shh), various fibroblast growth factor receptors (Fgfr’s), and various homeobox-containing genes of the D cluster (Hoxd) (Haraguchi et al. 2001, 2000; Klonisch et al. 2004; Perriton et al. 2002; Miyagawa et al. 2009; Cohn 2011; Seifert et al. 2009; Gredler et al. 2015; Sanger et al. 2015; Infante et al. 2015). In fact, the outgrowth of external genitalia shares many genetic pathways with the development of digits (Warot et al. 1997; Tschopp et al. 2014; Lonfat et al. 2014; Cobb and Duboule 2005; Haraguchi et al. 2000; Kondo et al. 1997; Perriton et al. 2002; Dollé et al. 1991; Zákány and Duboule 1999) or gut (Cohn 2011), although genital-specific enhancers modify their expression (Lonfat et al. 2014). Hoxd-13 mutant mice have smaller bacula than wild-type mice (Hérault et al. 1997; Zákány and Duboule 1999), and differently shaped bacula (Peichel et al. 1997), and Hoxd genes have enhancers that drive expression specific to external genitalia (Dollé et al. 1991; Lonfat et al. 2014). When Hoxd-13, Hoxd-11, and Hoxd-12 were all experimentally made nonfunctional, no genitals developed (Kondo et al. 1997). However, all these genetically manipulated individuals show extreme dysmorphia in body plan, and many are lethal before birth. Importantly, none of the QTL that we observed here overlap with Shh, Fgfr’s, or various Hoxd genes, suggesting the variation we observe is not due to genetic variation in these conserved pathways, although it is formally possible that something under our QTL affects expression of these genes in trans. Candidate genes Our study offers the first set of candidate genes to explain variation in baculum size and shape. Although systematically testing which (if any) of these protein-coding genes explain baculum variance remains outside the scope of the current study, we prioritized genes based on three criteria: 1) they fell under QTL identified above (Figure 3), 2) they were highly expressed or differentially expressed in 5-wk-old bacula, or had at least one nonsynonymous variant between parental strains, and 3) literature searching revealed a potential link to bone or genital morphogenesis. The baculum forms after birth, through both direct ossification and ossification of cartilage anlage as the animal approaches sexual maturity (Murakami and Mizuno 1984), so bone-related phenotypes are potentially interesting. By these criteria, we identified 16 candidate genes (Table S3). Four genes under the size QTL of chromosome 1 are potentially interesting candidates. Kif14 and Aspm have been shown to cause cranial deformities and reduced body size (Fujikura et al. 2013; Pulvers et al. 2010). A third, Ptprc, alters bone morphology when knocked out (Shivtiel et al. 2008). Finally, mice deficient in Mgat2 exhibit reduced body size and reduced bone density due to hyperactive osteoclasts (Wang et al. 2001). Under the size QTL of chromosome 12, Dact1 has been shown to negatively regulate Wnt signaling (Wen et al. 2010), which in turn affects the development of bone and genitals (Zheng et al. 2012; Berendsen and Olsen 2015; Hu et al. 2005; Haraguchi et al. 2000, 2001; Yamaguchi et al. 1999). Wnt signaling has also been implicated in the divergence of male genitalia in flies (Tanaka et al. 2015). Mice missing a functional Dact1 allele lacked external genitals and displayed numerous skeletal abnormalities (Wen et al. 2010). Another gene under this QTL, Hif1a, also affects bone density (Wang et al. 2007). Two genes under the shape QTL of chromosome 2 influence bone development. Rbl1 is a key regulator of cell development which, if inactivated, results in shortened limbs, defective endochondral ossification, altered chondrocyte growth (Cobrinik et al. 1996), and reduced body size (Scimè et al. 2005). A second gene, Tgm2, plays a role in chondrogenesis and ultimately bone formation (reviewed in Iismaa et al. 2009). Conclusions Our study provides the first candidate genes to explain variation in baculum size and shape, a structure with an astonishing rate of evolutionary divergence. Our study provides evidence that baculum variation is explained by a few QTL of large effect that independently affect size and shape with apparently minimal pleiotropic side-effects, and may help reconcile the paradox of rapid morphological divergence coupled with conserved developmental pathways. Future experiments should focus on testing candidate genes identified here as a means to genetically dissect its function, and to provide a deeper understanding of the rapid evolution of this amazing structure. Acknowledgments The University of Southern California’s (USC’s) Center for High-Performance Computing and Communications provided computational resources. David Threadgill (Texas A & M University), Allison Knoll (USC), and Pat Levitt (USC) provided specimens. Micro-CT scanning was performed at the USC Small Animal Imaging Core; we are especially thankful to Anitha Priya and Tautis Skorka. Joe Dunham and Brent Young assisted with molecular wrangling. All RNA sequencing was done at the USC Epigenome Center; we are especially thankful to Selene Tyndale and Charlie Nicolet. Karl Broman, Peter Carbonetto, Peter Ralph, Kevin Casey, and Erik Otárola-Castillo provided many useful discussions of methodology. We thank the students of the spring 2015 semester of BISC 313, whose term papers helped reveal some of the candidate genes discussed here. Jeff Good, Josh Akey, D. J. de Koning, and one anonymous reviewer offered comments on the manuscript. This work was supported by National Institutes of Health grant no. GM098536 (to M.D.D.). Footnotes Supplemental material is available online at www.g3journal.org/lookup/suppl/doi:10.1534/g3.116.027888/-/DC1 Communicating editor: D. J. de Koning Literature Cited Adams D C , Otárola-Castillo E, 2012   Package ‘geomorph’: Geometric morphometric analysis of 2d/3d landmark data. R package version 1.0. The Comprehensive R Network , CRAN . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Adler N T , Toner J P, 1986   The effects of copulatory behavior on sperm transport and fertility in rats. Ann. N. Y. Acad. Sci. 474 ( 1 ): 21 – 32 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Arnqvist G , 1998   Comparative evidence for the evolution of genitalia by sexual selection. Nature 393 ( 6687 ): 784 – 786 . Google Scholar Crossref Search ADS WorldCat Arnqvist G , Danielsson I, 1999   Copulatory behavior, genital morphology, and male fertilization success in water striders. Evolution 53 : 147 – 156 . Google Scholar Crossref Search ADS PubMed WorldCat Arnqvist G , Rowe L, 2005   Sexual conflict , Princeton University Press , Princeton, New Jersey . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Berendsen A D , Olsen B R, 2015   Bone development. Bone 80 : 14 – 18 . Google Scholar Crossref Search ADS PubMed WorldCat Bookstein F , 1991   Morphometric tools for landmark data: geometry and biology , Cambridge University Press , New York . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Bookstein F J , 1997   Landmark methods for forms without landmarks: morphometrics of group differences in outline shape. Med. Image Anal. 1 : 225 – 243 . Google Scholar Crossref Search ADS PubMed WorldCat Bookstein F L , 1996   Biometrics, biomathematics and the morphometric synthesis. Bull. Math. Biol. 58 : 313 – 365 . Google Scholar Crossref Search ADS PubMed WorldCat Brennan P L , Prum R O, McCracken K G, Sorenson M D, Wilson R E et al. , 2007   Coevolution of male and female genital morphology in waterfowl. PLoS One 2 ( 5 ): e418 . Google Scholar Crossref Search ADS PubMed WorldCat Brennan P L , Clark C J, Prum R O, 2010   Explosive eversion and functional morphology of the duck penis supports sexual conflict in waterfowl genitalia. Poc. R. Soc. Lond. B Biol. Sci. 277 : 1309 – 1314 . Google Scholar OpenURL Placeholder Text WorldCat Breseño R D , Eberhard W G, 2009   Experimental modifications imply a stimulatory function for male tsetse fly genitalia, supporting cryptic female choice theory. J. Evol. Biol. 22 ( 7 ): 1516 – 1525 . Google Scholar Crossref Search ADS PubMed WorldCat Broman K W , Sen S, 2009   A guide to QTL mapping with R/qtl , Springer , New York . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Burt W , 1960   Bacula of North American mammals. Misc Publ Mus Zool Univ Mich 113 : 1 – 75 . Google Scholar OpenURL Placeholder Text WorldCat Butterworth D , 2013  https://github.com/dbworth/minimum-area-bounding-rectangle. Carroll S B , 2008   Evo-devo and an expanding evolutionary synthesis: a genetic theory of morphological evolution. Cell 134 ( 1 ): 25 – 36 . Google Scholar Crossref Search ADS PubMed WorldCat Cheng R , Lim J E, Samocha K E, Sokoloff G, Abney M et al. , 2010   Genome-wide association studies and the problem of relatedness among advanced intercross lines and other highly recombinant populations. Genetics 185 ( 3 ): 1033 – 1044 . Google Scholar Crossref Search ADS PubMed WorldCat Cheng R , Abney M, Palmer A A, Skol A D, 2011   QTLRel: an R package for genome-wide association studies in which relatedness is a concern. BMC Genet. 12 ( 1 ): 66 . Google Scholar Crossref Search ADS PubMed WorldCat Chessel D , Dufour A B, Thioulouse J, 2004   The ade4 package-I-One-table methods. R News 4 ( 1 ): 5 – 10 . Google Scholar OpenURL Placeholder Text WorldCat Cheverud J M , Routman E J, Duarte F, van Swinderen B, Cothran K et al. , 1996   Quantitative trait loci for murine growth. Genetics 142 ( 4 ): 1305 – 1319 . Google Scholar Crossref Search ADS PubMed WorldCat Cobb J , Duboule D, 2005   Comparative analysis of genes downstream of the Hoxd cluster in developing digits and external genitalia. Development 132 ( 13 ): 3055 – 3067 . Google Scholar Crossref Search ADS PubMed WorldCat Cobrinik D , Lee M H, Hannon G, Mulligan G, Bronson R T et al. , 1996   Shared role of the pRB-related p130 and p107 proteins in limb development. Genes Dev. 10 ( 13 ): 1633 – 1644 . Google Scholar Crossref Search ADS PubMed WorldCat Cohn M J , 2011   Development of the external genitalia: Conserved and divergent mechanisms of appendage patterning. Dev. Dyn. 240 ( 5 ): 1108 – 1115 . Google Scholar Crossref Search ADS PubMed WorldCat Cox A , Ackert-Bicknell C L, Dumont B L, Ding Y, Bell J T et al. , 2009   A new standard genetic map for the laboratory mouse. Genetics 182 ( 4 ): 1335 – 1344 . Google Scholar Crossref Search ADS PubMed WorldCat Darvasi A , 1998   Experimental strategies for the genetic dissection of complex traits in animal models. Nat. Genet. 18 ( 1 ): 19 – 24 . Google Scholar Crossref Search ADS PubMed WorldCat Darvasi A , Soller M, 1995   Advanced intercross lines, an experimental population for fine genetic mapping. Genetics 141 ( 3 ): 1199 – 1207 . Google Scholar Crossref Search ADS PubMed WorldCat Dean M D , Ardlie K G, Nachman M W, 2006   The frequency of multiple paternity suggests that sperm competition is common in house mice (Mus domesticus). Mol. Ecol. 15 : 4141 – 4151 . Google Scholar Crossref Search ADS PubMed WorldCat Diamond M , 1970   Intromission pattern and species vaginal code in relation to induction of pseudopregnancy. Science 169 ( 949 ): 995 – 997 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Dines J P , Otárola-Castillo E, Ralph P, Alas J, Daley T et al. , 2014   Sexual selection targets cetacean pelvic bones. Evolution 68 ( 11 ): 3296 – 3306 . Google Scholar Crossref Search ADS PubMed WorldCat Dollé P , Izpisúa-Belmonte J C, Brown J M, Tickle C, Duboule D, 1991   HOX-4 genes and the morphogenesis of mammalian genitalia. Genes Dev. 5 ( 10 ): 1767 – 1776 . Google Scholar Crossref Search ADS PubMed WorldCat Eberhard W , 2009   Evolution of genitalia: theories, evidence, and new directions. Genetica 138 ( 1 ): 5 – 18 . Google Scholar Crossref Search ADS WorldCat Eberhard W G , 1985   Sexual selection and animal genitalia , Harvard University Press , Cambridge, Mass. Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Eberhard W G , 1996   Female control: sexual selection by cryptic female choice , Princeton University Press , Princeton, New Jersey . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Eberhard W G , 2000   Criteria for demonstrating postcopulatory female choice. Evolution Int J Org Evolution 54 ( 3 ): 1047 – 1050 . Google Scholar Crossref Search ADS WorldCat Firman R C , Simmons L W, 2008   The frequency of multiple paternity predicts variation in testes size among island populations of house mice. J. Evol. Biol. 21 : 1524 – 1533 . Google Scholar Crossref Search ADS PubMed WorldCat Fujikura K , Setsu T, Tanigaki K, Abe T, Kiyonari H et al. , 2013   Kif14 mutation causes severe brain malformation and hypomyelination. PLoS One 8 ( 1 ): e53490 . Google Scholar Crossref Search ADS PubMed WorldCat Glucksmann A , Ooka-Souda S, Miura-Yasugi E, Mizuno T, 1976   The effect of neonatal treatment of male mice with antiandrogens and of females with androgens on the development of the os penis and os clitoridis. J. Anat. 121 ( Pt 2 ): 363 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Good J M , Demboski J R, Nagorsen D W, Sullivan J, 2003   Phylogeography and introgressive hybridization: chipmunks (genus Tamias) in the northern rocky mountains. Evolution 57 ( 8 ): 1900 – 1916 . Google Scholar Crossref Search ADS PubMed WorldCat Goodale H , 1938   A study of the inheritance of body weight in the albino mouse by selection. J. Hered. 29 ( 3 ): 101 – 112 . Google Scholar Crossref Search ADS WorldCat Gredler M L , Seifert A W, Cohn M J, 2015   Tissue-specific roles of Fgfr2 in development of the external genitalia. Development 142 ( 12 ): 2203 – 2212 . Google Scholar Crossref Search ADS PubMed WorldCat Gunz P , Mitteroecker P, Bookstein F, 2005   Semilandmarks in three dimensions , pp. 73 – 98 in Modern morphometrics in physical anthropology , edited by Slice D E Kluwer Press , New York . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Haraguchi R , Mo R, Hui C-c, Motoyama J, Makino S et al. , 2001   Unique functions of Sonic hedgehog signaling during external genitalia development. Development 128 ( 21 ): 4241 – 4250 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Haraguchi R , Suzuki K, Murakami R, Sakai M, Kamikawa M et al. , 2000   Molecular analysis of external genitalia formation: the role of fibroblast growth factor (Fgf) genes during genital tubercle formation. Development 127 ( 11 ): 2471 – 2479 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Hérault Y , Fraudeau N, Zákány J, Duboule D, 1997   Ulnaless (Ul), a regulatory mutation inducing both loss-of-function and gain-of-function of posterior Hoxd genes. Development 124 ( 18 ): 3493 – 3500 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Herdina A N , Hulva P, Horácek I, Benda P, Mayer C et al. , 2014   MicroCT imaging reveals morphometric baculum differences for discriminating the cryptic species Pipistrellus and P. pygmaeus. Acta Chiropt. 16 ( 1 ): 157 – 168 . Google Scholar Crossref Search ADS WorldCat Higginson D M , Miller K B, Segraves K A, Pitnick S, 2012   Female reproductive tract form drives the evolution of complex sperm morphology. Proc. Natl. Acad. Sci. USA 109 ( 12 ): 4538 – 4543 . Google Scholar Crossref Search ADS WorldCat Hosken D J , Stockley P, 2004   Sexual selection and genital evolution. Trends Ecol. Evol. 19 ( 2 ): 87 – 93 . Google Scholar Crossref Search ADS PubMed WorldCat House C M , Simmons L W, 2003   Genital morphology and fertilization success in the dung beetle Onthophagus taurus: an example of sexually selected male genitalia. Proc. Biol. Sci. 270 ( 1514 ): 447 – 455 . Google Scholar PubMed OpenURL Placeholder Text WorldCat House C M , Simmons L W, 2005   The evolution of male genitalia: patterns of genetic variation and covariation in the genital sclerites of the dung beetle Onthophagus taurus. J. Evol. Biol. 18 ( 5 ): 1281 – 1292 . Google Scholar Crossref Search ADS PubMed WorldCat Hu H , Hilton M J, Tu X, Yu K, Ornitz D M et al. , 2005   Sequential roles of Hedgehog and Wnt signaling in osteoblast development. Development 132 ( 1 ): 49 – 60 . Google Scholar Crossref Search ADS PubMed WorldCat Iismaa S E , Mearns B M, Lorand L, Graham R M, 2009   Transglutaminases and disease: lessons from genetically engineered mouse models and inherited disorders. Physiol. Rev. 89 ( 3 ): 991 – 1023 . Google Scholar Crossref Search ADS PubMed WorldCat Infante C R , Mihala A G, Park S, Wang J S, Johnson K K et al. , 2015   Shared Enhancer Activity in the Limbs and Phallus and Functional Divergence of a Limb-Genital cis-Regulatory Element in Snakes. Dev. Cell 35 ( 1 ): 107 – 119 . Google Scholar Crossref Search ADS PubMed WorldCat Keane T M , Goodstadt L, Danecek P, White M A, Wong K et al. , 2011   Mouse genomic variation and its effect on phenotypes and gene regulation. Nature 477 ( 7364 ): 289 – 294 . Google Scholar Crossref Search ADS PubMed WorldCat Klaczko J , Ingram T, Losos J, 2015   Genitals evolve faster than other traits in Anolis lizards. J. Zool. 295 ( 1 ): 44 – 48 . Google Scholar Crossref Search ADS WorldCat Klonisch T , Fowler P A, Hombach-Klonisch S, 2004   Molecular and genetic regulation of testis descent and external genitalia development. Dev. Biol. 270 ( 1 ): 1 – 18 . Google Scholar Crossref Search ADS PubMed WorldCat Kondo T , Zakany J, Innis J W, Duboule D, 1997   Of fingers, toes and penises. Nature 390 ( 6655 ): 29 . Google Scholar Crossref Search ADS PubMed WorldCat Kramer M G , Vaughn T T, Pletscher L S, King-Ellison K, Adams E et al. , 1998   Genetic variation in body weight gain and composition in the intercross of Large (LG/J) and Small (SM/J) inbred strains of mice. Genet. Mol. Biol. 21 ( 2 ): 211 – 218 . Google Scholar Crossref Search ADS WorldCat Lemaître J-F , Ramm S A, Jennings N, Stockley P, 2012   Genital morphology linked to social status in the bank vole (Myodes glareolus). Behav. Ecol. Sociobiol. 66 ( 1 ): 97 – 105 . Google Scholar Crossref Search ADS WorldCat LeVasseur-Viens H , Moehring A J, 2014   Individual Genetic Contributions to Genital Shape Variation between Drosophila simulans and D. mauritiana. Int. J. Evol. Biol. 2014 : 808247 . Google Scholar Crossref Search ADS PubMed WorldCat Liu J , Mercer J M, Stam L F, Gibson G C, Zeng Z B et al. , 1996   Genetic analysis of a morphological shape difference in the male genitalia of Drosophila simulans and D. mauritiana. Genetics 142 ( 4 ): 1129 – 1145 . Google Scholar Crossref Search ADS PubMed WorldCat Lonfat N , Montavon T, Darbellay F, Gitto S, Duboule D, 2014   Convergent evolution of complex regulatory landscapes and pleiotropy at Hox loci. Science 346 ( 6212 ): 1004 – 1006 . Google Scholar Crossref Search ADS PubMed WorldCat Long C A , Frank T, 1968   Morphometric variation and function in the baculum, with comments on correlation of parts. J. Mammal. 49 ( 1 ): 32 – 43 . Google Scholar Crossref Search ADS WorldCat MacArthur J W , 1944   Genetics of body size and related characters. I. Selecting small and large races of the laboratory mouse. Am. Nat. 78 : 142 – 157 . Google Scholar Crossref Search ADS WorldCat Macdonald S J , Goldstein D B, 1999   A quantitative genetic analysis of male sexual traits distinguishing the sibling species Drosophila simulans and D. sechellia. Genetics 153 ( 4 ): 1683 – 1699 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Mackay T F C , Richards S, Stone E A, Barbadilla A, Ayroles J F et al. , 2012   The Drosophila melanogaster Genetic Reference Panel. Nature 482 ( 7384 ): 173 – 178 . Google Scholar Crossref Search ADS PubMed WorldCat Maga A M , Navarro N, Cunningham M L, Cox T C, 2015   Quantitative trait loci affecting the 3D skull shape and size in mouse and prioritization of candidate genes in-silico. Front. Physiol. 9(92) : doi: 10.3389/fphys.2015.00092 . Google Scholar OpenURL Placeholder Text WorldCat Matthews M K Jr, Adler N T, 1978   Systematic interrelationship of mating, vaginal plug position, and sperm transport in the rat. Physiol. Behav. 20 ( 3 ): 303 – 309 . Google Scholar Crossref Search ADS PubMed WorldCat McNeil C L , Bain C L, Macdonald S J, 2011   Multiple Quantitative Trait Loci Influence the Shape of a Male-Specific Genital Structure in Drosophila melanogaster. G3: Genes, Genomes . Genetics 1 ( 5 ): 343 – 351 . Google Scholar OpenURL Placeholder Text WorldCat Mitteroecker P , Gunz P, 2009   Advances in geometric morphometrics. Evol. Biol. 36 : 235 – 247 . Google Scholar Crossref Search ADS WorldCat Miyagawa S , Moon A, Haraguchi R, Inoue C, Harada M et al. , 2009   Dosage-dependent hedgehog signals integrated with Wnt/β-catenin signaling regulate external genitalia formation as an appendicular program. Development 136 ( 23 ): 3969 – 3978 . Google Scholar Crossref Search ADS PubMed WorldCat Murakami R , Mizuno T, 1984   Histogenesis of the Os Penis and Os Clitoridis in Rats. Dev. Growth Differ. 26 ( 5 ): 419 – 426 . Google Scholar Crossref Search ADS WorldCat Nikolskiy I , Conrad D F, Chun S, Fay J C, Cheverud J M et al. , 2015   Using whole-genome sequences of the LG/J and SM/J inbred mouse strains to prioritize quantitative trait genes and nucleotides. BMC Genomics 16 ( 1 ): 415 . Google Scholar Crossref Search ADS PubMed WorldCat Norgard E A , Lawson H A, Pletscher L S, Wang B, Brooks V R et al. , 2011   Genetic factors and diet affect long-bone length in the F34 LG, SM advanced intercross. Mamm. Genome 22 ( 3–4 ): 178 – 196 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Parker C , Cheng R, Sokoloff G, Lim J, Skol A et al. , 2011   Fine-mapping alleles for body weight in LG/J × SM/J F2 and F34 advanced intercross lines. Mamm. Genome 22 ( 9–10 ): 563 – 571 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Parker C C , Cheng R, Sokoloff G, Palmer A A, 2012   Genome‐wide association for methamphetamine sensitivity in an advanced intercross mouse line. Genes Brain Behav. 11 ( 1 ): 52 – 61 . Google Scholar Crossref Search ADS PubMed WorldCat Parker C C , Carbonetto P, Sokoloff G, Park Y J, Abney M et al. , 2014   High-Resolution Genetic Mapping of Complex Traits from a Combined Analysis of F2 and Advanced Intercross Mice. Genetics 198 ( 1 ): 103 – 116 . Google Scholar Crossref Search ADS PubMed WorldCat Patterson B D , 1983   Baculum-body size relationships as evidence for a selective continuum on bacular morphology. J. Mammal. 64 ( 3 ): 496 – 499 . Google Scholar Crossref Search ADS WorldCat Patterson B D , Thaeler C S Jr., 1982   The mammalian baculum: hypotheses on the nature of bacular variability. J. Mammal. 63 ( 1 ): 1 – 15 . Google Scholar Crossref Search ADS WorldCat Peichel C L , Prabhakaran B, Vogt T F, 1997   The mouse Ulnaless mutation deregulates posterior HoxD gene expression and alters appendicular patterning. Development 124 ( 18 ): 3481 – 3492 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Peirce J L , Lu L, Gu J, Silver L M, Williams R W, 2004   A new set of BXD recombinant inbred lines from advanced intercross populations in mice. BMC Genet. 5 : 7 . Google Scholar Crossref Search ADS PubMed WorldCat Peluffo, A.E., I. Nuez, V. Debat, R. Savisaar, D.L. Stern et al., 2015 A Major Locus Controls a Genital Shape Difference Involved in Reproductive Isolation Between Drosophila yakuba and Drosophila santomea. G3 (Bethesda) 5(12):2893–2901. Perriton C L , Powles N, Chiang C, Maconochie M K, Cohn M J, 2002   Sonic hedgehog Signaling from the Urethral Epithelium Controls External Genital Development. Dev. Biol. 247 ( 1 ): 26 – 46 . Google Scholar Crossref Search ADS PubMed WorldCat Perry J C , Rowe L, 2012   Sexual conflict and antagonistic coevolution across water strider populations. Evolution 66 : 544 – 557 . Google Scholar Crossref Search ADS PubMed WorldCat Pulvers J N , Bryk J, Fish J L, Wilsch-Bräuninger M, Arai Y et al. , 2010   Mutations in mouse Aspm (abnormal spindle-like microcephaly associated) cause not only microcephaly but also major defects in the germline. Proc. Natl. Acad. Sci. USA 107 ( 38 ): 16595 – 16600 . Google Scholar Crossref Search ADS WorldCat R Core Team, 2014 R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Rai M F , Schmidt E J, Hashimoto S, Cheverud J M, Sandell L J, 2015   Genetic loci that regulate ectopic calcification in response to knee trauma in LG/J by SM/J advanced intercross mice. J. Orthop. Res. 33 ( 10 ): 1412 – 1423 . Google Scholar Crossref Search ADS PubMed WorldCat Ramm S , Khoo L, Stockley P, 2010   Sexual selection and the rodent baculum: an intraspecific study in the house mouse (Mus musculus domesticus). Genetica 138 ( 1 ): 129 – 137 . Google Scholar Crossref Search ADS PubMed WorldCat Ramm S A , 2007   Sexual selection and genital evolution in mammals: a phylogenetic analysis of baculum length. Am. Nat. 169 ( 3 ): 360 – 369 . Google Scholar Crossref Search ADS PubMed WorldCat Ripley B D , 2002   Modern applied statistics with S , Springer-Verlag, New York . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rodriguez E , Weiss D A, Yang J H, Menshenina J, Ferretti M et al. , 2011   New insights on the morphology of adult mouse penis. Biol. Reprod. 85 ( 6 ): 1216 – 1221 . Google Scholar Crossref Search ADS PubMed WorldCat Rohlf, F. J., and F. L. Bookstein, 1990 Proceedings of the Michigan Morphometrics Workshop. Ann Arbor: University of Michigan Museum of Zoology. Rönn J , Katvala M, Arnqvist G, 2007   Coevolution between harmful male genetalia and female resistance in seed beetles. Proc. Natl. Acad. Sci. USA 104 : 10921 – 10925 . Google Scholar Crossref Search ADS WorldCat Rowe L , Arnqvist G, 2012   Sexual selection and the evolution of genital shape and complexity in water striders. Evolution 66 ( 1 ): 40 – 54 . Google Scholar Crossref Search ADS PubMed WorldCat Sanger T J , Gredler M L, Cohn M J, 2015   Resurrecting embryos of the tuatara, Sphenodon punctatus, to resolve vertebrate phallus evolution. Biol. Lett. 11 ( 10 ): 20150694. Google Scholar OpenURL Placeholder Text WorldCat Sasabe M , Takami Y, Sota T, 2007   The genetic basis of interspecific differences in genital morphology of closely related carabid beetles. Heredity 98 ( 6 ): 385 – 391 . Google Scholar Crossref Search ADS PubMed WorldCat Schäfer M A , Routtu J, Vieira J, Hoikkala A, Ritchie M G et al. , 2011   Multiple quantitative trait loci influence intra-specific variation in genital morphology between phylogenetically distinct lines of Drosophila montana. J. Evol. Biol. 24 ( 9 ): 1879 – 1886 . Google Scholar Crossref Search ADS PubMed WorldCat Schulte-Hostedde A I , Bowman J, Middel K R, 2011   Allometry of the baculum and sexual size dimorphism in American martens and fishers (Mammalia: Mustelidae). Biol. J. Linn. Soc. Lond. 104 ( 4 ): 955 – 963 . Google Scholar Crossref Search ADS WorldCat Scimè A , Grenier G, Huh M S, Gillespie M A, Bevilacqua L et al. , 2005   Rb and p107 regulate preadipocyte differentiation into white vs. brown fat through repression of PGC-1α. Cell Metab. 2 ( 5 ): 283 – 295 . Google Scholar Crossref Search ADS PubMed WorldCat Seifert A W , Yamaguchi T, Cohn M J, 2009   Functional and phylogenetic analysis shows that Fgf8 is a marker of genital induction in mammals but is not required for external genital development. Development 136 ( 15 ): 2643 – 2651 . Google Scholar Crossref Search ADS PubMed WorldCat Shivtiel S , Kollet O, Lapid K, Schajnovitz A, Goichberg P et al. , 2008   CD45 regulates retention, motility, and numbers of hematopoietic progenitors, and affects osteoclast remodeling of metaphyseal trabecules. J. Exp. Med. 205 ( 10 ): 2381 – 2395 . Google Scholar Crossref Search ADS PubMed WorldCat Simmons L W , Firman R C, 2014   Experimental evidence for the evolution of the mammalian baculum by sexual selection. Evolution 68 ( 1 ): 276 – 283 . Google Scholar Crossref Search ADS PubMed WorldCat Slice D E , 2007   Geometrics morphometrics. Annu. Rev. Anthropol. 36 : 261 – 281 . Google Scholar Crossref Search ADS WorldCat Stockley P , Ramm S A, Sherborne A L, Thom M D, Paterson S et al. , 2013   Baculum morphology predicts reproductive success of male house mice under sexual selection. BMC Biol. 11 : 66 . Google Scholar Crossref Search ADS PubMed WorldCat Takahara B , Takahashi K H, 2015   Genome-Wide Association Study on Male Genital Shape and Size in Drosophila melanogaster. PLoS One 10 ( 7 ): e0132846 . Google Scholar Crossref Search ADS PubMed WorldCat Tanaka K M , Hopfen C, Herbert M R, Schlötterer C, Stern D L et al. , 2015   Genetic Architecture and Functional Characterization of Genes Underlying the Rapid Diversification of Male External Genitalia Between Drosophila simulans and Drosophila mauritiana. Genetics 200 ( 1 ): 357 – 369 . Google Scholar Crossref Search ADS PubMed WorldCat Tasikas D , Fairn E, Laurence S, Schulte-Hostedde A, 2009   Baculum variation and allometry in the muskrat (Ondatra zibethicus): a case for sexual selection. Evol. Ecol. 23 ( 2 ): 223 – 232 . Google Scholar Crossref Search ADS WorldCat Taylor B A , Heiniger H J, Meier H, 1973   Genetic analysis of resistance to cadmium-induced testicular damage in mice . Proc. Soc. Exp. Biol. Med. 143 : 629 – 633 . Google Scholar Crossref Search ADS PubMed WorldCat Taylor B A , Wnek C, Kotlus B S, Roemer N, MacTaggart T et al. , 1999   Genotyping new BXD recombinant inbred mouse strains and comparison of BXD and consensus maps. Mamm. Genome 10 ( 4 ): 335 – 348 . Google Scholar Crossref Search ADS PubMed WorldCat Thomas O , 1915   XXXIV.—The penis-bone, or “Baculum,” as a guide to the classification of certuin squirrels. J. Nat. Hist. 15 ( 88 ): 383 – 387 . Google Scholar OpenURL Placeholder Text WorldCat Toner J P , Adler N T, 1986   The pre-ejaculatory behavior of male and female rats affects the number of sperm in the vagina and uterus. Physiol. Behav. 36 ( 2 ): 363 – 367 . Google Scholar Crossref Search ADS PubMed WorldCat Toner J P , Attas A I, Adler N T, 1987   Transcervical sperm transport in the rat: the roles of pre-ejaculatory behavior and copulatory plug fit. Physiol. Behav. 39 ( 3 ): 371 – 375 . Google Scholar Crossref Search ADS PubMed WorldCat Trapnell C , Pachter L, Salzberg S L, 2009   TopHat: discovering splice junctions with RNA-Seq. Bioinformatics 25 ( 9 ): 1105 – 1111 . Google Scholar Crossref Search ADS PubMed WorldCat True J R , Liu J, Stam L F, Zeng Z-B, Laurie C C, 1997   Quantitative genetic analysis of divergence in male secondary sexual traits between Drosophila simulans and Drosophila mauritiana. Evolution 51(3) : 816 – 832 . Google Scholar OpenURL Placeholder Text WorldCat Tschopp P , Sherratt E, Sanger T J, Groner A C, Aspiras A C et al. , 2014   A relative shift in cloacal location repositions external genitalia in amniote evolution. Nature ; advance online publication . Google Scholar OpenURL Placeholder Text WorldCat Wang, X., A.K. Pandey, M.K. Mulligan, E.G. Williams, K. Mozhui et al., 2016  Joint mouse–human phenome-wide association to test gene function and disease risk. Nat. Commun. DOI: 10.1038/ncomms10464. Wang Y , Tan J, Sutton-Smith M, Ditto D, Panico M et al. , 2001   Modeling human congenital disorder of glycosylation type IIa in the mouse: conservation of asparagine-linked glycan-dependent functions in mammalian physiology and insights into disease pathogenesis. Glycobiology 11 ( 12 ): 1051 – 1070 . Google Scholar Crossref Search ADS PubMed WorldCat Wang Y , Wan C, Deng L, Liu X, Cao X et al. , 2007   The hypoxia-inducible factor α pathway couples angiogenesis to osteogenesis during skeletal development. J. Clin. Invest. 117 ( 6 ): 1616 . Google Scholar Crossref Search ADS PubMed WorldCat Warot X , Fromental-Ramain C, Fraulob V, Chambon P, Dolle P, 1997   Gene dosage-dependent effects of the Hoxa-13 and Hoxd-13 mutations on morphogenesis of the terminal parts of the digestive and urogenital tracts. Development 124 ( 23 ): 4781 – 4791 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Waterston R H , Lindblad-Toh K, Birney E, Rogers J, Abril J F et al. , 2002   Initial sequencing and comparative analysis of the mouse genome. Nature 420 ( 6915 ): 520 – 562 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Wen J , Chiang Y J, Gao C, Xue H, Xu J et al. , 2010   Loss of Dact1 disrupts planar cell polarity signaling by altering dishevelled activity and leads to posterior malformation in mice. J. Biol. Chem. 285 ( 14 ): 11023 – 11030 . Google Scholar Crossref Search ADS PubMed WorldCat Yamaguchi T P , Bradley A, McMahon A P, Jones S, 1999   A Wnt5a pathway underlies outgrowth of multiple structures in the vertebrate embryo. Development 126 ( 6 ): 1211 – 1223 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Yang H , Wang J R, Didion J P, Buus R J, Bell T A et al. , 2011   Subspecific origin and haplotype diversity in the laboratory mouse. Nat. Genet. 43 ( 7 ): 648 – 655 . Google Scholar Crossref Search ADS PubMed WorldCat Zákány J , Duboule D, 1999   Hox genes in digit development and evolution. Cell Tissue Res. 296 ( 1 ): 19 – 25 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Zeng Z B , Liu J, Stam L F, Kao C H, Mercer J M et al. , 2000   Genetic architecture of a morphological shape difference between two Drosophila species. Genetics 154 ( 1 ): 299 – 310 . Google Scholar PubMed OpenURL Placeholder Text WorldCat Zheng H F , Tobias J H, Duncan E, Evans D M, Eriksson J et al. , 2012   WNT16 influences bone mineral density, cortical bone thickness, bone strength, and osteoporotic fracture risk. PLoS Genet. 8 ( 7 ): e1002745 . Google Scholar Crossref Search ADS PubMed WorldCat © 2016 Schultz et al. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) © 2016 Schultz et al. TI - The Genetic Basis of Baculum Size and Shape Variation in Mice JF - "G3: Genes, Genomes, Genetics" DO - 10.1534/g3.116.027888 DA - 2016-05-01 UR - https://www.deepdyve.com/lp/oxford-university-press/the-genetic-basis-of-baculum-size-and-shape-variation-in-mice-spmglImLYO SP - 1141 EP - 1151 VL - 6 IS - 5 DP - DeepDyve ER -