Research ArticleEVOLUTIONARY BIOLOGY

Genetic basis of species-specific genitalia reveals role in species diversification

See allHide authors and affiliations

Science Advances  26 Jun 2019:
Vol. 5, no. 6, eaav9939
DOI: 10.1126/sciadv.aav9939

Abstract

The diversity of genital morphology among closely related animals with internal fertilization is well known, but the genetic backgrounds are unclear. Here, we show that, in Carabus (Ohomopterus) beetles showing correlated evolution of male and female genital parts, only a few major quantitative trait loci (QTLs) determine differences in genital dimensions between sister species, and sequence divergence is pronounced in the genomic regions containing genital QTLs. The major QTLs for male and female genital dimensions reside in different locations within the same linkage group, implying that coevolution between the sexes is only loosely constrained and can respond to sexually antagonistic selection. The same genomic regions containing the major QTLs show elevated divergence between three pairs of parapatric species with marked differences in genital parts. Our study demonstrates that species diversification can follow coevolution of genitalia between the sexes, even without tight linkage of loci affecting male and female genital dimensions.

INTRODUCTION

The genital morphology of animals with internal fertilization shows diversity among species, often more so than their external morphology (1). The evolution of genital morphology is thought to be chiefly promoted by sexual selection (1, 2). Many studies have explored the causes and consequences of genital evolution (3). However, the genetic basis of species-specific genitalia remains unclear. In Drosophila flies, the genetic backgrounds of the species-specific shape and size of male external genitalia among allopatric species have been studied (47). However, a more comprehensive study of genital genetics in the context of speciation is still lacking.

Genital genetics has rarely been studied in animals other than Drosophila, with one exception being ground beetles of the genus Carabus, subgenus Ohomopterus (8, 9). This group has corresponding male and female genital parts, the copulatory piece in the male endophallus and the vaginal appendix in the female vagina; the former is inserted into the latter, resulting in a secure genital coupling during insemination (Fig. 1A) (1012). The shape and size of these genital parts match and are species specific, and their lengths show correlated evolution between the sexes toward exaggeration (Fig. 1B), which is especially accelerated in the Carabus iwawakianus–Carabus insulicola group (Fig. 1C) (13). Sexual selection, including sperm competition and sexual conflicts, may be involved in the correlated evolution of Ohomopterus genitalia (1416), and the resultant differences in genital morphology incur large costs in interspecific mating and contribute to mechanical reproductive isolation (17). The external morphology, along with the ecology of Ohomopterus beetles, is highly homogeneous, except for body size, in which divergence also causes mechanical reproductive isolation and promotes species coexistence (18). Thus, the species diversity of Ohomopterus largely depends on the divergent evolution of genital morphology and body size (13).

Fig. 1 Genital structure, diversity, and phylogeny of Carabus (Ohomopterus) species.

(A) Structure and coupling of male and female genitalia. (B) Correlated divergence of the male copulatory piece and female vaginal appendix length (VAL) [data from (9)]. Both lengths were divided by the pronotum width to correct for body size differences, multiplied by 103, and the common logarithm was taken. (C) Phylogeny and divergence time based on RAD sequences. Photographs of males and drawings of the male genital copulatory piece are shown at the right. (Photo credit: Teiji Sota, Kyoto University).

The genetic background for species-specific genital morphology can be investigated with genomic analyses of closely related species. Despite their diverged genitalia, closely related parapatric species can hybridize (13). A well-known example is the hybridization of C. iwawakianus and Carabus maiyasanus, which exhibit contrasting dimensions in the key genital parts (C. maiyasanus has an elongated copulatory piece and vaginal appendix, whereas C. iwawakianus has short ones) but have hybrid zones containing various intermediate genital morphologies (19). Laboratory interspecific crossing of these species has been used to study the genetic background of species-specific genital morphology (8, 9).

In this study, we explored the evolutionary history of genital diversification and speciation by reconstructing phylogenetic relationships among Ohomopterus species based on restriction site–associated DNA (RAD) sequences. Second, we obtained a reference whole-genome sequence of Carabus uenoi, which is closely related to both C. iwawakianus and C. maiyasanus, and performed quantitative trait locus (QTL) mapping of male and female genital dimensions (length and width of the copulatory piece and vaginal appendix) using F2 hybrids resulting from interspecific hybridization between C. iwawakianus and C. maiyasanus. Then, focusing on genomic sequence scaffolds harboring genital QTLs, we searched for single-nucleotide polymorphisms (SNPs) and genes related to species-specific genital morphology. We also compared gene expression in larvae and pupae between the two species to identify genes involved in differentiation of species-specific morphology. We found only a few QTLs for genital dimensions, and the major QTLs for the length of the copulatory piece and vaginal appendix resided within different regions of the same linkage group. Genomic regions involving genital QTLs diverged more than other regions among the three species. We also analyzed the sequence divergence of other pairs of species with diverged genital morphology and found that the same region showed high divergence in parallel with the C. iwawakianus and C. maiyasanus pair.

RESULTS

Phylogeny and genomic sequence of Ohomopterus

