Sex-specific and sex-chromosome regulatory evolution underlie widespread misregulation of inter-species hybrid transcriptomes

Abstract When gene regulatory networks diverge between species, their dysfunctional expression in inter-species hybrid individuals can create genetic incompatibilities that underlie the developmental defects responsible for intrinsic post-zygotic reproductive isolation. Divergence in cis- and trans-acting regulatory controls evolve despite stabilizing selection on gene expression, being hastened by directional selection with adaptation, sexual selection, and inter-sexual conflict. Dysfunctional sex-biased gene expression, in particular, may provide an important source of genetic incompatibilities, with more severe misregulation expected for the heterogametic sex. Here, we characterize and compare male and female transcriptome profiles for sibling species of Caenorhabditis nematodes, C. briggsae and C. nigoni, and allele-specific expression in their F1 hybrids to deconvolve features of expression divergence and regulatory dysfunction. Despite evidence of widespread stabilizing selection on gene expression, we find broad misregulation of sex-biased genes in F1 hybrids that is most pronounced for the X-chromosome, supporting a “large-X” effect, and that counters expectations by disproportionately affecting hybrid females. Hybrid male misexpression, however, is greater in magnitude, with spermatogenesis genes especially prone to high divergence in both expression and coding sequences that may explain elevated sterility of hybrid males, consistent with “faster male” and “fragile male” models for Haldane’s rule. Regulatory and coding divergence overall correlate only weakly, however, and downregulation of male-biased genes in females implicates trans-acting modifiers in the evolutionary resolution of inter-sexual conflicts. This work identifies important differences between the sexes in how regulatory networks diverge that contributes to sex-biases in how genetic incompatibilities manifest during the speciation process. Author’s summary Many mutations that affect traits as species diverge do so by altering gene expression. Such gene regulatory changes also accumulate in the control of static traits, due to compensatory effects of mutation on multiple regulatory elements. Theory predicts many of these changes to cause inter-species hybrids to experience dysfunctional gene expression that leads to reduced fitness, disproportionately affecting the sex chromosomes and sex-biased gene expression. Our analyses of genome-wide expression data from Caenorhabditis nematode roundworms support these predictions. We find widespread rewiring of gene regulation despite extensive morphological stasis, and conserved overall expression profiles, that is a hallmark of these animals. Misregulation of expression in both sexes is most severe for genes linked to the X-chromosome, sperm genes show distinctive signatures of divergence, and differences between the sexes in regulatory evolution implicate resolved historical sexual conflicts over gene expression. This work clarifies how distinct components of regulatory networks evolve and contribute to sex differences in the manifestation of genetic incompatibilities in the speciation process.

