The effects of the sex chromosomes on the inheritance of species-specific traits of the copulatory organ shape in Drosophila virilis and Drosophila lummei

The shape of the male genitalia in many taxa is the most rapidly evolving morphological structure, often driving reproductive isolation, and is therefore widely used in systematics as a key character to distinguish between sibling species. However, only a few studies have used the genital arch of the male copulatory organ as a model to study the genetic basis of species-specific differences in the Drosophila copulatory system. Moreover, almost nothing is known about the effects of the sex chromosomes on the shape of the male mating organ. In our study, we used a set of crosses between D. virilis and D. lummei and applied the methods of quantitative genetics to assess the variability of the shape of the male copulatory organ and the effects of the sex chromosomes and autosomes on its variance. Our results showed that the male genital shape depends on the species composition of the sex chromosomes and autosomes. Epistatic interactions of the sex chromosomes with autosomes and the species origin of the Y-chromosome in a male in interspecific crosses also influenced the expression of species-specific traits in the shape of the male copulatory system. Overall, the effects of sex chromosomes were comparable to the effects of autosomes despite the great differences in gene numbers between them. It may be reasonably considered that sexual selection for specific genes associated with the shape of the male mating organ prevents the demasculinization of the X chromosome.


Introduction
Reproductive isolation contributes to the evolutionary process by allowing diverging species to accumulate genetic variation independently, including for adaptively important traits. 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 nonadditive interactions of the alleles at various loci. In contrast, each allele does not exert a deleterious effect on the gene pool of its original population [1][2][3][4][5][6]. Prezygotic isolation results mainly from the variation driven by sexual selection and prevents the mating of individuals from different populations possessing different adaptations. In allopatric speciation, the selection of traits acting at the copulation stage starts when postzygotic incompatibility has already formed. Reinforcement [7,8] and the sexual selection or sexual conflict [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 findings 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 confirmed a higher accumulation rate for prezygotic isolation. Selection affects sex chromosome variation differently during the formation of prezygotic versus postzygotic isolation. The fixation of new alleles, which is essential for incompatibility to arise, implies the effect of positive selection only. Postzygotic incompatibilities are not selected directly, but, functional associations of particular genes with viability and fertility traits mediate their choice. The sex chromosomes are enriched in sexually antagonistic genes [15,16] 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 act as targets in the case of prezygotic isolation. In species with XY sex determination, the X chromosome, as a rule, undergoes demasculinization [16,17], while postzygotic isolation mechanisms act more often in males, in agreement with Haldane's rule. The dominance theory [6,18,19] 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 can be explained by nonadditive interactions of alleles and, primarily, recessive negative epistasis. Epistatic interactions usually 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 fitness [16].
Following the drive hypothesis, selection may additionally act indirectly, by suppressing sex chromosome meiotic drive, which results from sex chromosome evolution and permanent internal conflict between sex determination systems [20][21][22][23]. For example, various systems of sex chromosome segregation distorters have been identified in the X chromosome and autosomes in Drosophila species by an incompatibility locus analysis; the systems include pseudogenes, repeat regions [21,22,24], and heterochromatin blocks responsible for the failure of sister X chromosomes to segregate during mitosis in the hybrid genome [25,26].
The role that the X 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, Xchromosome variation is mostly restricted to female choice traits [27]. 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 [28,29], while not excluding an accumulation of alleles associated with particular characteristics of male mating behavior. Asymmetry of male success in reciprocal interspecific crosses of related species [30][31][32] formally resembles Haldane's rule. However, it 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 [33] and do not suggest linkage with the sex chromosomes.
In Drosophila, the X chromosome has been shown to play a role in several species-specific features that are observed in male mating behavior and prevent heterospecific mating, the set including various elements of the courting ritual and various modalities [34][35][36][37][38][39][40][41][42][43]. The results have shown that all chromosomes, including the sex chromosomes, contribute to the speciesspecific variation in mating song elements, while the set of chromosomes or loci associated with the differences vary among related species pairs. The variation due to interspecific differences displayed no predominant association with the sex chromosomes [44].
Copulation, as a final stage of mating behavior, determines the efficiency of the insemination of the female and the contribution of the male genotype to the gene pool of the next generation. The copulation efficiency depends on the copulation duration and the insemination reaction. Male accessory gland secreted proteins (Acps), which evolve at a far greater rate as compared with hemolymph proteins [45][46][47][48][49][50], play a crucial role in the latter. 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 heterospecific crosses of D. virilis and D. americana is associated with three D. americana and two D. virilis recessive autosomal genes [51]. The finding agrees with the idea that isolating barriers arise rapidly 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 [52] observed that in crosses between the closely related species D. melanogaster and D. simulans, mating behavior was directly associated with the species-specific 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 final stage of tight genital coupling was 2-5 times shortened in heterospecific crosses, while the preceding unstable coupling stage was prolonged. The findings indicate that female insemination may be incomplete in heterospecific crosses and that a particular mechanism exists to allow mating flies to be sensitive to the specificity 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-specific 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 [53][54][55][56]. 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 [57]. However, the shape of the epandrium, which is the posterior lobe of the genital arch, has been used as the only model to study the genetic basis of species-specific 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 species-specific 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 interspecific 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.

Methods
Strains D. lummei 200 (Serebryanyi Bor, Moscow) and D. virilis 160 were from the Drosophila collection of the Koltzov Institute of Developmental Biology. The latter strain carried the following recessive autosomal markers: broken (b) on chromosome 2, gapped (gp) on chromosome 3, cardinal (cd) on chromosome 4, peach (pe) on chromosome 5, and glossy (gl) on chromosome 6. The first two markers are expressed as gaps in the second transverse and L2 veins, respectively; the other markers determine different eye colors, which are well-identifiable visually. The eye color in D. virilis wild-type is middle red wine, and it is darker than the eye color of D. melanogaster wild-type. The recessive mutation cardinal generates in homozygotes of D. virilis vermilion-like eyes. The recessive mutation peach generates in homozygotes of D. virilis yellowish pink eyes. The combination of homozygotes cardinal and peach produces reddish-brown eyes that darken to brown with age. The most complicated eye mutation is the recessive mutation glossy that produces eyes with the surface with an oily luster surface in reflected light and irregular ommatidia. The eye color is slightly redder and brighter. The combination of a homozygote for glossy with the other mutations or their combinations leads to the appearance of an oily sheen effect and a violation of the ommatidia shape. F 1 hybrid males were obtained from direct crosses of D. virilis females with D. lummei males and reciprocal crosses of D. lummei females with D. virilis males. Also, F b males were obtained by backcrossing hybrid males with D. virilis females. The crossing scheme is shown in Fig 1. 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 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 flies per tube. All males taken for the copulatory organ dissection were virgins and seven days old.

Analysis of morphological structures
The mating organ was dissected using thin steel needles in a drop of water under a binocular microscope at a magnification of 12 × 8. Preparations were incubated in boiling 10% NaOH to remove adipose structures.  As opposed to discretely identifiable phenotypes, such as eye color or vein gaps, the object of quantitative genetics is phenotypes that vary continuously: every character that can be measured (the length or width of the organ or its part) or counted (the bristle number). Frequently, quantitative genetics is referred to as the genetics of complex traits and is based on a model in which many genes influence a trait in a complicated way when they interact within themselves (allele interactions: dominant-recessive relationship of alleles), between themselves (epistatic relationships), and with environmental factors (genotype-environmental interaction) as the model usually accounted for non-genetic factors. According to the current findings, the genetic underpinning of such complex traits is even more complicated as noncoding DNA sequences like enhancers, silencers, and insulators affect the quantitative variability or heritability of quantitative traits. Moreover, such non-genetic effects as epigenetic factors (CpGislands) may also affect quantitative trait variability. There are different approaches to deal with morphology (traditional morphology when a complex trait is described by a set of measurements or counts, or geometric morphometrics when a set of landmarks or semi-landmarks describes a complex trait) at the stage of collecting data and with the total variability of a complex trait at the stage of statistical analysis. The latter includes a set of different statistical methods beginning from the methods invented by Fisher and Wright based on the decomposition of the total variance, like the analysis of variance, etc., and ending with methods of multivariate statistics and methods of data dimensionality reduction.
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 magnification of 300-500×. The sagittal view of the phallus was conventionally divided into four areas: the aedeagus body, gonites, apodeme, and hook. 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 1. The declinations of the hook 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 intraspecific 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". Trait dominance in interspecific crosses was evaluated by comparing the trait variance in fly samples from the parental and hybrid strains; backcross samples included only males with a fully heterozygous autosome set. Differences were tested for significance by multivariate ANOVA; groups displaying homogeneity of variance were identified for each trait by post-hoc comparisons, using the Gabriel test and Tukey's HSD test with the Bonferroni correction.
The effects of the sex chromosomes, autosomes, paternal genotype, and their interactions on the traits were evaluated using samples from the parental, F 1 , and F b populations, including both homozygotes and heterozygotes for the autosome set. A linear combination of weighted genetic factors was used as a model. The formal model is as follows: 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 coefficients 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 chromosome on the recessive autosomal alleles; Y (or X)!AUT(dom.epist), the epistatic effect of the Y (or X) chromosome on the autosomal dominant alleles; and P!X, the possible nongenetic 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 nongenetic effects were expected only for a combination of this factor with the offspring genotype. The genetic factors were treated as categorical variables; their significance 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 significance 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 Distance from the MP1 axis to the lower outline as measured 75% of the MP1 length away from point 1

7
Aedeagus height measured 75% of the MP1 length away from point 1

8
Distance from the MP1 axis to the lower outline as measured 50% of the MP1 length away from point 1

9
Aedeagus height measured 50% of the MP1 length away from point 1

10
Aedeagus height measured 25% of the MP1 length away from point 1

11
Distance from the MP1 axis to the lower outline as measured 25% of the MP1 length away from point 1

12
Aedeagus height at the base Distance from the MP1 axis to the paramere outline as measured at the base

24
Paramere length at the ventral curvature

25
Paramere length in the central part (at the MP23 midpoint)

26
Distance from the MP23 midpoint to the outline of the basal part of the paramere (characterizes the paramere bending angle)

27
Paramere height in the central part (at the MP26 midpoint)

29
Apodeme width measured 25% of the MP28 length away from point 1

31
Apodeme width measured 50% of the MP28 length away from point 1

33
Apodeme width measured 75% of the MP28 length away from point 1

34
Apodeme bending parameter (75% of the MP28 length away from point 1) https://doi.org/10.1371/journal.pone.0244339.t001 indicator variables were collinear. The values of the corresponding vectors were summed to combine several effects in one pooled effect. 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. 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 to allow for possible spermatogenesis defects and subsequent epigenetic marking of the chromosomes in interspecific hybrids. The resulting matrix of genotype-dependent vector values is shown in Table 2.
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 first two components were used in visualization.
All stages of the study are reflected in the pipeline of all analyses (Fig 3).

Results
Following the pipeline (Fig 3) results are presented in the corresponding order with a short introduction to each analysis describing the aim of the particular type of statistical analysis and the significance of the results obtained.

Correlated variability of traits of the copulatory organ shape
Traits that provide one function are linked by that function, and they should be correlated among them. The aim of the factor analysis is to reveal such correlation pleiads or, in other words, sets of correlated traits and the number of sets corresponding to the number of new latent traits that are called Factors.
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 factor analysis. Seven factors were taken based on the scree plot and together accounted for 64.4% of the total variance ( Table 3). Characteristics of the variation in the traits under study are shown in supplementary materials (S1 Table).
Based on the traits with the greatest factor loadings, the factors were classified as follows.  (Table 3, highly specific).

Generalized estimate of the dependence of the revealed variability of the phenotype of F 1 and F b hybrids on the genotype
In our case, redundancy analysis (RDA) was carried out to relate the set of variables that determine the factor structure values (a set of explanatory variables) to the variables that The eigenvalue and percent variance accounted for by a respective factor are shown in parentheses. Traits included in factor loadings with weights higher than 0.5 are shown. � indicates the traits that had weights higher than 0.7.
characterize the genotype (a set of response variables). RDA is a direct gradient analysis technique which summarizes linear relationships between components of response variables that are "redundant" with (i.e. "explained" by) a set of explanatory variables. To do this, RDA extends multiple linear regression (MLR) by allowing regression of multiple response variables on multiple explanatory variables (Fig 4). A matrix of the fitted values of all response variables generated through MLR is then subject to principal components analysis (PCA). In the analysis, the ordination of the genotypes and factors was performed in the space defined 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 seen from Fig 4, a two-dimensional space is defined by the orthogonal factor pair Fa3-Fa4 (X) and Fa6, Fa5-Fa2, Fa7 (Y). The genotypes showed the following arrangement in the two-dimensional factor space: +X: X Vi Y Vi Aut Vi Vi/Lu, X Vi Y Lu Aut Vi Vi/Lu, +Y(− X):   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-specific 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 4. The logic of estimating the degree of dominance for a trait has been described previously [58]. 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 significant difference (ns) in evaluating the degree of dominance. The variants f1(  Table 2. The IMPs that showed significant correlations were grouped. We describe the most general results of the analysis of variance. First, traits at which the D. virilis phenotype dominated (43) prevailed over traits with the dominance of the D. lummei phenotype (32) and traits with intermediate dominance (26) in all four samples. The dominance of the same phenotype regardless of the sex chromosome combination was observed for only 5 (IMPs 7, 16, 27, 28, and 30) out of the 30 IMPs included in the analysis, suggesting a substantial role of the autosome combination and the identity of the male parent in trait 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. 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-specific 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 the following. Firstly, 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 (Table 4, category "ns"). It is of interest to note that the most remarkable difference in the total number of traits showing dominance of one of the parental phenotypes was observed in backcross flies heterozygous for the autosomes. Thus, F b X Vi Y Vi males had eight traits at which the D. virilis phenotype dominated and ten traits at which D. lummei phenotype dominated, while F b X Vi Y Lu males displayed an opposite pattern and had ten and five such traits, respectively.
As expected, the number of traits at which the D. virilis phenotype dominated increased in backcross males homozygous for the D. virilis autosomes (S2 Table); the set included 19 traits in F b X Vi Y Vi males and 23 traits in F b X Vi Y Lu males. Substitution of the Y chromosome changed the character of dominance at nine traits, of which seven again demonstrated an adverse effect of the species identity of the Y chromosome on the phenotype.

Factor
Sign Traits are grouped according to their maximal weights in the factors extracted isolated by the maximum likelihood method (Table 3) Factorial MANOVA was carried out to evaluate the effects of the sex chromosomes, autosomes, paternal genotype, and their interactions on the phenotypic traits. Significant effects were observed for each of the four factors taken alone and the interaction of the Y chromosome and autosomes (S3 Table).
To study the effect on particular phenotypic traits for each of the factors, factorial MAN-OVA was performed with a forced incorporation of all four predictors and the effect of the ChrY � Aut interaction. The results are summarized in S4 Table. The majority of the traits showed an association with the identities of the X chromosome and autosomes at a significance level p < 0.05; half of the traits were affected by the Y chromosome and the species identity of the Y chromosome (a male parent identity); and one-third of the traits, by the interaction of the Y chromosome and autosomes.
Criteria for the size of the effect show the greatest influence of the autosomes on the total variability of the traits, and the least influence of the Y-chromosomes ( Table 5).
Note that a dependence on the hereditary factors was not confirmed for the apodeme width-related traits (F1). The apodeme declination and curvature (F2) depend on the identities of the Y chromosome and the species identity of the Y chromosome (a male parent). Distinct species-specific 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 hook and the outline bend over the hook (F6) showed a similar dependence on the Y chromosome, autosomes, and a male parent identity; an additional dependence on the X chromosome is specific 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.
An independent effect of the Y chromosome is expected to be insignificant based on the composition and functional activity of its coding sequences. It is possible to assume that interactions of the two factors play a crucial role in the phenotype.

The role of the components of variability and their combinations in the inheritance of traits of the shape of the male copulative organ
Post-hoc tests were used to evaluate the probability of the formation of homogeneous groups in trait variance according to a published model [59] (S5 Table) to detail the strength and direction of the effects exerted by the hereditary factors and their combinations on the dominance of the D. virilis phenotype. Criteria for the size of the effect of these factors show an approximately equal contribution of each of them to the total variability of the traits (S6 Table). The equal contribution indicates similar participation of the sex chromosomes, Chr X-chromosome X, Chr Y-chromosome Y, AUT-autosomes, Pmale-nongenetic effect of the parent species (D. virilis or D. lummei), which Y chromosome was inherited by a male; Wilks' α-a test statistic that is reported in results from factorial MANOVA, F-Fisher's F-statistic, df-degrees of freedom, p-probability value, Partial eta squared is the ratio of variance associated with an effect, plus that effect and its associated error variance, λ-one of the two parameters of the noncentral chisquare distribution, observed power is an estimate of the power of a study based on the observed effect size in a study. https://doi.org/10.1371/journal.pone.0244339.t005 autosomes, and the species identity of the Y chromosome (a male origin factor) in the formation of the analyzed interactions. The results obtained for each particular trait are provided in Supplementary; data on latent and highly specific traits are summarized in Table 6. Table 6. Paired permutation test with the Bonferroni correction to evaluate the optimal genotype partitioning into groups homogeneous in latent variables and highly specific variables (up to 100 000 iterations in each case). Designations of the factors-as in Table 3 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 6.

X!AUT (dom.epist) +X +AUT (dom) [
The results of analyzing how the primary and latent traits varied under the influence of the hereditary factors show a dependence on the hereditary factors that was observed for the majority of the latent traits and the corresponding primary traits with a high commonality. Trait F1, which includes traits related to apodeme width, was an exception; i.e., its dependence on any of the hereditary factors was not confirmed.
The epistatic effects exerted by the conspecific sex chromosomes and by the interaction between the Y chromosome and dominant alleles of the D. virilis autosomes positively affect the expression of the D. virilis phenotype in the majority of the morphological structures (F2 and F4-F7, IMP 3, 18, 19, 24 with high specificity).
A significant influence of the additive effect of recessive autosomal alleles (AUT(add)) was demonstrated for latent traits F3 and F5-F7 and confirmed for the majority of the corresponding primary traits and highly specific IMPs 3,18,19,24,26,and 27. Traits of the genotypes that were heterozygous for the autosomes and had an intermediate value of the indicator variable displayed differentiation by the dominance of one of the parental phenotypes. Wherein, latent trait F3, which characterizes the most distinct species-specific traits, and highly specific IMPs 3, 18, 26, and 27 showed the dominance of the D. lummei phenotype in heterozygotes for the autosomes. The additive effect of the autosomes significantly influences the dominance of the D. virilis phenotype and incorporates a dominant component that determines the expression of one of the parental phenotypes, depending on the trait in question. The dominant effect of the D. lummei autosomes is manifested for pronounced species-specific traits.
The epistatic effects of the sex chromosomes on recessive autosomal alleles (factor Y!AUT (rec.epist)+X!AUT(rec.epist)) significantly influence latent traits F3, F5, and F7 and highly specific primary IMPs 3,18,19,23,24,26,and 27. Wherein, the epistatic effects that the D. virilis sex chromosomes exert on recessive alleles of the conspecific autosomes lead to the dominance of the D. virilis phenotype.
A combined effect of autosomal dominant alleles, the X chromosome, and their epistatic interactions (factor X!AUT(dom.epist)+X+AUT(dom)) was confirmed for the majority of the primary traits having a high commonality; latent traits F2-F4, F6, and F7; and highly specific IMPs 19, 23, 24, and 27. However, the direction of dominance differed among the traits. Latent traits F3 and F4 showed a dominance of the D. lummei phenotype in males (genotype F1 X Lu Y Vi ), suggesting a predominant effect for the D. lummei X chromosome and autosomes. In contrast, traits related to the apodeme bend and declination (F2) showed the possibility of superdominance of the D. virilis phenotype in the intermediate genotype and a leading role of dominant alleles of the D. virilis autosomes. A similar effect of D. virilis autosomal dominant alleles was confirmed for latent trait F7 and highly specific IMPs 19 and 24 at high significance. Latent trait F6 and highly specific IMP 23 displayed incomplete D. virilis phenotype dominance of genotype X Lu Y Vi Aut Lu/Vi .
The influence of the male parent identity on the effect of the X chromosome and autosomal dominant alleles (factor P!X+P!AUT(dom)) was demonstrated for latent traits F2, F3, F4, and F6 and highly specific primary IMPs 24 and 28. Latent traits F3 and F6 nonlinearly depended on the interaction of the male parent identity and the X chromosome and autosomes. The dominance of the D. lummei phenotype was observed in the genotype having the intermediate value of the indicator variable for latent trait F2; and highly specific primary IMP 28. A significant effect of this factor was not demonstrated for latent traits F1, F5, and F7. In contrast, individual primary traits incorporated in the latent factor displayed, in the intermediate genotype, the dominance of the D. virilis phenotype in the case of F1 and F5 and the D. lummei phenotype in the case of F7. With F7, the mean values gradually increased in the order of genotypes with indicator variable values 1-0-2, suggesting superdominance of the D. lummei phenotype.
The effect of the father's origin on the action of the Y chromosome and recessive alleles of autosomes (factor AUT (add) _P + Y_P) is shown for all latent variables, except the first, and all primary characters with a high characteristic, except IMP 26. The Y chromosome, under the influence of the father's origin, makes a maximum contribution to the phenotype manifestation for the latent traits F3 and F4. Wherein, for the F3 latent trait, a mixed effect of the Y chromosome and autosomes is observed, leading to the dominance of the D. lummei phenotype in genotypes with intermediate values of indicator variables. The primary traits defining the latent trait F4 confirm the noted dependence. The result of the non-random distribution of genotypes into heterogeneous groups for the latent trait F2 and the primary trait with high characteristic IMP 27 does not allow us to unambiguously relate the distribution data to the influence of the estimated factor. In both cases, a key role is played by the independently clustered male genotype from F1 X Lu Y Vi Aut Lu/Vi , which has the D. lummei X chromosome and, accordingly, all its effects, which suggests a possible contribution of this chromosome. The primary characters that make up F2 also show different variants of dependence on autosomes and Y-and X-chromosomes. Latent traits F5, F6, and F7, as well as primary traits with a high characteristic IMP 3, 18, 19, 23, 24, and 28, show a significant contribution of autosomes modified by the effect of the father's origin.
Thus, the estimates confirm that epistatic interactions of the sex chromosomes and autosomes and effects of the male parent origin from interspecific crosses influence the expression of D. virilis species-specific 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 efficiency of female insemination in insects depends on how well the female and male genitalia match each other [52]. Incomplete insemination must make repeated mating more likely, lead to displacement of sperm from the previous mating, and facilitate efficient selection against the incomplete insemination-associated genotype in Drosophila. Hoikkala and colleagues [60] have analyzed the within-population variation of mating duration in D. montana and associated the duration of the first 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 significant parts of the Drosophila copulatory system with the duration of mating [61].
Both the copulation duration 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 [62,63]. 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 specifically expressed as predictors of sexual reproduction efficiency. When a new adaptation forms and the adaptive norm changes rapidly, selection changes to directional or disruptive, but 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.
Two models describe the variation maintained by stabilizing selection equally well: a jointeffect model and an infinitesimal model. In the former, the variation is associated with the effects of moderate-frequency alleles, which are nearly neutral in terms of fitness and substantially influence the trait in question. In the latter, the variation involves many genes that each exert a minor effect on fitness and act additively. Modeling of comprehensive experimental data has shown that well reproducible results are obtained with both of the models, neither of them is preferable to the other [63]. 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 [64][65][66]. 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 [67,68]. Moreover, opposite epistatic effects may occur between tightly linked genes in a QTL [69] 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 significant traits has other genetic dynamics. Fisher [70] has noted that adaptation is not selection, meaning fitness-based selection that eliminates the least fit individuals. Adaptation is characterized by the effect of positive selection, which facilitates the progress of a population to an optimal phenotype. Based on Maynard Smith's model of adaptive walk in sequence space [71] and Gillespie's theory of fixation of rare beneficial alleles [72,73], it is thought now that adaptations arise in bursts alternating with periods of slow evolution [74]. Adaptation is associated with a few effects, and selection rapidly decreases in strength with each subsequent step; i.e., an exponential distribution is characteristic of the selection coefficients of consecutively fixed mutations [75]. An allele that expands the ecological niche for the species will be fixed rapidly in a new QTL. Experiments on QTL mapping support the conclusions based on model analyses [76][77][78][79]. It is safe to say that, along with additive effects, dominant interactions of alleles gain principal importance. A detailed analysis of the genetic architecture and ontogenetic mechanisms of species-specific traits in two Labeotropheus species from the Malawi Lake has made it possible to evaluate the role that regulatory, or epistatic, interactions of genes involved in the Tgfβ signaling pathway play in the formation of foraging adaptations [80].
Additive, dominant, epistatic, and nongenetic components of variation underlie the phenotypic differences in our experiment. The experiment was not designed to determine the variance components precisely. 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 the nongenetic 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 the sex chromosomes and autosomes were specific 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 autosomal dominant 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. One of these was the dominance of the D. lummei phenotype in latent trait F3, which combined the most distinct species-specific traits; certain primary traits that were incorporated in latent traits F4 and F5 and similarly showed species specificity [81]; and paramere width (IMP 27). The dominance of the D. virilis phenotype was characteristic of latent trait F7; several primary traits were 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 [82], 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, the 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 nongenetic influence on autosomal dominant 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 confirms that the dominant component of the autosomes contributes to the species-specific variation.
The effect that epistatic interactions between the X and Y chromosomes exert on speciesspecific traits was confirmed for half of the primary traits and latent traits F3 and F5, which incorporate taxonomically significant 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 conspecific 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 [83] 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 interspecific hybrid male sterility as a model of interactions between autosomes and the 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 [84,85], 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 heterospecific X and Y chromosomes leads to decondensation of Y-chromosomal heterochromatin, dramatically distorts the expression of autosomal genes, and causes sterility [24,86]. 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 confirmed 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 . 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-specific male courtship song pattern in Drosophila species of the montana phylad [38]. There is evidence that the sex chromosomes exert a significant regulatory effect on the expression activity of autosomal genes. Expression of X-chromosomal genes in hemizygous D. simulans males depends on trans regulations to a substantial extent [87], 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 [88].
The important role that the Y chromosome plays in regulating the phallus shape-related traits in interspecific hybrids is possible to associate with trans-regulatory activity directly. 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 [89][90][91][92]. Additive variation of the ten other genes should make a vanishingly small contribution to the phenotype at the 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-specific lethal proteins/roX1,2 RNA complex, which is responsible for the dosage compensation of X-linked 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 [93]. The viability was low in males with the paternal Y chromosome and high in males with the maternal X chromosome. The result indicates that the Y chromosome modifies dosage compensation through roX1, roX2-mediated modification of heterochromatin [94] and/or recognition of the X chromosome by the entire male-specific lethal proteins/roX1,2 RNA complex. The results support our finding of a substantial paternal effect, which acts independently in the majority of traits or combination with the X-chromosome effect in some other cases. Finally, Lemos and colleagues [92,95] have directly demonstrated the regulatory activity of the Y chromosome in a study of differential genome expression activity under the Y-chromosome influence.
Traits of the apodeme, a muscle attachment site, and an internal part of the copulatory system, are the least associated with evolutionary variation and chromosome effects. It is noteworthy that chromosome effects are mediated by the nongenetic effect of male parent identity in the apodeme traits. Note that significant nongenetic effects were similarly observed for all other groups of correlated traits. A change in the type of dominance of half of the analyzed characters in males with an identical genotype, differing only in the status of the father's origin, is the basis for identifying these effects as an independent factor. Neither the mitochondrial genome nor the genotype for the nuclear genome per se is related to this effect. It can be assumed that it may be associated with the Meiotic Sex Chromosome Inactivation mechanism (MSCI) or, in a broader sense, meiotic silencing of unpaired chromatin/DNA (MSUC/MSUD) [21,96], which leads to heterochromatinization and 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 have been the matter of extensive discussion [21,23,96,97]. MSCI is typically compensated for by the higher transcriptional activity of male-specific genes in Drosophila because their strong promoters are overrepresented in the X chromosome [98,99].
On the contrary, the delay in forming of heterochromatin blocks is often associated with sterility in hybrids from interspecific crosses and the death of hybrid progenies; defects in mitotic chromosome segregation are observed hybrids during their embryo development [24,25,96]. Meiotic silence of unpaired heterochromatin and delayed condensation of heterochromatin blocks, leading to lethality during mitosis, are opposite processes, but the marking of such blocks may be of a general nature. Epigenetic signatures are preserved in chromosomes at post meiotic stages and may be inherited through generations [92,100]. Although the regulation of genome-wide epigenetic states is often associated with the Y chromosome [92,101,102], a formal role of the interspecific hybrid status is noteworthy in our case, the manifestation of which may be associated with a violation of the diverging chromosomes leading to a violation of the synapse of divergent chromosomes and causing the formation and further preservation of noncanonical heterochromatin regions with altered expression activity. Thus, Argyridou and Parsch note a decrease in the average expression activity of genes of the X chromosome in somatic tissues compared with autosomal genes [103]. Although the authors conclude that this decrease is associated with the X chromosome being depauperate in genes with very high expression, even after removing overexpressed genes from the pool of analyzed genes, significant differences in the expression activity persist for intestinal tissues in males and females. The X chromosome enrichment with male-biased genes has also been shown for gene expression in the Drosophila brain [104], confirming the noted effect of selection as a factor in the evolution of sexual dimorphism.
Another potential mechanism should be noted, the action of which determines the effect of the father associated with the action of micro-RNA. Although the work was performed on a transcriptome from human tissue, this is the first study to show the accumulation of specific intergenic repeat-derived RNAs in spermatozoa [105], which can influence embryonic genome activation and determine differentiation processes at subsequent stages of development.

Conclusions
All types of genetic interactions between the chromosomes were revealed in our study, and most of the shape traits of the male copulatory organ demonstrated the dominance of the D. virilis type. Epistatic interactions between the X and Y-chromosomes, as well as between the sex chromosomes and autosomes, affect species-specific traits. The two sex chromosomes shifted phenotype dominance in the conspecific direction.
The effects of sex chromosomes on the shape of the male mating organ is comparable with the effect of autosomes, even though the effect sizes of the former were inferior to the latter. The X-chromosome contains several times fewer genes compared with autosomes, whereas the Y-chromosome contains no more than 20 genes. This insignificant number of genes is incommensurable with the effect sizes of the sex chromosomes, which leads to the assumption that natural selection supported the accumulation of structural and regulatory elements on the sex chromosomes. Also, it can be assumed that sexual selection for specific male-biased genes associated with the shape of the male mating organ might prevent the loss of these genes on the X chromosome.  Table 2. The chromosomes and paternal genotype are indicated in the following order: X chromosome, Y chromosome, autosomes, male parent identity. (DOCX) S1 Table. Variation of the morphometric traits in the shape of the male copulatory system. The sample size is indicated in parentheses. M, mean; σ, standard deviation. Superscripts in the first column indicate that the trait is incorporated with a high weight in the respective factor structure.  Table. Effects of the sex chromosomes, autosomes, and parental genotypes on trait expression. Traits are grouped according to their maximal weights in the respective factors (Table 3). F, Fisher's test; p, the significance of effects of independent variables, including X, the X chromosome; Aut, the autosomes; ♂P, the paternal genotype; and ChrY � Aut, a combined effect of the Y chromosome and autosomes. HC, a group of the traits that were not incorporated in the factor structures with weights higher than |0.5|. Significance values p < 0.05 are in bold. In each group of traits determining the respective factor structure, the lowermost row shows the estimated effects of the independent variables on the given factor. (DOCX) S5 Table. Permutational ANOVA to separate the genotypes into groups homogeneous in trait values. (+), a significant effect of a hereditary factor that is in linear relationship with the given genotype classes; (+ � ), a significant effect of a hereditary factor that is in nonlinear relationship with the given genotype classes; D Vi (Lu) , dominance of the D. virilis (D. lummei) phenotype in genotypes with the intermediate value of the indicator variable; P Y(Aut) , an epigenetic effect of the male parent identity on the Y chromosome (autosomes); the subscripts Aut, X-chr, and Y-chr indicate that the autosome set, X chromosome, or Y chromosome determines phenotype dominance; F#, latent trait number as in Table 3; F# � , permutation test values for the factors (latent traits) used as variables; CH, highly characteristic trait; FDR, false discovery rate for the primary traits.