The phylogenetic tree based on RAD sequences (Fig. 1C and table S1) resolved the relationships among species distinguished by morphology, especially of the genitalia, with the exception of those among C. iwawakianus, C. maiyasanus, and C. uenoi, which are parapatric or sympatric (Fig. 2A). The former two species are parapatric and have hybrid zones; introgressive hybridization likely has occurred repeatedly (13, 19). The time of the most recent common ancestor (tMRCA) of Ohomopterus was estimated to be approximately 2.6 million years (Ma) ago (the beginning of the Pleistocene), and the speciation rate [ln(no. of species)/divergence time] for the 17 species was 1.09 per Ma. The tMRCA of the C. iwawakianusC. insulicola group, which contains species with highly exaggerated genital parts, was 1.3 Ma ago, and the speciation rate was higher (1.50 per Ma) than in the other groups (Carabus albrechti group, 0.7 per Ma; Carabus japonicus group, 0.33 per Ma; Carabus dehaanii group, 0.64 per Ma), indicating an elevated speciation rate in the lineage exhibiting a marked diversification of genital morphology.

Fig. 2 Three species used for genome sequencing and interspecific hybridization, F2 phenotypes, and the results of QTL mapping.

(A) Distribution ranges and sampling locations. C. maiyasanus, 1 to 4, 8, and 9; C. iwawakianus, 5 to 7, 10, and 11; C. uenoi, 6. For each species, pictures of the male genitalia with an everted endophallus (left) and female genitalia (right) are shown. (B) Measurement of the genital dimensions of the copulatory pieces and vaginal appendix in C. maiyasanus and C. iwawakianus. (C and D) Genital dimensions for F2 males (C) and F2 females (D) (closed circles) and those of the parental species (C. iwawakianus, squares; C. maiyasanus, triangles; eight examples for each). (E and F) LOD (logarithm of the odds) scores of QTL mapping analysis of the genital dimensions of males (E) and females (F). Broken lines indicate the thresholds for significant QTLs. cM, centimorgan. (Photo credit: Teiji Sota, Kyoto University).

To explore the genomic basis of species-specific genital morphology, we focused on the marked divergence among C. iwawakianus, C. maiyasanus, and C. uenoi, the last of which shows the most exaggerated genitalia of all Ohomopterus species (Fig. 2A). We used F2 individuals from intercrossing of C. iwawakianus and C. maiyasanus (9) for revised QTL analysis in this study. Furthermore, we selected C. uenoi for the whole-genome sequencing because this species is confined to a small area, and its genome was expected to show relatively low heterozygozity, making it suitable as a reference (table S2). The draft genomic sequence of C. uenoi was 186,290,795 base pair (bp) in total length (48,829 scaffolds; N, 6.1%) and 3815 bp in mean sequence length; scaffold N50 was 1,370,541 bp, with a maximum scaffold length of 7,513,851 bp and GC (guanine-cytosine) content of 34.8%. Using the total evidence approach, we obtained 15,607 protein-coding genes, of which 70% were assigned to known genes in the Tribolium castaneum genome. The completeness of the genome assembly and predicted gene set was assessed using the Benchmarking Universal Single-Copy Orthologs (BUSCO) (20). Among the 2442 Endopterygota (holometabolous insects) universal orthologs, 88% (2150 genes; 2137 single copy, 13 duplicated) were complete, 5.4% (133 genes) were fragmented, and 6.6% (159 genes) were missing in the predicted protein set. The completeness was 96% (2351 genes; 2338 single copy, 13 duplicated) when BUSCO was run on the assembled genome.

QTLs for genital dimensions

The F2 individuals of the C. iwawakianus × C. maiyasanus cross showed various intermediate values with respect to the genital dimensions (Fig. 2, C and D). In males, copulatory piece length (CPL) was correlated with copulatory piece width (CPW) (r = −0.4142, n = 54, P = 0.0018; Pearson’s correlation). However, in females, vaginal appendix length (VAL) and vaginal appendix width (VAW) were not correlated with each other (r = −0.1556, n = 24, P = 0.4680; Pearson’s correlation). Thus, the covariation pattern between the length and the width of genital parts within the sexes differed between the sexes.

Using RAD sequence data from F2 individuals, we obtained 20 linkage groups, including one for the X chromosome (fig. S1). In QTL mapping for male genital dimensions, we found four significant QTLs for CPL and three for CPW in three autosomal linkage groups (LG1, LG7, and LG19) and the X chromosome (Fig. 2E and Table 1). These loci showed large additive effects. In particular, the QTLs for CPL and CPW were located at almost the same position on LG1 and explained the highest proportion of the variance (48.6% for CPL; 30.5% for CPW). The sum of the proportions of the variance explained by the detected QTLs was 119.6% for CPL and 61.0% for CPW. Regarding female genital dimensions, two significant QTLs for VAL were detected on LG1, explaining 60.3 and 55.4% of the variance (Fig. 2F and Table 1). However, their positions were not close to the QTLs for CPL and CPW. The VAL QTLs showed large dominance effects in addition to additive effects. Three QTLs for VAW were detected on LG2, LG9, and LG17. These QTLs also accounted for a relatively high proportion of the variance (42.4 to 74.2%), although this was apparently overestimated probably because of the small sample size of F2 individuals (Beavis effect) (21). Thus, the linkage group LG1 harbored QTLs with large effects on the lengths of the corresponding male and female genital parts, as well as on the width of the male genital part.

