The Effects of the Sex Chromosomes on the Inheritance of Species Speci c Traits of the Shape of the Copulatory Organ in Drosophila virilis and Drosophila lummei

Alex M Kulikov Koltzov Institute of Developmental Biology of Russian Academy of Sciences Svetlana Yu Sorokina Koltzov Institute of Developmental Biology of Russian Academy of Sciences Anton I Melnikov Koltzov Institute of Developmental Biology of Russian Academy of Sciences Nick G Gornostaev Koltzov Institute of Developmental Biology of Russian Academy of Sciences Dmitriy G Seleznev Papanin Institute for Biology of Inland Waters of Russian Academy of Sciences Oleg E. Lazebny (  oelazebny@gmail.com ) Koltzov Institute of Developmental Biology of Russian Academy of Sciences https://orcid.org/0000-0003-3136-0415


Background
Reproductive isolation contributes to the evolutionary process by allowing diverging species to accumulate genetic variation at adaptively important traits independently. Isolation is ensured by pre-and postzygotic isolating mechanisms, which differ in evolutionary origin and physiological basis. Postzygotic isolation arises as independent genetic variations accumulate in isolated populations and cause sterility and/or non-viability of hybrid progenies because of non-additive interactions of the alleles at various loci, while each allele does not exert a deleterious effect in the gene pool of its original population [1][2][3][4][5][6]. Prezygotic isolation results from the variation driven by sexual selection and prevents mating of individuals from different populations possessing different adaptations. In allopatric speciation, selection of traits acting at the copulation stage starts when postzygotic incompatibility has already formed. Reinforcement [7,8] and the sexual selection or sexual con ict [9][10][11][12] mechanisms can be involved in this case. In sympatric speciation, a prezygotic barrier should start forming before postzygotic isolating mechanisms are involved. It is noteworthy that experimental ndings agree with the assumption that prezygotic isolation accumulates at a higher rate as compared with postzygotic isolation. Coyne and Orr [13,14] estimated the parameters in experiments with many closely related Drosophila species and con rmed a higher accumulation rate for prezygotic isolation.
To discuss the formation rate and genetic basis of the two types of isolating barriers, it is necessary to understand the similarities and differences of their underlying processes. Postzygotic isolation is neutral or random for proper isolation of related species. An independent accumulation of variation in different evolutionary lineages leads to xation of the alleles that show incompatibility and, in the case of sex chromosome linkage, gives origin to Haldane's rule [15,16], which postulates a reduction in tness components for the heterogametic sex [17]. Low dominance values at incompatibility loci are required for Haldane's rule to occur concerning viability, while more complex interactions of dominance and sexspeci c incompatibility are required for Haldane's rule to occur for fertility [16]. The role that the sex chromosomes play as a main genetic structure element responsible for hybrid incompatibility is additionally illustrated by two other empiric regularities, which are known as the substantial X effect and sex-limited asymmetry of hybrid viability and fertility in reciprocal crosses [13,18].
Prezygotic isolation results from a selection at sex-limited traits, mostly those associated with mating behavior. Compared with autosomes, the X and W chromosomes are enriched in sexually antagonistic genetic variation because alleles exerting alternative effects in different sexes are more likely to be xed [19]. X-and W-linked variation is consequently more likely to be targeted by selection at sex-limited traits.
To explain the role that the sex chromosomes play in the evolution of isolating mechanisms, it is essential to consider the effects of all three evolutionary factors: mutagenesis, selection, and genetic drift. Chromosome variation due to mutation depends on the total evolutionary time that a given chromosome spend in the male genome and shows different extreme magnitudes in the case of the sex chromosomes, being the highest in the Y and Z chromosomes and the lowest in the X and W chromosomes [20]. Thus, the action of this evolutionary factor does not explain the substantial X effect but implies that the Y and Z chromosomes accumulate substantial genetic burden and degrade in the absence of recombination [21].
The conclusion is based on the assumption that all chromosomes accumulate variation due to mutation at the same rate in males. Studies of meiotic sex chromosome inactivation (MSCI) in males still make it possible to assume that heterochromatin forming in a major part of the sex chromosomes in pachytene [22] prevents repair mechanisms from fully restoring heterochromatin sequences, thus substantially increasing the mutation rate of the two sex chromosomes relative to that of the autosomes. If so, the circumstance provides a plausible explanation for the substantial X effect, faster-X theory [23], and faster male theory [24].
In contrast, the signi cance of genetic drift is well explained by the chromosome population effective size, which is smaller in the sex chromosomes compared with autosomes [25]. Alleles associated with pre-and postzygotic isolation are accordingly more likely to be xed on the sex chromosomes under the in uence of random genetic drift processes.
Selection differently affects sex chromosome variation during the formation of prezygotic vs. postzygotic isolation. The xation of new alleles, which is essential for incompatibility to arise, implies the effect of positive selection only. Postzygotic incompatibilities are not selected directly, but their choice is mediated by functional associations of particular genes with viability and fertility traits. The sex chromosomes are enriched in sexually antagonistic genes [26,27] and may more often be targeted by selection at differentially expressed traits. Genes involved in gametogenesis act as targets in the case of postzygotic isolation; and genes involved in mating behavior, in the case of prezygotic isolation. In species with XY sex determination, the X chromosome undergoes demasculinization [27,28], while postzygotic isolation mechanisms act more often in males, in agreement with Haldane's rule. The dominance theory [6, 15,29] postulates that recessive variations of X-chromosomal incompatibility loci are expressed in hemizygous males and explains well a viability decrease unrestricted to sex-linked traits. The discrepancy between X-chromosome demasculinization and fertility impairment in males are possible to explain by nonadditive interactions of alleles and, primarily, recessive negative epistasis. Epistatic interactions normally ensure homogametic sex-limited expression and are altered in a hybrid genome, leading to aberrant expression in the heterogametic sex to cause infertility and a drop in tness [27].
Following the drive hypothesis, selection may additionally act indirectly, by suppressing sex chromosome meiotic drive, which results from sex chromosome evolution and permanent internal con ict between sex determination systems [30][31][32][33]. For example, various systems of sex chromosome segregation distorters have been identi ed in the X chromosome and autosomes in Drosophila species by an incompatibility locus analysis; the systems include pseudogenes, repeat regions [31,32,34], and heterochromatin blocks responsible for failure of sister X chromosomes to segregate during mitosis in the hybrid genome [35,36].
The role that the Y chromosome plays in the selection-driven formation of prezygotic isolation depends to a substantial extent on the sex wherein a particular trait is expressed. Given that the X chromosome experiences demasculinization and a direct effect of positive selection, X-chromosome variation is mostly restricted to female choice traits [38]. However, demasculinization of the X chromosome is not absolute, and the hemizygous status of the X chromosome dramatically increases the effect of selection in males [20,39], while not excluding an accumulation of alleles associated with particular characteristics of male mating behavior. Asymmetry of male success in reciprocal interspeci c crosses of related species [40][41][42] formally resembles Haldane's rule but differs in arising directly at the stage of choosing a partner from a related species, being based on the differences accumulated in the respective genomes. Unequal partner preferences in reciprocal crosses depend on the reaction norm of each of the species concerning partner choice traits [43] and do not suggest linkage with the sex chromosomes.
In Drosophila, the X chromosome has been shown to play a role in several species-speci c features that are observed in male mating behavior and prevent heterospeci c mating, the set including various elements of the courting ritual and various modalities. When cuticular hydrocarbons (CHCs) are utilized as cues in D. virilis and D. americana texana, a leading role in choosing the mating partner is played by the male, which discontinues tapping and loses interest in the female if the female is heterospeci c. The capabilities of recognizing the female type and discontinuing the courting ritual are associated mostly with the X chromosome [44]. Sharma and colleagues [45] observed signi cant sex dimorphism in CHC pro le distribution in isofemale D. simulans lines, thus con rming that the sex chromosomes are involved in regulating expression of the genes belonging to the metabolic cascades of hydrocarbon synthesis. Characteristics of the male mating song were studied in crosses between closely related Drosophila species of the sophophora, obscura, and virilis groups [46][47][48][49][50][51][52][53]. The results have shown that all chromosomes, including the sex chromosomes, contribute to the species-speci c variation in mating song elements, while the set of chromosomes or loci associated with the differences vary among related species pairs. A role of the X chromosome was demonstrated for the majority of the variants tested. Arbuthnott [54] analyzed a large set of reproductively isolating traits with known genetic architectures and concluded that the majority (25 of 36) of the traits are controlled by several loci of relatively large effects. The conclusion was that the majority (70%) of traits responsible for intraspeci c and interspeci c differences are controlled by a few loci of relatively large effects. The variation due to interspeci c differences displayed no predominant association with the sex chromosomes and substantially differed from intraspeci c variation. Epistatic interactions were found to play a signi cant role in the intraspeci c and interspeci c variation of the male courtship song.
Copulation as a nal stage of mating behavior determines the e ciency of insemination of the female and the contribution of the male genotype to the gene pool of the next generation. The copulation e ciency depends on the copulation duration and the insemination reaction. A crucial role in the latter is played by male accessory gland secreted proteins (Acps), which evolve at a far greater rate as compared with hemolymph proteins [55][56][57][58][59][60]. A total of 46 genes are known to code for Acps in D. melanogaster, and only six of them are on the X chromosome. A decrease in female fertility in heterospeci c crosses of D. virilis and D. americana is associated with three D. americana and two D. virilis recessive autosomal genes [61]. The nding agrees with the idea that isolating barriers rapidly arise from a divergence of a few genes and demonstrates that the X chromosome is weakly involved in the formation of prezygotic isolating barriers.
The copulation duration directly depends on how well the male copulatory system matches the female genitalia in shape. Jagadeeshan and Singh [62] observed that in crosses between the closely related species D. melanogaster and D.
simulans, mating behavior was directly associated with the species-speci c shapes of the posterior lobes of the genital arch and the cercus, which are external structures of the Drosophila male copulatory system. The time course of copulation stages depended on the mechanic coupling of the male and female genitalia, and the coupling depended on the male external genital structures. The nal stage of tight genital coupling was 2-5 times shortened in heterospeci c crosses, while the preceding unstable coupling stage was prolonged. The ndings indicate that female insemination may be incomplete in heterospeci c crosses and that a particular mechanism exists to allow mating ies to be sensitive to the speci city of their contact. Laurie and colleagues analyzed the genetic basis of species divergence to the shapes of male copulatory system elements in Drosophila. A series of studies was carried out with D. simulans and D. mauritiana to genetically analyze the species-speci c distinctions in the shape of the epandrium, which is an external male genital structure, and at least eight linkage groups involved were localized to the X chromosome and chromosomes 2 and 3 [63-66]. Additive interactions were mostly observed for the respective loci.
The shape of the male genitalia is the most rapidly evolving among all morphological characters, and its control system is a target of sexual selection [67]. However, the shape of the epandrium, which is the posterior lobe of the genital arch, has actually been used as the only model to study the genetic basis of species-speci c differences in the shape of the copulatory system in D. simulans and D. mauritiana. Do the X chromosome and autosomes play the same role in speciesspeci c traits related to the shape of the copulatory system in other species groups? Does the Y chromosome contribute to the inheritance of these traits? In this study, we used interspeci c crosses between D. virilis and D. lummei and backcrosses of F1 hybrid males with D. virilis females to evaluate the contribution of the sex chromosomes and the total contribution of the autosomes to the degree of dominance of the D. virilis and D. lummei phenotypes in traits related to the shape of the male copulatory system.