Introduction necessarily be abundant on the X-chromosome (40-42). Thus, distinguishing between abnormal 114 expression in hybrids for X-linked genes overall and for sex-biased autosomal genes is 115 important to decipher the genetic mechanisms that underpin Haldane's rule in particular and 116 hybrid dysfunction in general. 117 118 Caenorhabditis nematode roundworms provide an especially tractable system to study 119 speciation genomics (34). The growing number of Caenorhabditis species known to science 120 conform to the biological species concept, with a few cases where sibling species can produce 121 some viable and fertile adult hybrid offspring (43,44). The C. briggsae × C. nigoni species pair is 122 one such case, where recent divergence (~3.5 Ma (45)) allows them to form hybrids of both 123 sexes. Haldane's rule is fulfilled: F1 male hybrids exhibit both greater infertility and inviability in 124 crosses of C. nigoni females to C. briggsae males, with nearly total hybrid male inviability in the 125 reciprocal cross (46-48). The disproportionate loss of genes with male-biased expression in the 126 C. briggsae lineage has the potential to disrupt genetic networks in a sex-biased way to 127 disproportionately compromise hybrid male fitness (45,(49)(50)(51). In an effort to identify and 128 locate incompatibility loci between these two species, Bi et al. (52) associated hybrid male 129 inviability and sterility with most of the X-chromosome, suggesting a large-X effect. Moreover, 130 analysis of X-chromosome introgression lines revealed that X-autosome incompatibility 131 involving misregulation of the 22G class of small RNAs led to down-regulation of 132 spermatogenesis genes as a contributor to hybrid male sterility, in addition to autosomal 133 factors responsible for hybrid inviability (29,53). Here, we analyze mRNA transcriptome 134 expression for each sex from each of C. briggsae, C. nigoni, and their F1 hybrids. Using ASE 135 profiles, we then quantify cisand trans-acting regulatory causes of expression divergence to 136 link genomic change to sex-biased expression, chromosomal features, and hybrid dysfunction. 137

138
Extensive expression divergence between species and between the sexes involve 139 the X-chromosome 140 Each species and sex show distinctive overall transcriptome profiles that are further 141 distinguishable from each sex of F1 inter-species hybrids (Fig. 1A). In the contrast of pure C. 142 briggsae and C. nigoni transcriptomes, we found more differentially expressed genes for 143 females and hermaphrodites, hereafter jointly referred to as females for brevity, than for 144 males. For the 12,115 one-to-one orthologs analyzed, females had a total of 66% (7,903) of 145 genes differentially expressed between species, compared to 53% (6,391) for males (Fig. 1B). 146 The X-chromosome shows the most extreme differences in expression between species, with 147 both males and females showing significantly higher expression of X-linked genes in C. briggsae 148 than in C. nigoni (Fig. 1C). Autosomes, by contrast, showed greater abundance of genes with 149 higher expression in C. nigoni, albeit only significantly only for Chromosome I in females (Fig. 150 1C; Fisher's exact test; P < 0.05). 151 152 Within each species, approximately half of genes exhibited significant sex biases in expression, 153 which, in turn, are roughly evenly divided between male-and female-biased expression: 28% 154 male-biased genes in C. briggsae (3,353 genes), 26% male-biased genes in C. nigoni (3,188 155 genes), 26% female-biased genes in C. briggsae (3,205 genes), and 24% female-biased genes in 156 C. nigoni (2,878 genes). Male-biased genes did not exhibit strong enrichment for any 157 chromosome in either species, whereas female-biased genes in C. briggsae were 1.6-fold 158 enriched on Chromosome I and 3.3-fold depleted on the X-chromosome (Fig. 1E), consistent 159 with previous studies in Caenorhabditis. The X-chromosome is significantly enriched, however, 160 for non-sexually differentiated genes in both C. briggsae and C. nigoni (Fisher's exact test, P < 161 0.05) (Fig. 1E). Overall, expression profiles for C. briggsae hermaphrodites were masculinized 162 relative to C. nigoni females (Fig. 1A), consistent with hermaphrodite production of both sperm 163 and oocytes in an otherwise morphologically female soma. 164 165 Expression dominance in F 1 hybrids differs between males and females 166 We contrasted expression profiles of F1 hybrids with their parent species to infer the expression 167 inheritance of genes, i.e., to identify genes that exhibited additive, dominant (C. briggsaeand 168 C. nigoni-like expression), or transgressive (overdominant and underdominant) expression 169 patterns for each sex ( Fig. 2A). Gene sets with distinct expression inheritance profiles revealed 170 substantial differences between the sexes in terms of expression distance (Fig. 2B), number of 171 genes (Fig. 2C), and enrichment across the genome (Fig. 2D). 172

173
The sexes differed most strikingly in their total number of transgressive genes (Fig. 2C). 174 Transgressive genes are thought to be associated with hybrid dysfunction, as they represent 175 expression phenotypes that are distinct from expression levels in either parent (28,54). Despite 176 our expectation of an especially high number of misregulated genes in F1 males due to their 177 pronounced sterility, we found that just 21% (2,552) of genes show transgressive profiles in 178 males (1,387 overdominant and 1,165 underdominant) compared to 55% (6,729) in females 179 (3,284 overdominant and 3,445 underdominant) ( Fig. 2A,C). However, we also found that 180 Euclidean expression distances from the centroid of expression space are higher on average in 181 males than in females for both overdominant (ordinary least-squared [OLS] regression, t = 6.28, 182 P < 0.001) and underdominant genes (OLS, t = 2.14, P < 0.05), indicating especially deviant 183 expression magnitude for those fewer transgressive genes in F1 males ( Fig. 2A,B). By contrast, 184 genes with simple dominance or additive expression in F1 males had consistently lower 185 expression distance from the centroid than did F1 females ( Fig. 2A,B; OLS, P < 0.05). These 186 results are consistent with our multidimensional scaling analysis (Fig. 1A) that showed shorter 187 expression distances for F1 males to parental males, in contrast to more dissimilar expression 188 profiles of F1 females to parental females. 189 190 F1 hybrids of both sexes had a relatively low percentage of genes with additive expression (6% 191 or 698 genes in females; 7% or 872 genes in males) (Fig. 2C). By contrast, approximately 30% of 192 genes expressed by each sex showed simple dominance, either matching C. briggsae or C. 193 nigoni expression (27% or 3,878 genes in females; 31% or 3,297 genes in males). Hybrids of 194 both sexes consistently showed a higher number of genes with expression dominance matching 195 C. briggsae (20% or 2,408 genes in males and 16% or 1,904 genes in females), compared to 196 expression dominance matching C. nigoni ( Fig. 2C; 12% or 1470 genes in males and 11% or 1393 197 in females). In both sexes, the ratio of genes with simple dominance matching either C. 198 briggsae or C. nigoni expression is close to 1:1 consistently across all five autosomes (mean 199 ratio = 1.15, sd = 0.13). In F1 females, the X-chromosome, however, was 4-fold biased towards 200 expression dominance matching C. briggsae (Fig. 2C,D). 201 202 Genes and traits with dysfunctional expression are often associated with the X-chromosome, 203 and differences between the sexes are expected due to Haldane's rule (24,26,27,38). Consistent 204 with these expectations, we found expression heritability profiles in F1 males and females to 205 vary across the genome, and to differ most conspicuously for the X-chromosome. The X-206 chromosome was enriched for underdominant genes in both F1 males and females (Fig. 2C,D; 207 Fisher's exact test, P < 0.05). The X-chromosome was also enriched for overdominant genes in 208 F1 males, whereas females had significant depletions of such genes on the X (autosomal 209 enrichments: V in males, I and III in females) (Fig. 2C,D). These data show clear differences in 210 expression heritability across chromosomes between the sexes and reflect distinct hybrid 211 expression dynamics between autosomal and X-linked genes. 212 213 Regulatory divergence is modulated by differences in cis and trans effects between 214 sexes 215 Identifying the spectrum of changes in cisand trans-acting regulators is important to 216 understand how selection influences the evolution of gene expression and its effects on hybrid 217 phenotypes. Correspondingly, we classified the types of regulatory changes and examined gene 218 misregulation in F1 hybrids for each sex. We found that most expression divergence between 219 species results from cis-only, trans-only, and enhancing cis-trans gene regulatory divergence 220 (Fig. 2E,F), consistent with other ASE studies in flies (16,23), mice (24), plants (17,18) and yeast 221 (19,25). Comparing regulatory divergence between sexes, we found double the number of 222 trans-only changes involving genes expressed in females (20% or 2,461 genes) compared to 223 males (10% or 1,230 genes) (Fig. 2G). However, the sexes showed a reciprocal pattern for genes 224 with cis-only divergence being more prevalent in males than in females (21% or 2,559 genes in 225 males; 17% or 2,089 genes in females). Genes with cis-only and trans-only effects were not 226 significantly enriched on any autosome for either sex, but genes expressed in males showed 227 strong enrichment on the X-chromosome for trans-only regulatory changes as well as a 228 depletion for cis-only changes ( Fig. 2H; Fisher's exact test, P < 0.05). While genes with C. 229 briggsae expression dominance were commonly associated with cis-only effects, especially in F1 230 males, genes with C. nigoni dominant expression tended to show trans-only effects, especially 231 for autosomes in F1 females (Fig. 3). In contrast, X-linked genes in F1 females rarely showed C. 232 nigoni expression dominance displayed by either cis-only or trans-only effects (Fig. 3). In 233 contrast, X-linked genes in both sexes have up to 3-fold higher proportion of genes showing C. 234 briggsae expression dominance more often associated with trans-only effects (432:106 C. 235 briggsae : C. nigoni in females; 329 : 274 in males) (Fig. 3). This skew is particularly notable for 236 F1 males given that they carry the C. nigoni X-chromosome, and therefore these genes would be 237 Consistent with the idea that cis and trans changes each play distinct roles in sex-biased 246 expression and sexual dimorphism (55), we find, on one hand, that trans-only changes in 247 females are more often associated with male-biased genes, and on the other hand, that female-248 biased genes show more cis-only regulatory changes in males (Fig. 4E) classified as additive to associate more often with significant cis-acting divergence in both sexes 260 (i.e., cis-only and enhancing cis-trans effects; Fig. 3; Supplementary Fig. S4). However, 261 expression additivity is not abundant in our analysis (Fig. 2C), suggesting that it is not a major 262 source of phenotypic dysfunction in hybrids of this system. 263

264
Stabilizing selection on expression level is thought to be common in many species, at the 265 molecular level enabled through the coevolutionary fine-tuning by changes to cisand trans-266 acting factors (7,12). cis and trans changes with opposing effects can interact epistatically in 267 hybrids to induce dysfunctional expression and allelic imbalance (21). Consistent with this idea, 268 we found that transgressive genes with overdominant effects in hybrids often are associated 269 with cis-trans regulatory changes in hybrids (43% or 539 genes in males and 42% or 1,312 genes 270 in females; Fig. 3). In F1 males, a higher fraction of changes was compensatory (26%) compared 271 to enhancing (17%). F1 females had the opposite pattern: a higher fraction of genes showed 272 enhancing compared to compensatory changes (26% versus 16%). In contrast to genes with 273 such overdominant expression profiles in F1 hybrids, underdominant genes consistently 274 exhibited more conserved regulatory controls in both sexes (53% or 577 genes in males and 275 50% or 1,541 genes in females identified as conserved or ambiguous), suggesting alternative 276 mechanisms of regulatory dysfunction. These results highlight how stabilizing selection can act 277 differently on sex-specific transcript abundance, leading to opposing cis and trans effects that 278 are dysfunctional in F1 hybrids. 279

280
To further assess the role that different regulatory controls play in the origin and maintenance 281 of divergent sex-biased expression, we contrasted expression heritability and patterns of cis-282 and trans-regulatory divergence for male-biased and female-biased genes (Fig. 4D,E). We found 283 that male-biased genes expressed in F1 females frequently show transgressive misexpression 284 (underdominant) or have expression levels matching those of the species with lower expression 285 ( Fig. 4A,D). In contrast, male-biased genes expressed in F1 males tend to show expression 286 dominance matching C. briggsae, regardless of which pure species had higher gene expression 287 ( Fig. 4A,D). 288 289 Interestingly, genes expressed in F1 males are more commonly underdominant when they 290 correspond to male-biased genes than to genes that have a female-biased expression pattern 291 (623 genes among male-biased genes vs 281 genes among female-biased genes), suggesting 292 male-biased regulatory networks are more fragile. This idea is also supported by male-biased 293 genes having a higher proportion of genes with enhancing or compensatory cis-trans changes 294 (835 genes in male-biased genes vs 398 genes in female-biased genes). By contrast, female-295 biased genes in F1 females are predominantly overdominant (Fig. 4D) and are more often 296 associated with cis and enhancing cis-trans changes, which suggest that female gene regulatory 297 networks can be more resilient to misexpression, which may translate to similar resilience for 298 traits such as fertility (46). 299 300 Faster regulatory and molecular evolution of male-biased and spermatogenesis 301 genes 302 Sexual selection and sexual conflict are predicted to drive faster rates of molecular evolution 303 and expression divergence (36,57,58). Consistent with these predictions, we found that male-304 biased genes have higher average expression divergence (OLS, P < 0.001) and faster rates of 305 molecular evolution (Ka/Ks, OLS, P < 0.05) than female-biased genes (Fig. 4B,C). Compared to 306 sex-neutral genes, however, the signal for faster sequence evolution was weak (Ka/Ks OLS, P = 307 0.3), despite remaining strong for elevated expression divergence (OLS, P < 0.001). 308 309 We observed the highest expression divergence as well as high average rates of molecular 310 evolution in the distinctive set of genes with male-biased expression, higher expression in C. 311 briggsae than C. nigoni, and with a species-by-sex interaction (M6) (Fig. 4B,C). The species-by 312 sex interaction in M6 indicated a masculinized expression profile for C. briggsae 313 hermaphrodites, implicating a role for them in sperm production. To test this idea, we looked at 314 C. elegans genes previously identified as spermatogenesis genes (42)  for "no change" and "ambiguous" categories). For F1 males, X-linked genes showing expression 325 divergence and dominance matching C. briggsae levels in hybrid males should be effectively 326 misexpressed, as regulation is controlled mostly by trans-acting factors coming from the C. 327 briggsae autosomal background that interact with the C. nigoni-derived X-chromosome, and 328 therefore represent misexpression (Supplementary Fig. S3; also see (26)). We observed a 1.4-329 fold enrichment, although not significant (Fisher's exact test, P = 0.31, odds ratio = 1.42), of 330 genes with these effects on the X-chromosome of hybrid males (Fig. 4H). Furthermore, 331 transgressive expression on the X-chromosome in F1 males should be driven by enhancing 332 and/or compensatory cis-trans changes, for which we observed significant enrichments on the 333 X ( Fig. 4H lower panel; Fisher's exact test, P < 0.05). Together, these results show that, although 334 there are fewer X-linked spermatogenesis genes compared to autosomes (Fig. 4G), they show 335 high levels of X-linked expression dysfunction (61%, 37 genes of 61 X-linked genes). In addition, 336 it suggests that modest number of key genes may dictate X-autosome incompatibilities and X-337 linked misregulation and misexpression with effects on hybrid male fertility. 338 Genome architecture and molecular evolution moderately affect regulatory 340 divergence 341 Given that protein-coding sequence evolution and gene composition vary non-randomly along 342 chromosomes in many Caenorhabditis species in association with the chromosomal 343 recombination landscape, we asked whether distinct chromosomal domains would also 344 associate with the degree of cis-regulatory divergence. We find higher molecular divergence 345 between the genomes of C. briggsae and C. nigoni in chromosomal arms compared to centers 346 in noncoding sequences upstream of protein-coding genes (Fig. 5A), in addition to protein-347 coding sequence divergence (Supplementary Fig. S5; also see (45)). These observations are 348 consistent with the idea of stronger purifying selection on mutations to genes and their cis-349 regulatory regions when linked to chromosome centers. Despite the elevated molecular 350 divergence in arm regions (Fig. 5A), we only found modest elevation of ASE divergence for 351 genes on arms, being strongest for chromosome V for both sexes (Fig. 5B). Moreover, we found 352 no significant differences in the magnitude of regulatory divergence, in either cis or trans, 353 between chromosome arms and centers ( year), including ~20% sequence divergence for synonymous sites, changes to genome size, and 380 disproportionate loss of short male-biased genes in C. briggsae since its transition in 381 reproductive mode to androdioecy (45,51). Despite this genomic divergence, hybridization 382 between these species yields viable and fertile F1 females, while hybrid males suffer complete 383 sterility and severe inviability depending on the cross direction (34,46). While in some systems, 384 such as fruit files and plants, hybrids can exhibit an expression bias towards the parental 385 species tending to show higher expression (16,18,23), we observe no such effect ( Fig. 1B;  386 Fisher's exact test, P value = 0.76). This symmetry in expression suggests that demographic 387 effects do not bias regulatory changes toward either increased or decreased expression in this 388 system, as could occur if regulatory changes fix more rapidly in species like C. briggsae with 389 lower effective population sizes. However, autosomes tend to have more genes with slightly 390 higher expression in C. nigoni whereas X-linked genes tend to have much higher expression in C. 391 briggsae. This chromosomal pattern implicates a disproportionate role for divergence of the X-392 chromosome in mediating misexpression in Caenorhabditis hybrids. 393 394 Despite our evidence of substantial expression divergence, 34% of genes in females (4,212 395 genes) and 47% of genes in males (5,724 genes) show no differential expression between 396 species. Stabilizing selection on transcript abundance is recognized as a common force acting to 397 conserve expression levels between species (8-10,12). Mechanistically, conservation of the 398 expression phenotype can occur, despite sequence evolution, when co-evolution changes both 399 cisand trans-regulatory elements: if a trans-acting mutation fixes due to a pleiotropic benefit 400 on other loci, selection would favor fixation of any subsequent compensatory mutation in cis 401 that returns expression to optimal levels at the focal locus (5,59). We find evidence of extensive 402 compensatory cis-trans divergence in gene regulation between C. briggsae and C. nigoni. Such 403 coevolution represents just one mechanism leading to "developmental system drift," in which differently and independently among diverging lineages. In hybrids, uncoupled regulatory 415 mechanisms from the two parental genomes can lead hybrids to experience misregulation and 416 therefore misexpression (6). Such a mechanism underlying misexpression could represent a 417 Dobzhansky-Muller incompatibility because genetic interactions untested by natural selection 418 will likely be detrimental (1). The clearest signal of misexpression in hybrids is the sharp 419 contrast in the fraction of sex-biased genes: ~90% in hybrids vs ~50 in each parental species 420 (Fig. 1D). Moreover, the usual depletion of female-biased gene expression from Caenorhabditis 421 X-chromosomes is even more extreme in F1s due in part to transgressive underdominance 422 effects and, unusually, the X-chromosome is highly enriched for male-biased expression in F1s 423 (Figs. 1E, 2C,D) (40,41). The strong downward misexpression (underdominance) observed for 424 the X in females, but not as strong for the X of males is likely to be responsible for this trend. In 48), so we expected more misexpression in hybrid males. In contrast, we find that it is hybrid 436 females that experience more extensive transgressive misexpression of genes across the 437 genome that exceed the expression extremes of either parental species ( Fig. 2A,C). Most 438 female-biased transgressive genes show overdominant misexpression whereas male-biased 439 transgressive genes tend toward underdominant misexpression (M8-M12 vs M3-M7, Fig. 4D). 440 Excluding spermatogenesis genes (most of which are male-biased, but also expressed in C. 441 briggsae hermaphrodites; Fig. 4F), it is plausible that cisand/or trans-regulatory changes 442 acquired after speciation favoring female-biased expression experienced selection to sustain 443 upregulation, behaving in an overdominant manner in hybrids. Indeed, overdominant genes 444 with cis-trans divergence in females have disproportionately evolved "enhancing" regulatory 445 changes (Figs. 2G, 3, 4E). The low magnitude of expression divergence among overdominant 446 and female-biased genes (Figs. 2B, 4B), however, together with the fact that hybrid females are 447 fertile, suggests that overdominant expression does not impact fitness as negatively as does 448 regulatory divergence that leads to underdominance in hybrids. 449 450 Interestingly, we find that regulatory controls suppressing or enhancing male-biased expression 451 in F1 females is largely due to trans-only regulatory changes, particularly among genes that 452 show both male-biased expression and expression divergence between the parent species (Fig.  453   4E). Many of these trans-acting regulatory changes tend to be more strongly associated with C. 454 nigoni dominant expression in females among autosomes, contrasting with cis-regulatory 455 changes which are more strongly associated with C. briggsae dominant expression in males 456 (Fig. 3). These results align well with observations of downregulation of spermatogenesis genes, 457 such as fog-1, by specific transcription factors (i.e., tra-1), and sperm-specific expression 458 depending more on upstream promoter regions than 3-UTRs in C. elegans (69,70). One 459 potential explanation involves the fixation of regulatory changes that facilitate resolution of 460 genomic inter-sexual conflict through sex-biased expression (i.e., Rice's hypothesis (71)). For 461 example, more sexual conflict is expected in outcrossing than selfing species, such as C. nigoni, 462 due to stronger sexual selection of male traits (72). To avoid traits that are detrimental to 463 females but improve male performance, genomic conflict resolution by means of sex-biased 464 expression may be attained faster through trans-regulatory changes, which are more 465 pleiotropic, downregulating male-biased genes in females. This logic aligns with the hypothesis 466 that sex-biased expression is partly driven by selection acting to resolve sexual conflict by 467 means of modifier alleles or regulators (55,71). However, the fact that trans-only regulatory 468 changes do not predominate in the control of female-biased genes in males (Fig. 4E), suggests 469 that regulatory mechanism to resolve genomic sexual conflict act in different ways for the two 470 sexes. 471 472 Several studies associate regulatory divergence in cis-acting factors with genome-wide 473 expression profiles that enhance sex-specific traits and sex-biased expression (55,73,74). Our 474 analyses are consistent with these observations in terms of the higher proportion of genes 475 under cis-only compared to trans-only regulatory divergence among genes with significant sex-476 biased expression in their respective sex (i.e., male-biased in males; Fig. 4E). When also 477 considering genes with cis-trans divergence, however, we see that a large portion of sex-biased 478 genes show significant trans-divergence ( Fig. 4E; Supplementary Fig. S6) Faster evolution of male-biased genes is the premise behind the "faster male" and "fragile 495 male" hypotheses to explain the high incidence of hybrid male sterility in XY and X0 496 heterogametic systems (31). While we find that male-biased genes collectively do not show a 497 strong signal of faster molecular evolution, the subset of male-biased genes that show 498 exceptionally high expression divergence do have faster evolving coding sequences (M4 and 499 M6; Fig. 4A-C). Additionally, these genes are implicated in spermatogenesis, based on C. 500 elegans orthologs, and show upregulated expression in sperm-producing C. briggsae 501 hermaphrodites as well as males of both species (M6; Fig. 4F). These findings accord with faster 502 molecular evolution of spermatogenesis and male germline genes of C. elegans (75-77). Their 503 rarity on the X-chromosome (Fig. 4G), however, suggests the "faster X" model does not provide 504 a compelling explanation for Haldane's rule on hybrid sterility (32). Nevertheless, our 505 transcriptome analyses support the idea that the X-chromosome plays an especially important 506 role in hybrid male sterility. This "large-X" effect arises despite just a few highly dysfunctional X-507 autosome incompatibilities between C. briggsae x C. nigoni potentially explaining hybrid male 508 sterility (52,53), in contrast to the numerous X-linked hybrid male sterility factors reported for 509 Drosophila (35,37). 510 511 Gene expression in hybrid males predominantly shows either simple dominance or no change 512 (Fig. 2C). While it is tempting to speculate that regulatory changes affecting males tend to be 513 generally more conserved as males of different species share the same reproductive role, 514 reduced sexual selection in C. briggsae males (77), genomic divergence (50,51), and clear sex 515 differences in hybrid fertility (46,47), suggest otherwise. If most transgressive expression occurs 516 in the gonad, then the small and defective gonad development of F1 males may have led to 517 their observed paucity of transgressive expression. Two non-mutually exclusive ways in which 518 hybrid male dysfunction (i.e., sterility) can arise are: 1) through misexpression and 519 misregulation of X-linked genes involved with male function, and 2) through negative epistatic 520 interactions (i.e., incompatibilities) between X-linked and autosome genes involved in male-521 specific pathways. Our results suggest that both cases are plausible. 522 523 First, the paucity of X-linked sex-biased genes in parental genotypes of Caenorhabditis species 524 suggests that any misregulation and misexpression on the X might exert little downstream 525 impact ( Fig. 1E; (40,41)). However, misexpression of X-linked genes in hybrids is relatively 526 common compared to autosomes in both sexes (Fig. 2D), with hybrid males having higher 527 relative incidence of effectively misregulated genes (trans-only, enhancing and compensatory 528 cis-trans changes) compared to female hybrids (Fig. 2H). We find that trans-acting factors often 529 contribute to misexpression in both sexes (Figs. 2H, 3). In hybrid females, trans-acting factors 530 largely drive the expression of X-linked genes with C. briggsae dominant and underdominant 531 expression, unlike autosomes (Fig. 3). In hybrid males, this effect is even more pronounced, in 532 part due to our inference that all X-linked genes with C. briggsae dominant expression arise 533 from trans-only effects. These findings are consistent with previous observations, particularly in 534 Drosophila, of trans-acting changes sex-specific causing misregulation of X-linked genes 535 (26,74,78). 536 537 Second, we find extensive expression dominance in F1 males that disproportionately matches 538 the C. briggsae expression level and that have strong cis effects (Fig. 3). Many of these genes 539 also have biased expression of C. briggsae alleles (Supplementary Fig. S7) and therefore have 540 the potential of disrupting gene networks as they may interact negatively with C. nigoni X-541 linked genes in hybrid males. Autosomal spermatogenesis genes, by contrast, tend to show C. 542 nigoni-dominant expression in F1 hybrid males (mean of 2.17-fold difference across autosomes 543 relative to genes with C. briggsae-dominant expression) (Fig. 4G), consistent with previous work 544 showing recessive effects of the C. briggsae autosomal portions of genetic incompatibilities 545 (52). In addition, this prior work also showed that sterility in C. nigoni x C. briggsae hybrid males 546 may not require many X-autosome incompatibilities (53). Despite their low abundance on the 547 X-chromosome, X-linked spermatogenesis genes are often enriched for both misexpression and 548 misregulation (Fig. 5C), potentially enhancing their role in hybrid dysfunction. Interestingly, we 549 find that at least three X-linked spermatogenesis genes that are effectively misregulated (i.e. 550 with trans-only, cis-trans enhancing and cis-trans compensatory regulatory changes) also have 551 high rates of molecular evolution (Supplementary Fig. S8) and lie near to an X-chromosome 552 segment implicated previously in hybrid male sterility (53). Our genome-wide transcriptome 553 analysis of cisand trans-regulatory divergence therefore reinforces some previous key 554 inferences about hybrid dysfunction associated with males, spermatogenesis, and the X- and polymorphism compared to central regions ( Fig. 5A; Supplementary Fig. S5; (45,80)), and 562 therefore have a greater potential to facilitate fixation of weakly beneficial cis-regulatory 563 mutations. We first predicted no effect of genomic region on differences in trans regulatory 564 divergence and confirmed this null expectation. cis-regulatory divergence, however, showed 565 only a weak elevation in chromosome arms compared to centers, with the signal being 566 somewhat stronger in females and for some autosomes (Fig. 5B,C). We further looked at 567 broader correlations for both coding and non-coding upstream sequence divergence with cis 568 regulatory divergence (Fig. 5D) and found only weak positive correlations. In line with previous 569 studies (81,82), these results indicate that rates of regulatory divergence due to cis-acting 570 elements are largely decoupled from rates of molecular evolution. 571