Table 1 QTLs for male and female genital dimensions.

Effect: a, additive; d, dominance. R2, proportion of variance explained.

View this table:

Genomic regions showing monophyly of species

We sequenced the whole genomes of C. iwawakianus and C. maiyasanus individuals sampled far from hybrid zones and searched for species-monophyly regions in 10 scaffolds with genital QTLs (total, 22,798,092 bp) and 11 without genital QTLs (total, 46,575,664 bp) using ARGweaver version 0.8.1 (22). Species-monophyly regions are genomic regions that contain SNPs sorted by species and show monophyly of two species in the local genealogy. Species-monophyly regions were detected in five scaffolds with genital QTLs (Table 2 and table S3); the total length was 5744 bp, which was only 0.03% of the total length of the scaffolds. Most species-monophyly regions (5665 bp, 98%) were located within genes or flanking regions of genes; only 2% of the regions were within noncoding sequences. Scaffolds without genital QTLs examined also contained species-monophyly regions, spanning even smaller proportion (0.007%) of the total scaffold length (Table 2).

Table 2 Species-monophyly regions in genome scaffolds.

Numbers of species-monophyly regions and genes within 2500 bp of segregating SNPs, as determined by ARGweaver analysis.

View this table:

By annotating 27 genes found within 2500 bp upstream and downstream of the species-monophyly regions in scaffolds with genital QTLs, we identified 21 genes (table S4). Of these, NOA1 (nitric oxide–associated protein 1) was located at the LOD (logarithm of the odds) peak position of the QTL for CPL on scaffold 25. Enrichment analysis of the 21 genes suggested that the species-monophyly regions characteristically contain cuticle pigmentation genes; these were pinstripe (pns) on scaffold 237 and Kinesin-like protein at 61F (Klp61F) on scaffold 406. Fixed SNPs between the two species were located in the intergenic region of pns and M6, as well as in an intron region of M6 and in an intron region of Klp61F. On scaffold 237, fixed SNPs were also found in and around the intergenic region of zinc finger FYVE domain–containing protein 1–like (ZFYVE1) and transmembrane protein 53 (TMEM53). Comparing the 21 genes with candidate genes for Drosophila genital morphology (6, 23), we found that one gene (CG7971) on scaffold 25 and three genes (cubilin, masquerade, and Klp61F) on scaffold 406 of the Ohomopterus genes are differentially expressed in the male genital discs of Drosophila mauritiana and Drosophila sechellia (table S4) (6).

Differentially expressed genes

To detect genes involved in the development of species-specific genital morphologies, we compared the gene expression of a total of 14,820 loci in third-instar larvae (5 days prepupae) and pupae (2 days after pupation) of C. iwawakianus and C. maiyasanus. By comparing gene expression between the two species, we found 873 and 676 differentially expressed genes (DEGs) for male and female larvae, respectively, and 408 and 1059 DEGs for male and female pupae, respectively. The DEGs were enriched significantly for one and two functional categories in male larvae and female pupae, respectively (table S5), but these functions appeared not to be related to genital development. Although not significant, three functional categories related to pigmentation showed relatively high enrichment rates in male larvae, compared to one in female larvae and none in pupae (table S5).

We evaluated the expression of developmental genes potentially involved in genital and appendage development and determination of body part size (table S6). We found that both Abdominal-B and doublesex (dsx), which indicated sex-specific genital development in Drosophila (24), were similarly expressed in both sexes of the larval and pupal samples in both species. Among other genes potentially involved in genital disc development and appendage patterning, the within-disc morphogen decapentaplegic (dpp) and the transcription factor rotund (rn) were up-regulated in male pupae and larvae of C. maiyasanus, respectively, and the transcription factor escargot was down-regulated in female larvae. We did not find DEGs for genes involved in the Fat/Hippo and insulin-like growth factor signaling pathways with the exception of 4E binding protein, which was down-regulated in female pupae. We also evaluated the expression of 27 genes associated with the species-monophyly regions detected by the ARGweaver analysis and found that six genes, including ZFYVE1 (up-regulated in male larvae) and TMEM53 (down-regulated in female pupae) on scaffold 237, were differentially expressed (table S4).

Genomic divergence patterns between species