Results
Correlated variability of traits of the shape of the copulatory organ To study the structure of correlations among the traits, and the latent factors responsible for a correlated variation, total variance over all samples was examined by a factor analysis. Seven factors were taken based on the scree plot and together accounted for 64.4% of the total variance ( Traits included in factor loadings with weights higher than 0.5 are shown; * indicates the traits that had weights higher than 0.7.
The eigenvalue and percent variance accounted for by a respective factor are shown in parentheses.
Based on the traits with the greatest factor loadings, the factors were classi ed as follows.
Factor 1 characterized the variation in the width of the apodeme outline; Factor 2 characterized the variation in the traits related to the apodeme bend; Factor 6 characterized the variation in the cook angle and the apical part of the aedeagus outline; and Factor 7 characterized the variation in traits related to the dorsal bend in the distal part of the aedeagus outline and the lowermost point of the ventral bend in the same aedeagus part. IMP 3,18,19,23,24,26,27, and 28 were highly speci c ( Table 1, Highly speci c).
Generalized Estimate of the Dependence of the Revealed Variability of the Phenotype of F 1 and F b Hybrids on the Genotype RDA was carried out to relate the set of the variables that determine the factor structure values to the variables that characterize the genotype. In the analysis, the ordination of the genotypes and factors was performed in the space de ned by correlations between factor weights, which were weighted sums of trait weights, and linear combinations of genotype scores, which were obtained for the genotypes characterizing the respective point in trait distribution. The center of gravity of the genotype distribution was brought into coincidence with the centers of gravity of the distribution of the factor structures. As is seen from Fig. 1, a two-dimensional space is de ned by the orthogonal factor pair ML3-ML4 (X) and ML6, ML5-ML2, ML7 (Y). The genotypes showed the following arrangement in the two-dimensional factor space: +X: AAAH, ABAH, +Y(− X): BBBLu, AAHH, −Y(− X): ABHLu, BAHVi (genotype abbreviations are as in Table 5). Genotypes ABHH and AAAVi occur in the region of medium values. The angle between vectors on the plot and between genotype positions is proportional to the correlation between them. Therefore, the variation in traits of backcross males homozygous for the autosomes (with the D. virilis X chromosome) positively correlates with Factor 3 and negatively correlates with  Table 5. The chromosomes and paternal genotype are indicated in the following order: X chromosome, Y chromosome, autosomes, male parent identity.
The results are supported by the distribution of the genotypes on a scatter plot (Suppl. Figure 1 Two approaches were used to more precisely evaluate the effect that the sex chromosomes exert on trait expression in the shape of the male copulatory system. First, ANOVA was performed to compare the dominance of parental species-speci c phenotypes in progenies from reciprocal crosses between D. virilis and D. lummei and backcrosses of F 1 males with D. virilis females. All genotypes had the same autosome set and differed in sex chromosome combination and the identity of the male parent (the original parental species or the F 1 hybrid). Second, ANOVA and post-hoc tests were performed using the genotype at the sex chromosomes, the genotype at the autosomes, and combinations of these factors, including the male parent identity, as independent variables.
The phenotypes of males obtained in direct and reciprocal crosses and backcross males heterozygous for all autosomes were compared with the phenotypes of males of the parental species. The results are summarized in Table 2. The logic of estimating the degree of dominance for a trait has been described previously [68]. Based on the distribution of the hybrid and parental genotypes over groups isolated by post-hoc comparisons, it is possible to identify the following variants: incomplete dominance, the dominance of the D. virilis or D. lummei phenotype, and lack of difference among all phenotypes. The resulting data clustering variants were ranged by the degree of phenotype dominance. All cases where a genotype in question appeared together with the parental genotypes in one group were considered to suggest no signi cant difference (ns) in evaluating the degree of dominance. The variants f1  Table 2 Dominance of the copulative system shape-related traits as dependent on the sex chromosome composition in D. virilis/D. lummei hybrid males, heterozygous for the autosomes.

Factor
Sign Traits are grouped according to their maximal weights in the factors extracted isolated by the maximum likelihood method (Table 5)

Factor
Sign Traits are grouped according to their maximal weights in the factors extracted isolated by the maximum likelihood method ( expression. The effect of the male parent identity was evaluated by comparing F 1 males with genotype X Vi Y Lu and backcross males with the same genotype, given that crossing over is absent in males and the males in question were entirely identical in chromosome composition. Differences in parental phenotype dominance were observed for 16 out of the 30 traits, and the dominance character changed to the opposite one in the case of apodeme declination angle. A substitution of the D. lummei Y chromosome for its D. virilis counterpart in backcross males similarly changed the dominance character in 16 traits. Of these, ten traits, which mostly loaded on Factors 1, 3, 4, and 5, showed a change to dominance of the opposite species-speci c phenotype relative to the species origin of the Y chromosome. Phenotypic comparisons of F 1 progenies showed that simultaneous substitution of the sex chromosomes changed the character of dominance at 12 traits, of which three (IMP 11, IMP20, and β) changed their dominance status to the opposite one, according to the species origin of the X chromosome. Other changes were less distinct and included rstly, intermediate dominance was shifted toward dominance of one of the parental species. Secondly, the phenotype was changed relative to that of the parental species so that a homogeneous variance group within one of the parental species and offspring with a particular genotype was replaced with a group that included offspring with an alternative genotype and both of the parental species ( To study the effect on particular phenotypic traits for each of the factors, MANOVA was performed with forced incorporation of all four predictors and the effect of the ChrY*Aut interaction. The results are summarized in Suppl. Table 4. The majority of the traits showed an association with the identities of the X chromosome and autosomes at a signi cance level p < 0.05; half of the traits were affected by the Y chromosome and male parent identities; and one-third of the traits, by the interaction of the Y chromosome and autosomes. A cooperative effect of these variables was mostly observed for the groups of traits extracted by the factor analysis. An analysis of the factor values con rmed the effect in the cases where the majority of the traits determining the respective factor structure signi cantly depended on the given variable. Note that a dependence on the hereditary factors was not con rmed for the apodeme width-related traits (F1).
The apodeme declination and curvature (F2) depend on the identities of the Y chromosome and male parent. Distinct species-speci c traits of the aedeagus and parameres (F3) depend on the X chromosome, autosomes, and male parent identity. The shape of the ventral bend in the proximal part of the aedeagus outline (F4) is determined by the interaction between the Y chromosome and autosomes. Traits related to the height of the dorsal bend in the aedeagus outline (F5) and the shape of the cook and the outline bend over the cook (F6) showed a similar dependence on the Y chromosome, autosomes, and male parent identity; an additional dependence on the X chromosome is speci c to F5-related traits. Traits describing the shape of the distal part of the outline (F7) depend on the autosomes, male parent identity, and their interaction.
The effect of the male parent identity is possibly mediated through the epigenetic marking of chromosomes in interspeci c hybrids during gametogenesis and a subsequent effect of the resulting signatures on the ontogeny of offspring. As mentioned above, an independent effect of the Y chromosome is expected to be insigni cant based on the  Table 5). The results obtained for each particular trait are provided in Supplementary; data on latent and highly speci c traits are summarized in Table 3. Table 3 Paired permutation test with the Bonferroni correction to evaluate the optimal genotype partitioning into groups homogeneous in latent variables and highly speci c variables (up to 100 000 iterations in each case Designations of the factors -as in Table 5 The model implies that each hereditary factor exerts a discrete effect on the traits in question, depending on the genotype. The expected effect of each factor was analyzed using the indicator variables listed in Table 5. Designations of the factors -as in Table 5 The model implies that each hereditary factor exerts a discrete effect on the traits in question, depending on the genotype. The expected effect of each factor was analyzed using the indicator variables listed in Table 5. Factor X→Y(epist) has two categories for the A 1 indicator variable, the effect being present or absent. A signi cant difference in these categories between samples con rms the signi cant effect of the interaction between the sex chromosomes. The effect of the conspeci c D. virilis sex chromosomes is examined versus the effects of all other combinations in this case.
Factor Y + Y→AUT(dom.epist) has two similar categories for the A 5 indicator variable. The independent effect of the Y chromosome alone on the traits is close to zero. Therefore, a signi cant difference in the categories between samples con rms the signi cant effect of the interaction between the Y chromosome and autosomes. The effect of the D. virilis Y chromosome and autosomes is examined versus the effects of all other combinations in this case.
With all subsequent factors, a signi cant difference between the extreme variants con rms that the respective factor signi cantly affects the trait in question. A signi cant difference between intermediate and other genotype groups in the absence of differences between the extreme groups indicates that the effect in question depends nonlinearly on the genotype.  (factor P→AUT(add) + P→Y), A6 (factor P→X + P→AUT(dom)), and A7 (factor X→AUT(dom.epist) + X + AUT(dom)). It is clear that the phenotypic features of this genotype are determined by combined effects of the three factors, and its clustering with other genotypes may therefore be distorted in the case of factor P→AUT(add) + P→Y.
Thus, the estimates con rm that epistatic interactions of the sex chromosomes and autosomes and epigenetic effects of the male parent origin from interspeci c crosses in uence the expression of D. virilis species-speci c traits in the shape of the male copulatory system.

Discussion
How do traits determining the shape of the male copulatory system become a target of selection? As mentioned in the Background section, the e ciency of female insemination in insects depends on how well the female and male genitalia match each other [62]. Incomplete insemination must make repeated mating more likely, lead to displacement of sperm from the previous mating, and facilitate e cient selection against the incomplete insemination-associated genotype in Drosophila. Hoikkala and colleagues [70] have analyzed the within-population variation of mating duration in D. montana and associated the duration of the rst mating with female resistance to repeated mating. We have shown that sensory microchaetae on the ventral surface of the aedeagus mediate the association between evolutionarily signi cant parts of the Drosophila copulatory system with mating behavior [68].
Both mating behavior and the shape of the copulatory system are not directly involved in adaptation but are maintained close to the adaptive optimum of a population as a result of apparent stabilizing selection [71,72]. In other words, the selection at these characters acts through adaptively valuable traits characteristic of a particular group of individuals. This selection takes the form of sexual selection because the characters are expressed speci cally as predictors of sexual reproduction e ciency. When a new adaptation forms and the adaptive norm changes rapidly, selection changes to directional or disruptive, but still remains indirect. It is, therefore, possible to expect that the variation observed experimentally will correspond to the variation in quantitative traits affected by one of the above selection types.
The variation maintained by stabilizing selection is equally well described by two models, a joint-effect model and an in nitesimal model. In the former, the variation is associated with effects of moderate-frequency alleles, which are nearly neutral in terms of tness and substantially in uence the trait in question. In the latter, the variation involves many genes that each exert a minor effect on tness and act additively. Modeling of ample experimental data has shown that well reproducible results are obtained with both of the models, none of them is possible to prefer over the other [72]. An analysis of QTLs for commercially valuable traits in farm animals supports well the conclusion that recessive variation maintained in a population plays an important role [73][74][75]. The majority of the traits are related to morphological and physiological characteristics affected by stabilizing selection. The facts that a variety of farm animal breeds formed rapidly on an evolutionary scale, that homologous haplotypes are responsible for similar traits in different breeds and are present in populations of founder species, and that epistatic interactions are observed between QTLs indicate that diversity at quantitative traits is maintained because neutral variation is preserved in populations.
Lower estimates could be obtained for genetic variation at traits under stabilizing selection when the effects of suppressing epistasis are underestimated, the fact being masked by overestimating the effect of stabilizing selection and/or underestimating the magnitude of variation due to mutation [76,77]. Moreover, opposite epistatic effects may occur between tightly linked genes in a QTL [78] to reduce their total effect observed, and opposite effects occurring between different QTLs substantially increase the total variance.
The formation of evolutionary new and usually adaptively signi cant traits has another genetic dynamics. Fisher [79] has noted that adaptation is not selection, meaning tness-based selection that eliminates the least t individuals. Adaptation is characterized by the effect of positive selection, which facilitates the progress of a population to an optimal phenotype. Additive, dominant, epistatic, and epigenetic components of variation underlie the phenotypic differences in our experiment. The experiment was not designed to exactly determine the variance components. However, given that the homozygous or heterozygous autosome sets were identical in the males examined, it is clear that changes in the dominance of the parental phenotypes are associated with epigenetic effects of the male parent identity and epistatic effects of the X and Y chromosomes. Neither primary nor latent traits remained phenotypically the same in genotypes with different combinations of the sex chromosomes and different male parent identities, while the autosome set was identical.
The effects of sex chromosomes and autosomes were speci c to different groups of traits. For example, the additive effect of recessive autosomal alleles determined the traits of the aedeagus (apart from those characterizing the ventral bend in the proximal part of its outline) and parameres. A combined effect of dominant autosomal alleles, the X chromosome, and its dominant epistatic interaction with the autosomes was observed for traits of the aedeagus and parameres and the bend and declination of the apodeme. It is of interest that the two factors exerted similar effects in genotypes with intermediate values of indicator variables. This was dominance of the D. lummei phenotype in latent trait F3, which combined the most distinct species-speci c traits; certain primary traits that were incorporated in latent traits F4 and F5 and similarly showed species speci city [90]; and paramere width (IMP 27). Dominance of the D. virilis phenotype was characteristic of latent trait F7; several primary traits incorporated in latent traits F5 and F6, and paramere width in the ventral bend region (IMP 24). The examples show that the dominant component of the autosomes contributed to the effect of additive interactions and that its contribution depended on the epistatic interaction with the X chromosome.
Following Huang and Mackay [91], the components of variance are impossible to strictly extracted isolate in the majority of cases, especially when the experimental design is not an orthogonal one, which at least formally ensures uncorrelated variance components. Our estimates are therefore total effects of several components of variance. Minimization of the autosomal effect in the case of dominance of the D. lummei phenotype and, oppositely, minimization of the effect of the X chromosome and its epistatic in uence on dominant autosomal alleles show that the X chromosome plays alternative roles in the variation of different traits. Similar dominance changes observed for genotypes heterozygous for the autosomes when estimating the role of the additive component of variance con rms that the dominant component of the autosomes contributes to the species-speci c variation.
The effect that epistatic interactions between the X and Y chromosomes exert on species-speci c traits was con rmed for half of the primary traits and latent traits F3 and F5, which incorporate taxonomically signi cant traits. It should be noted here that unequivocal interpretation is impossible for the results indicating that the D. virilis phenotype is enhanced in males with conspeci c sex chromosomes. Given that different variance components may contribute to this phenomenon, an independent effect can be expected for genes of the D. virilis X chromosome. For example, Carson and Lande [92] have shown that the formation of an evolutionarily new secondary sex characteristic (an additional row of bristles on the male tibia) is determined to the extent of 30% by sex-linked genes in a natural Drosophila silvestris race. As for the traits that contribute to isolating barriers, the most detailed data are available for sterility-related traits. The examples below illustrate interspeci c hybrid male sterility as a model of interactions between the autosomes and sex chromosomes. A cluster of Stellate sequences on the D. melanogaster X chromosome codes for a homolog of the β subunit of protein kinase CK2, and its overexpression causes male sterility. RNA interference mediated by Y-linked Su(Ste) repeats prevents Stellate overexpression [93,94], demonstrating that male fertility strongly depends on the interaction between the X and Y chromosomes. A similar model of fertility regulation utilizes the Odysseus-site homeobox protein (OdsH), which is encoded by an X-chromosomal locus in species of the melanogaster group. OdsH binds to heterochromatin sequences of the Y chromosome. In many cases, the presence of heterospeci c X and Y chromosomes leads to decondensation of Ychromosomal heterochromatin, dramatically distorts the expression of autosomal genes, and causes sterility [34,95]. The genetic system also illustrates the interaction between the X and Y chromosomes.
The effect of epistatic interactions between the sex chromosomes and recessive alleles of the autosomes was con rmed for the majority of latent and primary traits. A leading role of epistatic interactions of the X chromosome is evident from the observation that the D. virilis phenotype dominated at the majority of traits in males with genotype X Vi /Y Lu Aut Vi /Aut Vi .
A substantial role of epistatic interactions between the X chromosome and autosomes has similarly been demonstrated for another trait involved in prezygotic barriers, namely, a species-speci c male courtship song pattern in Drosophila species of the montana phylad [50]. There is evidence that the sex chromosomes exert a signi cant regulatory effect on expression activity of autosomal genes. Expression of X-chromosomal genes in hemizygous D. simulans males depends on trans regulations to a substantial extent [96], and trans effects and joint action of cis and trans effects of the X chromosome on genome-wide expression activity are the second most important to trans-regulatory effects of chromosome III in D. melanogaster males [97].
The important role that the Y chromosome plays in regulating the phallus shape-related traits in interspeci c hybrids is possible to directly associate with trans-regulatory activity. The Y chromosome harbors only 23 single-copy protein-coding genes, and 13 of them are strongly restricted to the testis in expression and are mostly associated with hybrid male sterility [98][99][100][101]. Additive variation of the ten other genes should make a vanishingly small contribution to the phenotype at quantitative traits of morphological structures. However, epistatic interactions of Y-chromosomal sequences with the X chromosome and autosomes have received experimental support. For example, experiments were performed to evaluate the activity of the male-speci c lethal proteins/roX1,2 RNA complex, which is responsible for dosage compensation of Xlinked genes in Drosophila males. The Y chromosome proved to affect the viability in roX1, roX2 mutant males, the effect depending on the Y-chromosome source [102]. The viability was low in males the paternal Y chromosome and high in males with the maternal Y chromosome. The result indicates that the Y chromosome modi es dosage compensation through roX1, roX2-mediated modi cation of heterochromatin [103] and/or recognition of the X chrome by the entire malespeci c lethal proteins/roX1,2 RNA complex. The results support our nding of a substantial paternal effect, which acts independently in the majority of traits or in combination with the X-chromosome effect in some other cases. Finally, Lemos and colleagues [101,104] have directly demonstrated the regulatory activity of the Y chromosome in a study of differential genome expression activity under the Y-chromosome in uence.
Traits of the apodeme, which is a muscle attachment site and internal part of the copulatory system, are the least associated together by evolutionary variation and chromosome effects. It is noteworthy that chromosome effects are mediated by the epigenetic effect of male parent identity in the case of the apodeme traits. Note that signi cant epigenetic effects were similarly observed for all other groups of correlated traits. The epigenetic effect may be related to the mechanism of meiotic sex chromosome inactivation (MSCI) or, in a broader sense, meiotic silencing of unpaired chromatin/DNA (MSUC/MSUD) [22,31], which inactivates the chromosome regions with altered meiotic synapsis in leptotene. The MSCI mechanism and its substantial evolutionary role in organizing the chromosomes and facilitating postzygotic isolation has been the matter of extensive discussion [22,31,33,105]. MSCI is typically compensated for by the higher transcriptional activity of X-linked genes in Drosophila because potent promoters are overrepresented in the X chromosome [106,107]. The formation of heterochromatin blocks is often associated with sterility in hybrids from interspeci c crosses and death of hybrid progenies; defects in mitotic chromosome segregation are observed hybrids during their embryo development [22,34,35]. Interspeci c differences cannot accumulate at once and may be seen as changes in gene expression activity and alter the trait magnitudes in early species divergence. Epigenetic signatures are preserved in chromosomes at postmeiotic stages and may be inherited through generations [108,101]. Although the regulation of genome-wide epigenetic states is often associated with the Y chromosome [101,109,110], a formal role of the interspeci c hybrid status is noteworthy in our case, leading to distorted synapsis of divergent chromosomes and causing the formation and further preservation of noncanonical heterochromatin regions with altered expression activity.

Conclusions
Thus, most of the traits demonstrated dominance of the D. virilis type. It can be comprehended, given the region of D.
virilis and its sibling species origin -Southeast Asia. Among the sibling species, D. virilis has the shape of the copulatory organ the closest to ancestral one. Concerning the chromosome effects on the expression of the shape traits, all types of genetic interactions between the chromosomes were revealed in this study, namely, additive, dominance, epistatic and even epigenetic. The effects of sex chromosomes and autosomes were speci c to different groups of traits. The X chromosome plays alternative roles in the variation of different traits. Epistatic interactions between the X and Ychromosomes exert on species-speci c traits. The two sex chromosomes act unidirectionally to shift phenotype dominance in the conspeci c direction. The effect of the Y chromosome is mediated to a substantial effect by its epistatic interactions with the X chromosome and autosomes. Epistatic interactions of the sex chromosomes and autosomes and epigenetic effects of the male parent origin from interspeci c crosses in uence the expression of species-speci c traits in the shape of the male copulatory system. It can be assumed that sexual selection for speci c genes associated with male traits implemented in the courtship ritual prevents the well-known effect of demasculinization of the X chromosome. Direct crosses

Methods
A Vi , A Lu , AH are the D. virilis, D. lummei, and heterozygous autosomes, respectively.
Backcross males heterozygous at all autosomes and fully homozygous backcross males were selected for further analysis. Backcrosses with heterozygous F 1 males were performed to exclude autosomal recombination events and to allow a rigorous analysis of the effects for variants of homozygous and heterozygous combinations of autosomes of the parental species.
All crosses were carried out at 25 °C; a standard food medium and glass tubes of 22 mm in diameter were used; the progeny density was 50-70 ies per tube.

Analysis of morphological structures
The mating organ was dissected using thin steel needles in a drop of water under a binocular microscope at a magni cation of 12 × 8. To remove adipose structures preparations were incubated in boiling 10% NaOH. Morphometry was carried out using organ images, which were obtained using a Jen-100C electron microscope in the scanning mode at an accelerating voltage of 40 kV and an instrumental magni cation of 300-500×. The sagittal view of the phallus was conventionally divided into four areas: the aedeagus body, gonites, apodeme, and cook. A coordinate grid was superimposed on the image. A conventional point at the junction of the aedeagus area, gonites, and apodema (referred to as the central point) was used as a landmark. Morphometric parameters (MPs) were obtained as distances between the intercrosses of coordinate lines with each other and the image outline. A scheme of measurements is shown in Fig. 2; MP characteristics are summarized in Table 4. The declinations of the cook and apodeme (axes MP 2 and MP 28) from axis MP 1 were expressed in radians and designated α and β, respectively. To exclude the dimensionality factor, MP indices (MPIs) were calculated according to a method used previously to evaluate the inter-and intraspeci c variations of the genitals in the virilis species group (1). MPIs were obtained as MP-to-MP 1 ratios and were numbered according to the numbers of respective MPs. The traits expressed in radians were included in the analyses without normalization.

Statistical analyses
Statistical analyses were carried out using the program IBM SPSS Statistics v. 23 and the R 3.3 statistical analysis environment with the packages "vegan", "lmPerm", "psych", and "rcompanion". where X ij is the trait value in the sample of a particular genotype; µ is the mean; ε is the fraction of variance unexplained; A n values are the weights or indicator variables corresponding to the regression coe cients E n ; and n indicates the respective genetic factor or factor combination: X and Y, independent effects of the respective sex chromosomes; AUT(add), the additive effect of recessive autosomal alleles; and AUT(dom), the dominant effect of the autosomes. Pairwise interactions of the factors were also considered, including X→Y(epist), the epistatic effect of the interaction of the X and Y chromosomes; AUT(add) (or AUT(dom))_P, the epigenetic effect of the paternal genotype on the additive or dominant effect of the autosomes; Y (or X)→AUT(rec.epist), the epistatic effect of the Y (or X) chromosome on the recessive autosomal alleles; Y (or X)→AUT(dom.epist), the epistatic effect of the Y (or X) chromosome on the dominant autosomal alleles; and P→X, the epigenetic effect of the paternal genotype on the effect of the X chromosome. A plus sign indicates a combination of the respective effects. The effect of the paternal genotype alone on the traits under study was not considered because epigenetic effects were expected only for a combination of this factor with the offspring genotype. The genetic factors were treated as categorical variables; their signi cance was assessed by permutational ANOVA (PermANOVA). A paired permutation test was employed in post-hoc comparisons, and the Bonferroni correction was used for multiple comparisons. The maximal number of iterations was 100 000; the signi cance level was 0.05.
Several effects were combined in one variable because the study was not designed to examine all possible combinations of the factors and because several vectors determining the A n indicator variables were collinear. To combine several effects in one joint effect, the values of the corresponding vectors were summed. Backcrosses to restore the D. virilis genotype serve to evaluate the effect of the D. virilis chromosomes on the resulting phenotype. The indicator variables were therefore taken to be 1 for the D. virilis genotype and 0 for the D. lummei genotype when considering effects of the X and Y chromosomes and a dominant effect of the autosomes and to be 1 for homozygotes for the D. virilis chromosomes, 0 for heterozygotes, and − 1 for homozygotes for the D. lummei chromosomes when considering additive effects of the autosomes. To allow for possible spermatogenesis defects and subsequent epigenetic marking of the chromosomes in interspeci c hybrids, the indicator variables for the factor Paternal Genotype were taken to be 1 for homozygous males of the parental species and 0 for heterozygous F1 males. The resulting matrix of genotype-dependent vector values is shown in Table 5.  A factor analysis was used to capture a covariance of traits; factors were extracted by the maximum likelihood method; factor loadings were examined by rotating principal component axes via the Promax procedure with the Kaiser normalization; estimates were obtained by regression analysis. The number of factors to be extracted isolated was determined using a scree plot of the eigenvalues of the trait correlation matrix.
Because a distinct elbow was not observed in the plot, we chose the factors whose eigenvalues were above a simulation curve obtained by a parallel analysis with 1000 permutations. The results of the factor analysis and their associations with genotypes were visualized by redundancy analysis (RDA), using the aggregate matrix of mean estimates for each factor. A principal component analysis (PCA) was performed for the initial MPs. The rst two components were used in visualization. Availability of data and materials

Abbreviations
The original datasets are available in the supplemental le "Kulikov_Raw Data.xlsx."

Competing interests
The authors declare that they have no competing interests.

Funding
We gratefully acknowledge funding provided by the Program "Biodiversity of Natural Systems" (Subprogram "Gene Pools of Living Nature and Their Conservation") of The Presidium of the Russian Academy of Sciences to AMK. The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript. Authors