572
We contrasted sex-specific transcriptomic profiles between C. briggsae and C. nigoni and their 573 hybrids to understand how the evolution of cisand trans-regulatory elements can drive F1 574 hybrid dysfunction. Such evolution may arise from divergent expression changes as well as with 575 stabilizing changes that lead overall expression to remain conserved between species. The 576 sharp contrast of Caenorhabditis morphological stasis and extensive expression conservation 577 between species with extensive misexpression in F1 hybrids indicates substantial 578 developmental system drift of regulatory networks that destabilize in hybrids to enforce 579 reproductive isolation between species. Despite more extensive transgressive expression in 580 hybrid females, they are fertile but unable to produce self-sperm compared to the entirely 581 sterile hybrid males, suggesting that hybrid females may represent "demasculinized" 582 hermaphrodites through the disruption of sperm-specific regulatory networks. Despite the 583 rarity of sex-biased genes on the X-chromosome, the X is home to disproportionate 584 misexpression in both sexes, with misregulation in hybrid males largely through trans-acting 585 factors. X-autosome incompatibilities in hybrid males likely result from the propensity for C. To obtain allele-specific read counts in F1 hybrids, we applied a competitive read mapping 624 approach using a custom Python script (https://github.com/santiagosnchez/CompMap) that 625 uses the PYSAM library (https://github.com/pysam-developers/pysam). We then compared the 626 alignment score (AS) and number of mismatches (nM) to both reference genomes, retaining the 627 best single read alignments and excluding ambiguous reads (i.e., alignments with the same 628 value in both parents). We have high power to detect ASE, given ~20% neutral sequence 629 divergence between C. briggsae and C. nigoni (45) that confers an expected ~5 nucleotide 630 differences for every 100 bp of coding sequence (0.2 divergence * 0.25 fraction of synonymous 631 sites * 100 bp). To account for potential mapping bias (85) and unaccounted ambiguous reads, 632 we subjected all samples to competitive read-mapping (hybrids and pure species) and retained 633 only unambiguously mapped reads. 634

635
Ortholog identification and read abundance quantification 636 We quantified gene expression abundance for a set of 12,115 genes that we inferred to be one-637 to-one reciprocal orthologs between C. briggsae and C. nigoni. To identify orthologs, we applied 638 a phylogenetic approach using ORTHOFINDER v2. Genes were then filtered based on the amount of expression using EDGER's filterByExpr function 660 (93). Expression counts were then log2-transformed using the voom mean-variance trend 661 method to ensure consistent, normalized read counts across samples (94). Before statistically 662 assessing differential expression, we summed the allele-specific counts from F1 hybrids to yield 663 a single count of transcripts per gene. We visualized the overall expression distance between 664 samples using a non-metric multi-dimensional scaling plot, which showed all three biological 665 replicates to cluster consistently within their corresponding treatment (Fig. 1A). We inferred 666 sex-biased gene expression by comparing differential expression profiles between males and 667 females (or hermaphrodites) in each genetic group (C. briggsae, C. nigoni, F1 hybrids). We also 668 quantified differential expression between the genetic groups in a pairwise manner (C. briggsae 669 vs F1, C. briggsae vs C. nigoni, C. nigoni vs F1) for each sex separately. We then contrasted 670 expression patterns between species (C. briggsae and C. nigoni) by looking at sex differences 671 (sex-biased expression) and their interaction (expression ~ species * sex). Linear regression 672 models were applied to make statistical inferences on differential expression with the lmFit and 673 both parental species, meaning that there were significant differences in expression between 683 all groups in a manner where expression levels in F1s fall in between both species. Genes with 684 dominant allelic effects were those with expression levels in F1s matching either one of the 685 parent species (i.e. no significant differential expression), but with significant differential 686 expression between species. Finally, genes with significant differential expression from both 687 parents, but that were either significantly underexpressed (overdominant) or overexpressed 688 (underdominant) compared to both species were regarded as transgressive. Genes falling 689 outside any these specific categories were considered ambiguous. 690

691
We also measured absolute Euclidean distances in expression relative to the centroid or origin 692 in expression space of F1 hybrids relative to both parent species. For example, for every gene 693 we took the expression difference between F1s and C. briggsae and between F1s and C. nigoni 694 as an xy coordinate system. Then, we measured the Euclidean distance from that point in 695 expression space to the origin (0,0), reflecting no change in expression: 696 697 = $(∆ !"/$%& ) ' + (∆ !"/$%& ) ' 698 Where ∆ !"/$%& and ∆ !"/$() are coefficients of differential expression between F1 hybrids and 699 each parent species. This metric allowed us to visualize the magnitude of expression distance 700 from a "conserved" expression profile. 701 702 cis-and trans-regulatory divergence 703 We also used ASE in F1s to quantify the extent and type of cisand trans-regulatory differences 704 between species. Expression divergence between parent species results from both cisand 705 trans-regulatory changes, whereas significant differential expression between alleles in F1s 706 results from cis-regulatory divergence only (16). To quantify the extent of trans effects, we 707 applied a linear model to test for differences in gene expression between parent species (P) and 708 between alleles in F1 hybrids (ASE) using the following model: expression ~ species/group, 709 where "group" represented categorical variables pointing to data from P and ASE. The division 710 operator of the function "/" measures expression ratios independently for each category in 711 "group". We then used a post-hoc Wald-type test (linearHypothesis from the CAR package) to 712 test for significant differences between both coefficients (P[C. nigoni/C. briggsae] = ASE[C. 713 nigoni/C. briggsae]). P values were considered significant after a 5% FDR analysis (95). 714 715 We inferred the influence of cisand trans-regulatory divergence on genes linked to autosomes, 716 as well as to the X-chromosome in females, following the criteria in McManus et al. (16). This 717 procedure allowed us to designate genes having undergone significant regulatory divergence 718 due to cis-only, trans-only, and compensatory cis-trans effects. Genes with either significant 719 synergistic (cis + trans) or antagonistic (cis x trans) cis-trans effects were grouped together as 720 representing changes with enhancing cis-trans effects. Genes expressed with no significant 721 differences between parents, ASE, or trans effects were deemed as conserved and those that 722 did not strictly fit into any of the previous groups were considered ambiguous. 723 724 Given the hemizygous condition of the X-chromosome in males, we cannot use F1 ASE of X-725 linked genes in males to assess cis and trans regulatory divergence. However, we devised a 726 scheme to assign different types of regulatory divergence to X-linked genes given differences in 727 expression between male F1 hybrids and parent species (Supplementary Fig. S3; Wayne et al. 728 2004). Given that F1 males' X-chromosome derives solely from their maternal C. nigoni, X-linked 729 genes that differ in expression between the parental species and showing C. nigoni dominant 730 expression in hybrid males were considered as having significant cis-only effects, as differences 731 in the C. briggsae trans autosomal environment did not lead to significant deviations from C. 732 nigoni expression. Alternatively, X-linked genes in hybrid males found to be C. briggsae 733 dominant reflect significant trans-only effects. This implies two things: 1) that X-linked cis 734 elements in C. nigoni are not sufficiently different from their counterparts in C. briggsae to 735 prevent C. briggsae trans regulators from acting on them; and 2) that C. nigoni trans-regulators 736 on those pathways are potentially recessive. Supporting this assignment scheme, X-linked 737 genes with cis-only effects on regulatory divergence have significantly higher rates of molecular 738 evolution compared to genes with trans-only effects (Supplementary Fig. S2; ordinary linear 739 regression P value = 0.05 for Ka/Ks and P value < 0.05 for the proportion of conserved 5 bp 740 windows 500 bp upstream). Consequently, we inferred compensatory cis-trans effects for X-741 linked genes in males where expression was not significantly different between parent species, 742 but significantly different from F1 hybrid males. Lastly, we inferred significant enhancing cis-743 trans effects for those X-linked genes with intermediate (additive) expression in F1s, or with 744 significant differential expression between parent species coupled to significantly higher or 745 lower expression in F1s than in both parent species (Supplementary Fig. S3; (26)).