To characterize the pattern of genomic sequence divergence between C. iwawakianus and C. maiyasanus, we obtained 10 kilo–base pair (kbp) average Fst between C. iwawakianus and C. maiyasanus, as well as between both of these and C. uenoi (Fig. 3, A to C). The 10-kbp average Fst between C. iwawakianus and C. maiyasanus (parapatric pair with moderate genital differences) was much smaller than that between C. iwawakianus and C. uenoi (sympatric pair with large genital differences) and that between C. maiyasanus and C. uenoi (allopatric pair with large genital differences). There were several regions with elevated Fst, which were generally common in all three comparisons, and some of these regions were located in scaffolds with genital QTLs (Fig. 3, A to C). The mean 10-kbp Fst was greater for scaffolds with QTLs than without genital QTLs in each species pair, and it was much smaller for the C. iwawakianusC. maiyasanus pair than in the other comparisons with C. uenoi (Fig. 3D and table S7). The 10-kbp Fst was negatively correlated with the local effective population size (Ne) of the 10-kbp sections, and this trend was most prominent for Fst between C. iwawakianus and C. uenoi, which are sympatric in the limited locality of the latter species but show large differences in genital morphology that would prevent hybridization (Fig. 3E and table S7).

Fig. 3 Genomic sequence divergence among three species in terms of the 10-kbp average Fst.

(A to C) Divergence of genomic sequences (Fst at 10-kbp intervals) among C. iwawakianus, C. maiyasanus, and C. uenoi. Scaffolds longer than 300 kbp are shown. In (A), scaffolds with genital QTLs are colored. (D) Comparison of Fst between scaffolds with and without genital QTLs and among species pairs. (E) Plot of Fst against the local effective population size of the 10-kbp sections (in log10 scale) for the C. iwawakianusC. maiyasanus (green), C. maiyasanusC. uenoi (blue), and C. iwawakianusC. uenoi (pink) pairs. Solid, broken, and dotted lines indicate regression lines for the C. iwawakianusC. maiyasanus, C. maiyasanusC. uenoi, and C. iwawakianusC. uenoi pairs, respectively.

Parallel sequence divergences associated with divergence of genital morphology

We compared the RAD sequence divergence of five scaffolds, in which candidate genes were found by QTL and ARGweaver analyses, among three closely related species pairs with markedly divergent genital part lengths (Fig. 4, A to E). We selected two pairs of parapatric species with contrasting genital dimensions, the Carabus komiyaiCarabus arrowianus and Carabus esakiiC. insulicola pairs, for comparison with the C. iwawakianusC. maiyasanus pair. SNPs were extracted from RAD sequences using C. uenoi scaffolds as the reference. We used the weighted mean Fst at 50-kbp intervals as a measure of sequence divergence. The Fst values for intervals with QTLs or species-monophyly regions showed coincident elevation in all or two of the tree species pairs. For example, in scaffold 237 (Fig. 4C, top), Fst was highest at the interval of 2150 to 2200 kbp, which contains the pns locus. This region also showed high Fst values for the C. arrowianusC. komiyai and C. esakiiC. insulicola pairs (Fig. 4C, bottom). Similarly, the species-monophyly region in the middle of scaffold 25, which harbored a QTL for CPL, showed high Fst values for all three species pairs (Fig. 4D).

Fig. 4 Sequence divergence based on RAD sequences in five scaffolds containing genital QTL candidate genes between parapatric species with diverged genital morphology.

(A) Distribution of three pairs of species and C. uenoi. (B) Neighbor-joining tree of species based on the Fst values between species using the RAD sequence data of scaffold 237. Illustrations of male copulatory pieces are shown. (C to G) For each scaffold, top panels show weighted mean Fst of 50-kbp genomic sequence windows (black circles) and the posterior probability of segregation between C. iwawakianus and C. maiyasanus (orange bars); green lines indicate the average Fst of the scaffolds; pale blue arrows indicate genital QTL positions; and pink arrows with numbers indicate intervals that contain candidate genes. The positions of candidate gene loci are indicated at the top (blue bars) with gene names mentioned in the text. Bottom panels show the correlation of Fst values from the C. iwawakianusC. maiyasanus pair (x axis) and C. arrowianusC. komiyai (blue) or between C. esakiiC. insulicola (red) pairs (y axis). Pale blue and pink arrows correspond to those in the top panels.

DISCUSSION

We explored the genetic basis of species-specific genital morphology using various complementary approaches. Although the results of our genital QTL mapping may not be sufficient enough because of limited sample sizes, we detected diverged genomic regions in terms of Fst values and species-monophyly regions, the latter being a stricter indicator of interspecific differentiation. Although the diverged genomic regions are not necessarily related to genital differences, we found correspondence among the evidences from Fst, species monophyly, and QTL analyses in some genomic regions, which should constitute the genetic basis of species-specific genitalia. However, we have yet to specify the genes and nucleotide substitutions responsible for interspecific differences in genital morphology. Therefore, additional analyses of candidate genomic regions and functions of candidate genes are needed in future studies.

The functional parts of male and female internal genitalia of Ohomopterus beetles show a polygenic inheritance pattern (8), but we found only a few genomic regions (QTLs) with major effects. For the copulatory piece of the male genitalia, the length (CPL) QTL with the largest effect resided on LG1, and the width (CPW) QTL resided at almost the same position. Together with the fact that these dimensions were negatively correlated in the F2 population, the length and width appear to be pleiotropically controlled by the same locus. In the corresponding female genital part, the vaginal appendix, the length (VAL) and width (VAW) were not correlated in F2, and the VAL and VAW QTLs were separated. QTLs for male and female genitalia existed in different genomic regions, although the VAL QTL was found in the same linkage group as the major CPL QTL. Thus, the sizes of the corresponding male and female genital parts appear to be almost independently controlled.

We have obtained candidate genes for species-specific genital morphology by identifying genes associated with genomic regions showing species monophyly in the scaffolds with genital QTLs. Although the roles of these candidate genes in genital formation are unknown, the candidate genes particularly contained pigmentation genes (pns and Klp61F). These genes affect abdominal pigmentation in Drosophila melanogaster adults (25). Klp61F was differentially expressed in male genital discs of D. mauritiana and D. sechellia (6) and may be related to the species-specific morphology of a cuticular projection of the genital arch (the posterior lobe) in the male external genitalia. Variations in the shape and size of sclerotized parts constitute the most obvious differences in male insect genitalia among species, and sclerotized genital parts are usually pigmented. Therefore, some common genes that regulate the formation of sclerotized parts, including pigmentation genes, may be involved in the divergence of male insect genitalia. The copulatory piece, a sclerotized part of the endophallus (Fig. 2A), is an innovative character of a group of the genus Carabus including Ohomopterus, called Digitulati, and is considered homologous to the pigmented patch on the membranous endophallus of another group, the subgenus Morphocarabus, which lacks a copulatory piece (26). Therefore, differential expression of pigmentation genes may be involved in the evolution and subsequent diversification of the copulatory piece.

In the analysis of DEGs, we found that the morphogen dpp and the transcription factor rn were up-regulated in C. maiyasanus male pupae and larvae, respectively. Higher expression of these genes may be involved in elongation of the copulatory piece in C. maiyasanus. The rn transcription factor particularly is necessary for the development of various appendages and wings in Drosophila, and we found that the rn gene knockdown via larval RNA interference (RNAi) suppressed the formation of the copulatory piece in C. maiyasanus (27). An elevated rn expression may be required for the extension of the copulatory piece. On the other hand, rn RNAi did not affect the development of female genitalia (27). This is consistent with the results of QTL mapping, which indicated that the genetic basis of the development of genital parts differs between the sexes.

As we have not provided any evidence for the involvement of aforementioned genes, except rn, in the formation of species-specific genitalia, we must study their functions by knocking down candidate genes via larval RNAi. Larval RNAi is generally effective in beetles including Carabus (27) and is one of the only available methods for assessing the functions of focal genes in nonmodel beetles, which are difficult to breed in large numbers.

Our study provides insight into the coevolution of mating traits between the sexes. Matching of male and female traits for mating is important for reproduction, and this is also true for genitalia. Compared to other body dimensions, genital size is less affected by disturbances because of nutritional conditions (28); thus, matching between the sexes can be maintained developmentally. In our study, the major QTLs for CPL and VAL, which show species-specific matching, were located on the same linkage group but at separated positions. This implies that the evolutionary changes in male and female genital dimensions are only weakly linked. The lack of tight linkage is rather unexpected because matching of reproductive traits between the sexes can be more easily achieved when the genes determining those traits are linked or when the traits are pleiotropically controlled by the same genes. Tight linkages between QTLs (29, 30) and single-gene control (31, 32) of male and female reproductive traits have been demonstrated in insects. Compared to matching between the sexes based on tight linkages of genes or single-gene control, the genital matching based on weakly linked loci appears to be unstable against external disturbances such as interspecific hybridization. However, strong selection against mismatching appears to have been exerted in natural populations because the occurrence of various intermediate genital morphologies is restricted to narrow hybrid zones (17).

One-gene or one-region control of sexual traits may restrict the coevolution of the sexes (33) and can hinder resolution of sexual conflict. However, the lack of a tight genetic link between male and female sexual traits would allow them to have independent responses to selection under conditions of sexual conflict (i.e., interlocus sexual conflict) and may facilitate antagonistic coevolution of sexual traits (34). In Ohomopterus beetles, genital matching between the sexes ensures insemination (11, 12) and is also favored in sperm competition (16). However, this selection alone cannot explain the correlated genital evolution between the sexes toward exaggeration. Because the vaginal appendix is unique to the subgenus Ohomopterus and is not shared by related subgenera with small male copulatory pieces, the acquisition of a vaginal appendix may have triggered the enlargement and diversification of the copulatory piece, followed by correlated evolution between the sexes. The elongation of the copulatory piece may have been promoted by sperm competition (15), but correlated evolution of female genitalia cannot be explained solely by sperm competition. A recent study showed that a longer copulatory piece incurs a cost for females (reduced fecundity), but a longer vaginal appendix counteracts this and reduces the cost (14). Thus, the correlated genital divergence between the sexes in Ohomopterus may have been driven by sexual conflict and may represent sexually antagonistic coevolution (35). The genetic basis of male and female genitalia in Ohomopterus appears to allow a response to antagonistic selection under sexual conflict, despite strong selection pressure for matched coupling parts between the sexes.

Aside from the diverged genomic regions responsible for characteristic differences among closely related species, differences in other genomic regions largely depend on the divergence time, gene flow between species, and background selection process. In the three Ohomopterus species, the local Fst was negatively correlated with the level of allelic fixation within species indicated by the local Ne (Fig. 3E), as was expected from linked selection (36). In general, it is difficult to distinguish regions under species-specific selection (here, divergent selection on genital dimensions) from the effects of background purifying selection (here, those unrelated to genital evolution) because both processes can reduce the local Ne and increase the local Fst (36). In our study, we separated the two effects by combining analyses of genomic divergence with QTL analysis. The divergent regions that we found in the vicinity of the genital QTLs are likely under the species-specific selection process.

The patterns of correlation between local Fst and Ne (Fig. 3E) seem to reflect the different diversification processes of the species pairs. Specifically, genomic divergence in terms of the local Fst appears to have been promoted by both allopatry (C. maiyasanus and C. uenoi) and extreme differences in genital morphology (C. iwawakianus and C. uenoi in sympatry); both could reduce gene flow. The steeper decline in Fst observed in the latter sympatric pair (Fig. 3E) indicates that the genomic regions under strong selection or with tight linkage (hence small local Ne) are more frequently involved in the divergence process driven by selection on genital morphology. In contrast, the parapatric C. iwawakianus and C. maiyasanus showed fairly low Fst values for most of the genomic regions, suggesting extensive gene flow. Gene flow, as indicated by mitochondrial introgression, has been detected in almost all parapatric and sympatric species pairs of Ohomopterus (13), but interspecific differences in genital morphology have been maintained except within hybrid zones. This suggests that species-specific genital morphologies have been maintained by selection against the large mismatch between the sexes (17).

The genitally divergent species in Ohomopterus indicate that speciation or species diversification is possible solely via differentiation of sexually selected traits, although this has been considered unlikely or rare (37, 38). We found parallel divergence of the same genomic regions among species of the C. iwawakianusC. insulicola species group that exhibit divergent genital morphology (Fig. 4). Although the shape and size of genital parts diverge among closely related species, mutations in the same genomic regions may have caused the evolution of genital morphology. The genomic divergence pattern among Ohomopterus beetles demonstrates how divergence of a key trait for reproductive isolation (i.e., genital morphology) may result in species richness despite occasional gene flow due to hybridization, which is a particular case of speciation with gene flow (39, 40). Our study demonstrates that genital matching between the sexes, which is the primary factor in copulatory compatibility, can be achieved and maintained by natural and sexual selection without strong physical linkage or pleiotropy. Natural selection against mismatching appears to be effective despite occasional disturbance by interspecific hybridization at secondary contacts. The antagonistic selection caused by sexual conflict can mediate coevolution between the sexes. In the study of speciation, the role of genital evolution in species diversification is one of the least explored topics (41). Future genomic studies in Ohomopterus beetles should identify the genes and mutations responsible for species-specific male and female genitalia and confirm their pivotal role in maintaining species differences via selection against mismatching between the sexes.

MATERIALS AND METHODS

RAD sequencing

We obtained RAD sequence data from 89 individuals representing all 17 species and major subspecies of Ohomopterus and from an outgroup species Carabus (Isiocarabus) fiduciarius (table S1). The methods of DNA extraction, library construction, and sequencing were the same as described previously (42). The data were processed using pyRAD version 3.0.5 (43). We used a sequence similarity of 70% and minimum number of taxa of 45 for data matrix construction (alignment length, 970,742 bp). Maximum likelihood analysis was performed using RAxML version 8.1.2 (44), and a rapid bootstrap analysis with 100 replications was conducted using a general time-reversible model with gamma-distributed rate variation. For divergence time estimation, we determined the tMRCA of Ohomopterus using the mitochondrial cytochrome oxidase c subunit I (COI) gene sequences (45) and BEAST version 1.8.4 (46) (strict clock model), applying the evolutionary rate of beetle COI (47) (0.0177 per Ma per lineage, SD = 0.0019; normal prior) as a prior. The estimated tMRCA of Ohomopterus COI sequences was 2.6 Ma [95% highest posterior density interval (HPDI), 1.9 to 3.5 Ma]. Note that mitochondrial gene genealogy is affected by introgressive hybridization, and node ages within Ohomopterus cannot be estimated reliably. Therefore, the COI tree was only used for estimating the tMRCA. The maximum likelihood tree of RAD sequences was made ultrametric to obtain the ages of different nodes based on the tMRCA of Ohomopterus (2.6 Ma) using the LSD program version 0.3beta (48). To account for the uncertainty of the Ohomopterus tMRCA, we conducted bootstrap analyses with LSD using the lower and upper limits of 95% HPDI (i.e., 1.9 and 3.5 Ma) and adopted the youngest and oldest ages of the 95% bootstrap confidence limits as the confidence intervals of node ages.

Genomic and mRNA sequences

We conducted whole-genome sequencing for males of C. uenoi (table S2). Total DNA was extracted from the testes of three males using Genomic Tip (Qiagen, Hilden, Germany) and used for Illumina paired-end and mate-pair sequencing with insert lengths of 170 bp, 2 kbp, and 5 kbp on a HiSeq 2000 sequencer (Illumina, San Diego, CA, USA). One lane of a HiSeq 2000 run each was used for the 170-bp paired-end library and the combined 2- and 5-kbp mate-pair libraries. In addition, we sequenced 170-bp paired-end libraries of an additional four C. uenoi, five C. iwawakianus, and six C. maiyasanus males (for details of libraries, see table S2). Library construction and sequencing were performed at Beijing Genomics Institute (Shenzhen, China).

We also obtained transcriptome sequence data for two third-instar larvae of C. uenoi 10 to 11 days after molting. The tissues were fixed in RNAlater and stored in a freezer until RNA extraction. Total RNA was extracted from the abdominal part containing genital part (primordial) and the head to mid-body using a Qiagen RNeasy kit. For each larva, two paired-end libraries with 170-bp insertions were constructed (one each for abdominal and head-mid body) and sequenced using an Illumina HiSeq 2000 at Beijing Genomics Institute, Shenzhen, China (for details of the libraries, see table S2).

Genome assembly, gene prediction, and annotation

The cleaned sequence reads of C. uenoi genome in three libraries (insert sizes: 170 bp, 2 kbp, and 5 kbp) were assembled using Platanus version 1.2.1 (49) with the default parameter settings. Methods of subsequent gene prediction and annotation with the assembled sequences are described in Supplementary Methods. The completeness of the predicted gene set and assembled genome was assessed using the BUSCO pipeline version 3.0.2 (20) with the Endopterygota single-copy orthologs (endopterygota_odb9).

QTL mapping

We reanalyzed F2 adult data (54 males and 24 females) for the C. iwawakianus × C. maiyasanus cross (9) using revised measurements of genital dimensions (table S8) and RAD sequencing for linkage map construction and QTL mapping. RAD sequencing with PstI digestion was performed as described previously (42). In constructing a linkage map, because multiple male and female parents were used, we first determined SNPs on orthologous RAD loci, which were segregated between the male and female parents, and determined genotypes of F2 individuals at the segregating loci. We discriminated SNPs on the X chromosome from those on autosomes based on the criterion that no heterozygote existed in males and no homozygote existed in females. Linkage maps were constructed for autosomes and the X chromosome separately using JoinMap version 4.1 (Kyazma V.B., Wageningen, The Netherlands).

We conducted composite interval mapping using the Windows QTL Cartographer version 2.5_011 (50) for autosomal chromosomes by sex because genital morphology was sex specific. QTL mapping for the X chromosome was performed only on male data because the linkage map was constructed from male data only. We conducted a composite interval mapping (51) with a walk speed of 1 centimorgan (cM; standard background control method with forward and backward regression; control marker number, 5; probability of into/out, 0.1; window size, 10 cM). LOD threshold values at α = 0.05 were obtained by 1000 permutations.

DEGs in larvae and pupae

For both species, we obtained abdominal tissues of third-instar larvae at 5 days after burrowing into the soil for pupation (i.e., prepupae, about 2 days before pupation) and the genital parts of pupae at 2 days after pupation. Development of genital discs was assumed to have started during the prepupal stage, and the copulatory piece and vaginal appendix were assumed to be formed from 2 days after pupation. The tissues of abdominal segments A9 to A11 of larvae and genital tissues of pupae were dissected, fixed in RNAlater, and stored in a freezer until RNA extraction. The sex of a larva was determined by the length difference of the polymerase chain reaction products of dsx gene due to the sex-specific isoforms. The sex of a pupa was determined by the external morphology of the abdominal tip. For each species and sex, three individuals were obtained (total, 24 individuals; table S2). Total RNA was extracted using a Qiagen RNeasy kit. For each sample, paired-end libraries with 170-bp insertions were constructed and sequenced on the Illumina HiSeq 2000 platform at Beijing Genomics Institute (Shenzhen, China).

To identify DEGs between species in each sex and stage, we performed a differential expression analysis (see Supplementary Methods for details). To detect the predominant functional categories of the DEGs, we performed an enrichment analysis for each sex and stage using DAVID Bioinformatics Resources (52) with the D. melanogaster data as the background.

We checked whether the candidate genes in the species-monophyly regions detected by ARGweaver analysis were DEGs. We also listed genes involved in the basic development of insect genitalia and appendage-patterning genes tested in the development of various appendages in beetles. We determined whether these genes were differentially expressed (see table S6 for details).

Genomic sequence divergence among three species

We studied the genomic divergence among C. iwawakianus (four males), C. maiyasanus (six males), and C. uenoi (four males) using whole-genome sequence data from these species (table S2). The sequence reads of the three species were mapped on the draft genome of C. uenoi, SNPs were identified, and the average Weir-Cockerham Fst values between two species at 10-kbp intervals were obtained to identify regions of genomic divergence (see Supplementary Methods for details).

We also compared the 10-kbp average Fst values in scaffolds with and without genital QTLs for 11 scaffolds with genital QTLs [total, 22.6 million base pairs (Mbp)] and 8 scaffolds without (total, 37.8 Mbp) to determine the differences in genomic divergence patterns among the three species pairs. We studied the relationship between the 10-kbp average Fst value and the Ne of each 10-kbp interval, which was estimated in the ARGweaver analysis as described in the next section. In general, Ne is correlated with the nucleotide diversity (π) of the region, and a smaller Ne or π is related to a reduced genetic diversity within species, due to linked selection, and a larger Fst between species. Thus, the relationship between Ne and Fst may reveal selective divergence of genomic regions between species.

Candidate regions of species-specific genital morphology

Using the whole-genome sequence data of C. iwawakianus (n = 5) and C. maiyasanus (n = 6), we determined diverged genomic regions between these species in scaffolds with genital QTLs using ARGweaver version 0.8.1 (22), which estimates an ancestral recombination graph from multiple genomic sequences. We call such a region of a genomic sequence showing monophyly of species in its local genealogy a “species-monophyly region.” ARGweaver infers local genealogies (local coalescent trees) along a genomic sequence using a hidden Markov model and Markov chain Monte Carlo (MCMC) sampling. The MCMC sampling of ARGweaver was run with the following initial parameters: mutation rate = 2.8 × 10−9 (53), recombination rate = 2.9 × 108 (54), and effective population size = 2.0 × 105 (genome-wide average π × 1/4 × mutation rate). The first half of the samples was discarded as burn-in; subsequently, samples were taken from the MCMC chain with 20-step intervals. Scaffolds longer than 5 Mbp were split into 5-Mbp subsequences with 100-kbp overlaps, and results were combined after MCMC runs. The local Ne and recombination rates along the genome were calculated from MCMC samples using scripts bundled in ARGweaver. This analysis was performed for 11 scaffolds with genital QTLs and for another 10 without QTL as controls. If the action of selection favors specific alleles for each species, then genomic regions around the alleles are nonrandomly sorted and the local coalescent trees are more likely to attain monophyly than the genomic background. We calculated the posterior probability of monophyly of three species at each genomic position by obtaining the frequency of monophyly in local coalescent trees for each genomic position in the MCMC samples. We classified the regions of species monophyly into the following four categories: exon, intron, flanking regions (2500 bp upstream and downstream of genes), and noncoding.

RAD sequence divergence associated with divergence of genital morphology

To assess whether the QTL scaffolds for C. iwawakianus and C. maiyasanus were also divergent in other species in the C. iwawakianus–C. insulicola group with various genital morphologies, we measured sequence divergence between two pairs of parapatric species with long and short genital parts, C. arrowianusC. komiyai and C. insulicolaC. esakii. We mapped the RAD sequence reads onto scaffolds 237, 25, 385, 406, and 41106 using bowtie2 and determined SNPs using SAMtools. The Weir-Cockerham weighted Fst values of RAD sequences between species were calculated using VCFtools. We obtained the average Fst values for RAD sequences involved in every 50-kbp segment of each scaffolds. We evaluated the divergence of scaffold sequences among species in the C. iwawakianus–C. insulicola group by constructing a neighbor-joining tree for each scaffold using the Fst values.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/6/eaav9939/DC1

Supplementary Methods

Fig. S1. Linkage map constructed on the basis of F2 genotypes of the C. iwawakianus × C. maiyasanus cross.

Table S1. RAD sequences used in the present study.

Table S2. Whole-genome and transcriptome sequence data for the three Carabus (Ohomopterus) species.

Table S3. Species-monophyly regions detected by ARGweaver analysis.

Table S4. Genes associated with species-monophyly regions detected by ARGweaver analysis and their expression patterns in larvae and pupae of C. maiyasanus and C. iwawakianus.

Table S5. DAVID enrichment analysis for DEGs by sex and stage.

Table S6. Comparison of expression patterns of developmental genes between C. maiyasanus and C. iwawakianus larvae and pupae.

Table S7. Multiple regression analyses for the effects of presence or absence of genital QTLs and of species pair and local effective population size (10-kbp window; log10 scale) of 10-kbp average Fst among three species pairs.

Table S8. Genital dimensions of F2 individuals from crossing between C. iwawakianus and C. maiyasanus.

References (5575)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank H. Sugawara for assistance in experiments and J. Konuma, T. Niimi, and S. Nomura for advice and discussion. Funding: This study was supported by JSPS KAKENHI (grant nos. 23370011, 23128507, 25128707, 26251044, and 18H04010). Author contributions: T.S. designed the study, conducted the experiments, and analyzed the data. T.F. designed and conducted the data analyses. M.S., Y.T., and N.N. conducted the sampling and experiments. T.S. and T.F. wrote the manuscript, and Y.T. revised the manuscript. All authors read and approved the final manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The raw sequence reads used in the present study are available from the DDBJ Read Archive (DRA) of the DNA Data Bank of Japan (DDBJ) (BioProject; PRJDB2148, PRJDB5474, and PRJDB5403). Other relevant data are available via figshare (doi: 10.6084/m9.figshare.c.4349672.v2). Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article