Alternative splicing shapes sexual dimorphism and erodes following the loss of sex in stick insects

7 Sexually dimorphic traits emerge from divergent selective pressures acting on males and 8 females and are facilitated by gene regulatory mechanisms. While gene expression 9 differences have been the subject of considerable focus in the evolution of sexual 10 dimorphism, other regulatory processes remain comparatively understudied. One such 11 mechanism is alternative splicing, whereby a single gene can encode multiple distinct 12 transcripts, which can lead to functional diversification of the proteome. Leveraging long-read 13 PacBio Iso -sequencing, we examine patterns of sex differences in gene expression and 14 alternative splicing under different selection regimes . Alternative splicing is widespread in 15 Timema, with the greatest isoform complexity found in the male gonad tissue. We observe 16 similar rates of sex differences in gene expression and splicing, and a core set of genes 17 involved in sexual differentiation and germline development that have conserved differential 18 splicing patterns across the Timema phylogeny. Beyond this, however, differentially spliced 19 genes undergo faster lineage -specific turnover and are subject to greater pleiotropic and 20 purifying constraints than differentially expressed genes. Following the transition from sexual 21 reproduction to asexuality, there is a systematic reduction in isoform diversity and erosion of 22 sex-specific splicing variation. Taken together, our findings suggest a role for sexual selection 23 in maintaining splicing complexity in sexual Timema species. 24 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 3

25 A central aim in evolutionary biology is to understand how a single genome can generate 26 multiple phenotypes, such as queen and worker castes in eusocial insects (Price et al., 2018), 27 morphologically and ecologically distinct larval, pupal and adult stages of holometabolous 28 insects (Simpson et al., 2011) , or different sexes, which represent one of the most extreme 29 and prevalent forms of intra -specific diversity (Mank, 2017) . Owing to their distinct 30 reproductive strategies, males and females are often subject to opposing selective pressures, 31 which can generate conflicts in the shared genome (Lande, 1980). Intralocus conflicts occur 32 when traits with a shared genetic basis have opposing fitness effects in different phenotypes, 33 as is the case for sexually antagonistic traits in males and females (Bonduriansky & 34 Chenoweth, 2009). These conflicts can potentially be resolved through regulatory evolution, 35 ultimately leading to the formation and maintenance of sexually dimorphic traits. 36 37 Differential gene expression has long been identified as a primary mechanism by which 38 selection mitigates the constraints of a shared genome (Connallon & Knowles, 2005; Mank, 39 2017) and a large body of research has investigated the role of sex-biased genes in sex-specific 40 adaptation (Darolti & Mank, 2023; Djordjevic et al., 2022; Harrison et al., 2015; Immonen et 41 al., 2017; Ingleby et al., 2014; Magnusson et al., 2011; Perry et al., 2014; Whittle & Extavour, 42 2019). Such studies have greatly advanced our understanding of sexually dimorphic traits by 43 identifying candidate loci and signaling pathways (Emlen et al., 2012; Galouzis & Prud'homme, 44 2021; Kopp et al., 2000; Toubiana et al., 2021) . However, gene expression levels represent 45 only one dimension of phenotypic diversification , and the contribution of other regulatory 46 mechanisms to facilitating the evolution of distinct male and female forms has been 47 comparatively understudied. 48 49 One such regulatory mechanism is post -transcriptional alternative splicing which can 50 generate structurally variable isoforms of a gene, for example through different combinations 51 of exon exclusion or intron retention events (Verta & Jacobs, 2022) . This process expands 52 proteome complexity (Nilsen & Graveley, 2010) and is known to regulate important biological 53 functions such as tissue development (Wang et al., 2008) , diverse physiological processes 54 (Kalsotra & Cooper, 2011) , and sex -determination (Schutt & Nothiger, 2000) . However, a 55 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 4 significant proportion of isoforms are likely non-functional, largely representing stochastic 56 mis-splicing (Benitiere et al., 2024; Wan & Larson, 2018) . Alternative splicing rates can vary 57 substantially across phenotypes (Jacobs & Elmer, 2021; Price et al., 2018), tissues (Naftaly et 58 al., 2021), developmental stages (Gibilisco et al., 2016), and species (Blekhman et al., 2010). 59 As such, the extent to which this splicing diversity reflects functional innovation and 60 phenotypic optimization as opposed to relaxed selection and noise remains an open question. 61 62 Different studies on the role of alternative splicing in sex -specific adaptation and the 63 resolution of sexual conflict demonstrate that sex differences in splicing extend far beyond 64 primary sex determination pathways and suggest that sex-specific selection can shape global 65 splicing patterns to facilitate sexual dimorphism (Blekhman et al., 2010; Chen et al., 2023; 66 Gibilisco et al., 2016; McIntyre et al., 2006; Naftaly et al., 2021; Rogers et al., 2021; Singh & 67 Agrawal, 2023). Despite these advances, we still lack a phylogenetically explicit understanding 68 of how sex -biased splicing evolves across lineages, how stable such patterns are through 69 evolutionary time, and to what extent they are actively maintained by sexual selection rather 70 than arising from neutral or stochastic processes. A powerful way to disentangle these forces 71 is to compare closely related species that differ in reproductive mode and, consequently, in 72 the presence of sexual selection and efficiency of purifying selection . The stick insect genus 73 Timema provides such a natural evolutionary experiment , whereby multiple independent 74 transitions from sexual reproduction with separate sexes to parthenogenesis have occurred 75 within the clade (Schwander et al., 2011) . Parthenogenetic Timema species are 76 morphologically and ecologically equivalent to their sexual sister species (Law & Crespi, 2002), 77 however, because they consist of only females, they experience no sexual selection or sexual 78 conflict (Parker et al., 2019b) . Moreover, many of the former sexually dimorphic traits are 79 expected to decrease in strength or disappear as they lose their adaptive value (Schwander 80 et al., 2013). Interestingly, as Timema species possess and XX:XO sex determination system, 81 rare males can be produced in asexual species through the accidental loss of an X 82 chromosome during oogenesis (Schwander et al., 2013) . Such asexual males may serve as 83 controls when investigating patterns of sex -biased transcriptomes, as they are expected to 84 evolve primarily under relaxed selection and genetic drift. 85 86 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 5 Here, we leverage multiple sexual and parthenogenetic Timema species pairs, together with 87 PacBio long -read Iso-sequencing (Iso-seq) data from somatic and reproductive tissues of 88 males and females , to provide an isoform -resolved analysis of sex -biased gene regulation. 89 Specifically, we compare sex differences in gene expression and alternative splicing across 90 tissues, species, and reproductive modes to test whether sex-specific splicing is conserved or 91 labile, whether it targets genes under distinct selective constraints from sex -biased 92 expression, and how these regulatory patterns change following the loss of sex and sexual 93 selection. By integrating long-read transcriptomics with a replicated evolutionary framework, 94 our study directly addresses the evolutionary forces shaping alternative splicing in the context 95 of sexual dimorphism and sexual conflict. 96

97 Biological samples, Iso-seq library preparation and sequencing 98 We conducted PacBio Iso-seq of RNA previously extracted for different projects (Parker et al., 99 2025). We used RNA extractions of three tissues (femur, gut, and gonad) of males and females 100 collected in the field as juveniles and reared to adulthood in the laboratory under standard 101 conditions (see Djordjevic et al., 2025 for details). Prior to dissection, all individuals were fed 102 with an artificial diet for two days, allowing their gut contents to be cleared of plant material 103 (their natural food source). Insects were then anesthetized with CO2, dissected, and tissues 104 individually flash -frozen in liquid nitrogen before being stored at −80°C. Extractions were 105 done using a TRIzol -Chloroform protocol (Djordjevic et al., 2025) . RNA quantity and quality 106 was measured using a fluorescent RNA-binding dye (QuantiFluor RNA System) and a Fragment 107 Analyzer (Agilent Technologies) . Equimolar amounts of RNA extractions from four different 108 individuals were then pooled to generate a single Iso -seq library per tissue, sex, and species 109 (for a total of 48 libraries, see Table S1). SMRTbell library preparation and Sequel II sequencing 110 using two SMRT cells pere sample were completed at the Lausanne Genomic Technologies 111 Facility (University of Lausanne). 112 113 Iso-seq read processing, isoform identification and filtering 114 SMRT Cell subread data was run through the PacBio Iso -Seq Analysis pipeline in SMRT Link 115 v.12, which generated full -length non -concatemer (FLNC) reads trimmed of primers and 116 poly(A) tails. Through the same pipeline, FLNCs were then clustered to produce high quality 117 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 6 consensus (HQ) transcripts, with 99% sequence accuracy and supported by a minimum of two 118 FLNCs. 119 120 HQ transcript sequences from all sexual and asexual samples in a pair (Table S1, Table S2) 121 were aligned to the genome assembly and annotation of the sexual species using pbmm2 122 v1.10 (https://github.com/PacificBiosciences/pbmm2), with the options --preset ISOSEQ, --123 sort, --unmapped. Redundant transcripts were collapsed into unique isoforms using isoseq3 124 collapse v3.8.2. To avoid a high rate of reference transcript redundancy, we set --max-5p-diff 125 1000000 and --max-3p-diff 1000000 and collapsed isoforms that match ed a reference 126 transcript at all splice junctions but had an alternative 5’ or 3’ end, which may reflect partial 127 cDNA synthesis or degradation (Tardaguila et al., 2018). 128 129 Isoform characterization and filtering was performed using SQANTI3 v5.1.1 (Pardo-Palacios 130 et al., 2024). For each species, we used previously generated male and female RNA-seq short-131 read data from the same three tissues (femur, gut, gonad) (Parker et al. 2025) to provide 132 information on isoform expression and splice junction coverage. Specifically, we aligned 133 short-reads to the non-redundant isoform sequences with kallisto v0.48 (Bray et al., 2016) to 134 obtain expression matrices and separately aligned short-reads to the reference genome and 135 non-redundant isoform annotation with STAR v2.7.11a (Dobin et al., 2013) to obtain splice 136 junction coordinates and coverage. Based on this information, we used SQANTI3 to remove 137 isoforms classified as intra -priming products, RT-switching artifacts, and isoforms with non-138 canonical junctions that have short -read coverage < 3. We further excluded mono-exonic 139 transcripts and collapsed isoforms with the same open reading frame, only preserving the 140 longest transcript for subsequent analyses. 141 142 Differential gene expression and differential transcript usage analysis 143 For each species and tissue, we first mapped the paired-end short-read RNA sequencing data 144 to the final filtered set of isoforms using Salmon v 1.10.2 (Patro et al., 2017) with 100 145 bootstraps and correcting for GC and sequence-specific biases with --seqBias and --gcBias. We 146 then extracted gene-level and transcript-level abundance matrices using the tximport v1.34 147 package (Soneson et al., 2015) in R v4.4.3 ( R Core Team ). Finally, for each tissue of each 148 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 7 species we removed lowly expressed genes and isoforms by imposing a minimum 2 TPM and, 149 respectively, 5 TPM expression threshold in at least half of the individuals of one sex. 150 151 We identified differentially expressed (DE) genes with the DESeq2 v1.46 package (Love et al., 152 2014) in R, with an adjusted false discovery rate < 0.05 and an absolute log2 fold change ≥ 1. 153 We also identified differentially spliced (DS) genes using a differential transcript usage that 154 detects changes in relative isoform abundance, irrespective of overall gene expression. For 155 this, we first transformed transcript abundance estimate s from Salmon to normalize for 156 transcript length and library size using the dtuScaledTPM method in tximport. W e then 157 modeled relative isoform proportions using the Dirichlet -multinomial model from the 158 DRIMSeq v1.34 R package (Nowicka & Robinson, 2016) with parameters 159 min_samps_feature_expr=2, min_samps_feature_prop=2, min_samps_gene_expr=2, 160 min_feature_expr=10, min_feature_prop=0.1, min_gene_expr=10. To reduce false positives 161 and control the overall false discovery rate (Love et al., 2018) , the resulting gene - and 162 transcript-level p values were subjected to stage -wise adjustment with the stageR v1.28 163 package (Van den Berge et al., 2017). A first screening stage was performed at the gene level, 164 with genes failing the adjusted p value threshold of 0.05 excluded from further analysis. 165 Secondly, for genes passing the first stage, a transcript -level significance testing was 166 performed. 167 168 Gene orthology inference 169 For each species, we used the nucleotide sequences for the filtered set of isoforms and 170 selected the longest isoform per gene. Orthologs between species were inferred using 171 OrthoFinder v2.5.5 (Emms & Kelly, 2019) with default parameters. 172 173 Estimating gene expression tissue specificity 174 Tissue specificity of gene expression is best estimated across many different tissues. In 175 addition to the short -read RNA -sequencing data for the three focal tissues with long -read 176 sequencing data (femur, gonad, and gut ), we therefore included data from antenna, brain, 177 defense glands, and fat body, and estimated gene expression as detailed above. For each 178 gene, we calculated a tissue specificity index (𝜏) (Yanai et al., 2005) based on the formula : 179 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 8 𝜏 = ∑ 1− 𝑒𝑥𝑝𝑟𝑖 𝑒𝑥𝑝𝑟𝑚𝑎𝑥 𝑁 𝑖=1 𝑁−1 , where N is the number of tissues, expri is the expression of the gene in 180 tissue i, and exprmax is the maximum expression of that gene across all tissues. 𝜏 values range 181 from 0 to 1, where higher values indicate a higher tissue -specific expression, while lower 182 values indicate a more even expression across tissues. 183 184 Estimating rates of coding-sequence evolution 185 For all single-copy orthologs identified through the gene orthology inference, we extracted 186 nucleotide coding-sequences as inferred by SQANTI3. We aligned sequences using macse 187 v2.07 (Ranwez et al., 2018) and removed alignments shorter than 150bp after gap removal. 188 We used the branch model test (model = 1, nssites = 0) in the CODEML package in PAML v4.8 189 (Yang, 2007) to estimate rates of coding-sequence evolution for each branch of the phylogeny 190 “((TpsTdi,TcmTsi),TceTms),(TpaTge,Tbi));“. To avoid divergence estimates being affected by 191 mutational saturation and to exclude cases of gene misalignment , we removed genes with 192 the rate of synonymous substitutions ( dS) > 2 as per previous studies (Axelsson et al., 2008). 193 For each category of genes ( DE and DS genes), we extracted from the PAML results the 194 number of synonymous ( DS) and non-synonymous (DN) substitutions, as well as the number 195 of synonymous (S) and non-synonymous (N) sites. For each category, we then calculated the 196 average rate of synonymous substitutions ( dS) and the average rate of non -synonymous 197 substitutions ( dN) as dS = DS / S and, respectively, dN = DN / N. Bootstrapping with 1 ,000 198 replicates was used to calculate the 95% confidence intervals for divergence estimates in each 199 group. Lastly, permutation tests with 1 ,000 replicates were used to test for significant 200 differences in dN, dS, and dN/dS between categories. 201

and Discussion 202 Alternative splicing is prevalent in Timema stick insects 203 We discovered very high levels of splicing variation in Timema stick insects, recovering a range 204 between 18,686 and 30,212 isoforms across the five species pairs (Table S2) . Using a 205 combination of PacBio Iso-seq and Illumina paired -end RNA-sequencing we found that 40-206 53% of curated Timema transcriptomes have two or more isoforms ( Fig. S1; Table S2). Our 207 estimates of splicing variation may yet be conservative as we imposed stringent filters to limit 208 false positive isoforms. The proportion of alternatively spliced genes that we observe across 209 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 9 Timema species is most directly comparable to those from recent long -read invertebrate 210 studies, such as those in planthoppers (58%) (Tong et al., 2022) and honey bees (32–68%) (Hu 211 et al., 2024; Zheng et al., 2023), which similarly highlight extensive transcriptomic complexity. 212 In contrast, earlier studies using less sensitive short -read sequencing approaches reported 213 generally lower rates, including in silkworms (27%) (Li et al., 2012) , pea aphids (34%) 214 (Grantham & Brisson, 2018) , Drosophila (37%) (Gibilisco et al., 2016) , and flatworms (42%) 215 (Wang et al., 2015). This distinction highlights that reported fractions of alternatively spliced 216 genes are often as much a reflection of sequencing sensitivity as they are of underlying 217 biology. 218 219 Transcriptional complexity varies across biological levels 220 Across all nine analyzed species, we observed similar patterns of spliceosome complexity 221 variation among tissues and between sexes. At the tissue level, we found that the gonad has 222 a significantly higher per gene isoform abundance than both the femur (Wilcoxon rank -sum 223 test, W = 9, p = 0.007) and the gut (Wilcoxon rank -sum test, W = 0, p < 0.001) (Fig. 1A, Fig. 224 S2). Of the genes expressed in all three tissues, 50% have more than one isoform in the gonad, 225 in contrast to 34% in the femur and the gut. A high splicing complexity in the gonad emerges 226 as a systematic pattern across multiple clades (Gibilisco et al., 2016; Mazin et al., 2021; Naftaly 227 et al., 2021), with brain and liver being the few somatic tissues with similarly high levels of 228 splicing (Gibilisco et al., 2016; Grosso et al., 2008; Mazin et al., 2021) . Compared to somatic 229 tissues, gonads express the largest number of splicing regulators and exhibit high differential 230 expression of splicing factor genes (Grosso et al., 2008) . Together with the extreme 231 heterogeneity in germ and somatic cells present in the gonad, this might underly the observed 232 isoform diversity in the reproductive tissue. 233 234 Gonad spliceosome complexity further differs between the sexes. Males have, on average, a 235 larger number of expressed genes (one -sample t -test, t(6) = 11.5, p < 0.001, 95% CI 236 [1.11,1.17]; Fig. S3), and a greater per gene isoform diversity (one-sample t-test, t(6) = 7.6, p 237 < 0.001, 95% CI [1.12,1.23]; Fig. 1B, Fig. S4) compared to females in the reproducti ve tissue. 238 By contrast, there is no difference in transcriptional complexity between the sexes for either 239 of the somatic tissues. 240 241 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 10 This observed male-bias in both isoform and gene expression diversity contrasts with some 242 previous reports in Drosophila species which suggested that males and females adopt 243 different strategies for diversifying their transcriptome, wit h testis expressing more genes 244 while a higher fraction of alternatively spliced genes being found in ovary (Gibilisco et al., 245 2016). Recent single-cell atlases have consistently identified the testis as one of the tissues 246 with the most distinct cell types (Han et al., 2018; Han et al., 2020; Li et al., 2022) , and the 247 increase in expressed genes and isoform diversity in males that we observe here could simply 248 reflect a greater cellular diversity in testis compared to ovaries. Increased male transcriptional 249 activity could also be due to a more open chromatin state in spermatocytes and spermatids 250 during meiosis and spermiogenesis (Soumillon et al., 2013) . Importantly, this process can 251 generate a substantial amount of stochastic splicing noise, as the permissive chromatin 252 environment may lead to a saturation of the splicing machinery and a reduced precision in 253 splice-site recognition (Mazin et al., 2021) . Our analysis shows, however, that the higher 254 isoform diversity in Timema male gonad does not result from a disproportionate increase in 255 non-coding relative to coding isoforms (Fig. S5). Future coupling of alternative splicing data 256 with single -cell transcriptomics or chromatin accessibility analysis could help disentangle 257 between these factors. Alternatively, sperm competition has been implicated in driving rapid 258 testis transcriptome evolution (Trost et al., 2023) and, to some extent, this might also 259 accelerate splicing rates in males (Rogers et al., 2021) , although comparative tests at the 260 isoform diversity level are lacking. 261 262 Although alternative splicing is widespread across eukaryotes, the extent to which such 263 events are functional and adaptive, as opposed to splicing noise, remains unclear (Wan & 264 Larson, 2018). Isoform diversity tends to scale with organismal complexity (Chen et al., 2014). 265 However, according to the “drift barrier” theory, deleterious mutations due to small effective 266 population size (Ne) can also expand isoform diversity by increasing splicing errors (Benitiere 267 et al., 2024) . Indeed, a n a nalysis across dozens of species revealed an inverse correlation 268 between rates of alternative splicing and proxies for Ne, a pattern preponderantly due to low 269 abundance and non-functional splicing variants (Benitiere et al., 2024). 270 271 Our findings in Timema argue against a purely drift-driven model of splicing complexity. The 272 autosomal effective population size varies more than tenfold across sexual Timema species, 273 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 11 ranging from ~162,000 in T. poppense to ~1,900,000 in T. podura. (Parker et al., 2022), yet we 274 find no significant correlation between Ne estimates and the rate of alternative splicing in 275 these species (Fig. S6). The lack of recombination in asexual Timema species is expected to 276 significantly reduce the efficacy of selection, leading to an accelerated accumulation of 277 deleterious mutations (Bast et al., 2018). This relaxed constraint can decrease the efficiency 278 of the spliceosome and increase mis-splicing events (Benitiere et al., 2024) in asexual 279 compared to sexual species. Contrary to this expectation, we find a significantly greater 280 isoform diversity in sexual Timema species compared to parthenogenetic species in both 281 femur (one-sample t-test, t(6) = 5.4, p = 0.002, 95% CI [1.11,1.31]) and gonad (one-sample t-282 test, t(5) = 2.7, p = 0.04, 95% CI [1.00,1.16]), but not in gut ( one-sample t-test, t(6) = 2.2, p = 283 0.07, 95% CI [0.99,1.22] ) (Fig. 1C, Fig. S7). Our combined results thus indicate that splicing 284 variation in sexual Timema is not primarily driven by genetic drift and stochastic events, and 285 may in fact reflect functional complexity. 286 287 Conservation and turnover of sex differences in splicing 288 Alternative splicing differences between the sexes have commonly been regarded as 289 secondary to sex-biased gene expression. Our analyses of global patterns of differential gene 290 expression and differential transcript usage between males and females across the five 291 Timema sexual species suggest instead a comparable prevalence of these two regulatory 292 mechanisms. Consistently across species, we find the gonad tissue to be highly 293 transcriptionally dimorphic at both the gene expression level and the splicing level, while 294 somatic tissues show only minimal differentiation (Fig. 2). Although we identify twice as many 295 differentially expressed (DE) genes as differentially spliced (DS) genes in the gonad, these 296 rates converge when accounting for the total number of testable genes in each category, with 297 32% sex-biased genes and 34% genes with differential isoform abundance. In stark contrast, 298 dimorphism in the femur and gut was minor, with only 1 -2% DE genes and 4-7% DS genes. 299 Our finding that gonad exhibits profound sexual dimorphism a cross both regulatory levels is 300 consistent with previous research (Blekhman et al., 2010; Gibilisco et al., 2016; Naftaly et al., 301 2021; Rogers et al., 2021; Singh & Agrawal, 2023) . However, while sex-biased splicing has 302 been viewed as a more minor component of transcriptional dimorphism compared to 303 differential gene expression (Gibilisco et al., 2016; Rogers et al., 2021), the fact that we find a 304 comparable prevalence between these two mechanisms suggests that alternative splicing can 305 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 12 also be a fundamental driver of sex -specific differentiation. Our detection of such a high 306 proportion of genes with differential transcript usage may stem from the application of long-307 read sequencing, which captures full -length isoforms that are often undetected with short -308 read RNA -sequencing pipelines. The functional contribution of these two regulatory 309 mechanisms to sex-specific differentiation warrants further validation. 310 311 Sex-biased gene expression can exhibit rapid evolutionary turnover across species and 312 lineages (Grath & Parsch, 2016; Harrison et al., 2015) . In contrast, the conservation versus 313 turnover of sex -biased alternative splicing remain comparatively unexplored. Here, we 314 present the first explicit phylogenetic test of these dynamics. The proportion of DS genes are 315 remarkably consistent across the five Timema sexual species (Fig. 2), which prompted us to 316 investigate the extent to which these patterns of gonadal differential regulation are 317 conserved across the phylogeny or represent independent lineage -specific events. We 318 identified 6,311 one-to-one orthologous genes across the five species. Based on this 319 orthologous gene set, we find a significantly high overlap across all species for both DE and 320 DS genes (SuperExactTest padj < 0.001; Fig. 3A, B), suggesting a conserved ancestral core set 321 of genes that may be involved in sexual differentiation. Indeed, among the 54 genes with a 322 conserved differential splicing pattern in all five species, we recover several major splicing 323 regulators and genes involved in germline development (Table 1). These results indicate that 324 sex-biased splicing in sexual Timema species is a functionally driven process rather than a by-325 product of stochastic noise. 326 327 Although we observe a slight trend where genes with a conserved differential splicing status 328 across four or five species generally possess a higher number of isoforms compared to genes 329 with species-specific differential splicing (Fig. S8A), we still recover a significant overlap in DS 330 genes across species when limiting the analysis to genes with a maximum of three isoforms 331 (Fig. S8B; SuperExactTest padj < 0.001, fold enrichment = 10). This confirms that the observed 332 conservation of DS genes is not simply due to high transcript complexity in specific loci, but 333 instead reflects a stable regulatory program. 334 335 Sex-biased genes, and in particular male -biased genes, frequently exhibit accelerated 336 evolutionary rates at both the coding-sequence and expression levels (Grath & Parsch, 2016). 337 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 13 This rapid divergence has been attributed to intense sexual selection acting in males (Proschel 338 et al., 2006; Sawyer et al., 2007), however, relaxed pleiotropic constraints may also facilitate 339 these transitions (Dapper & Wade, 2020; Djordjevic et al., 2022; Harrison et al., 2015) . 340 Consistent with these patterns, we observe that while a significant core of both male-biased 341 and female-biased genes is shared across Timema species (SuperExactTest padj < 0.001; Fig. 342 S9), a substantial fraction of the sex-biased transcriptome is still species-specific or subject to 343 rapid change. While high rates of expression turnover have been documented in diverse taxa, 344 including Drosophila (Whittle & Extavour, 2019) , Galliform birds (Harrison et al., 2015) , 345 mosquitoes (Papa et al., 2017), and fish (Lichilin et al., 2021), the evolutionary stability of sex-346 biased splicing has remain ed largely unexplored. Our results reveal a striking disparity in 347 conservation between these two regulatory mechanisms, where the enrichment for 348 conserved sex-biased gene expression (54-fold enrichment) is ten times greater than that for 349 differential splicing (5 -fold enrichment). This confirms that alternative splicing undergoes 350 faster lineage-specific turnover (Barbosa-Morais et al., 2012; Gibilisco et al., 2016; Merkin et 351 al., 2012; Rogers et al., 2021). 352 353 Gene splicing may be subject to greater evolutionary flexibility than the more constrained 354 gene-level expression landscape, facilitating rapid lineage-specific adaptation and proteomic 355 diversification. Alternatively, the elevated turnover of DS relative to DE genes may be driven 356 by differences in the mutational target size of these two regulatory mechanisms. While 357 differential gene expression is to a large degree specified by promoters and enhancers, the 358 fidelity of alternative splicing depends on a complex network of cis -regulatory signals, 359 including splice sites , enhancers and silencer s that are distributed across exons and introns 360 throughout the gene body (Wang & Burge, 2008; Yeo et al., 2007). 361 362 Differentially expressed and spliced genes are under distinct constraints 363 An emerging pattern from genome -wide analyses suggests that transcriptional and splicing 364 regulatory mechanisms are largely non -redundant as they target different gene sets and 365 affect distinct biological processes (Grantham & Brisson, 2018; Jacobs & Elmer, 2021; Rogers 366 et al., 2021; Singh & Agrawal, 2023) . In the context of sexual dimorphism, it has been 367 proposed that DE and DS genes may resolve sexual conflict through distinct molecular 368 strategies (Rogers et al., 2021; Singh & Agrawal, 2023) . Specifically, while sex -biased 369 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 14 expression shifts the overall quantity of a protein towards a sex-specific optimum, sex-biased 370 splicing enables the production of sexually dimorphic isoforms from a shared genomic locus. 371 Therefore, alternative splicing may offer a complementary way towards the resolution of 372 sexual conflict and the development of sexual dimorphism by acting on genes that are under 373 stronger constraints (Rogers et al., 2021). 374 375 Throughout the Timema phylogeny, although we find that the overlap between DE and DS 376 genes in the gonad tissue is not significantly different than expected by chance 377 (Hypergeometric tests, p > 0.05), on average 68% of DS genes have an unbiased expression 378 between males and females (Fig. 4A). This is consistent with previous findings that, to a large 379 extent, these two regulatory mechanism s impact different sets of genes. We also employed 380 two measures for testing whether DE and D S genes are subject to different selective 381 constraints. First, we used expression data across a range of tissues to obtain tissue specificity 382 estimates. The evolution of sex -biased gene expression can be limited by pleiotropic 383 constraints, as genes expressed across multiple tissues or developmental stages are often 384 unbiased in expression and subject to stabilizing selection to maintain their functions in both 385 males and females (Mank et al., 2008; Meisel, 2011) . Alternative splicing may circumvent 386 these pleiotropic constraints on gene expression level by generating distinct male and female 387 isoforms (Rogers et al., 2021) . In line with this, w e observe that DS genes are significantly 388 more broadly expressed than DE genes (Fig. 4B; Wilcoxon rank -sum tests, p < 0.05 in all 389 comparisons). Secondly, we calculated rates of coding-sequence evolution and found that DS 390 genes have lower rates of evolution than DE genes, attributable to a significantly lower rate 391 of non -synonymous substitutions (Fig. 4C,D). Taken together, our results suggest that 392 differential splicing is targeting genes that are subject to stronger purifying selection and 393 functional constraints. 394 395 Erosion of sex-biased splicing patterns in asexual species 396 The comparison between sexual and asexual species provides a natural experiment to 397 disentangle selection, drift, and developmental constraint. The transition from sexual 398 reproduction to asexuality significantly alters the selective landscape, which can lead to 399 genome-wide transcriptomic changes (Bast et al., 2018; Liu et al., 2014; Parker et al., 2019a). 400 Although reduced selective efficacy in asexual lineages might predict random drift of gene 401 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 15 expression away from its optima (Huylmans et al., 2021) , empirical evidence reveals that 402 more complex selective dynamics are at work. Asexual species are relieved from sexual 403 selection and sexual confl ict, and t his resolution should in theory allow asexual species to 404 invest more towards female -specific traits and functions, potentially manifesting as large -405 scale transcriptomic shifts, such as the upregulation of female -biased genes or 406 downregulation of formerly male -biased ones (Parker et al., 2019b) . Contrary to this 407 expectation, studies including those in Timema have found repeated patterns of 408 desexualization of gene expression in asexual females, characterized by a lowered expression 409 for female-biased genes, consistent with a change in female trait optima in asexual lineages 410 and the decay of sexual traits (Huylmans et al., 2021; Parker et al., 2019b) . Comparative 411 studies of sexual and asexual lineages can thus provide a powerful test for determining 412 whether patterns of sex differences in splicing are maintained by sexual selection. Moreover, 413 Timema asexual males, which can be produced sporadically via loss of an X chromosome in 414 XO systems and which are subject to relaxed selection pressures, offer a unique opportunity 415 to further disentangle the effect of selection from drift on splicing changes. 416 417 For the two sexual – asexual species pairs (T. poppense – T. dougla si; T. podura – T. 418 genevievae) for which we have gonadal isoform information for both sexes, we find that, on 419 average, 47% of DS genes between sexual males and females are also differentially spliced in 420 the asexual species. T hese sex differences in splicing likely represent broad developmental 421 constraints that persist even after the transition from sexual reproduction to 422 parthenogenesis. However, a large fraction of DS genes found in the sexual species do not 423 preserve their splicing pattern in the asexual lineages. To understand why that is the case, we 424 investigated the direction of change in isoform abundance between sexual and asexual 425 individuals based on three potential scenarios (Fig. 5A). Asexual females may generally retain 426 female-specific splicing patterns under developmental constraint. Consistent with the trends 427 found at the gene expression level (Parker et al., 2019b) , there may also be a shift in 428 differential splicing patterns towards a loss of female -abundant isoforms and an increase in 429 male-abundant transcripts. Alternatively, due to relaxed purifying selection and drift, asexual 430 lineages may exhibit an increased variance or noise in isoform abundance. 431 432 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 16 Consistently across comparisons, there is a positive correlation in isoform abundances 433 between sexual and asexual individuals (Fig. 5B, C ), in line with the notion that regulatory 434 conservation or broad developmental constraints preserve ancestral splicing patterns even 435 after the transition to asexuality. However, the magnitude of this correlation is attenuated in 436 males compared to females , highlighting the impact of increased genetic drift in asexual 437 males. Furthermore, we repeatedly identify a pattern of desexualization of isoform 438 abundance, wherein both asexual females and males exhibit a decreased abundance for some 439 sexual female-biased isoforms and, respectively, male-biased isoforms (Fig. 5B, C lower right 440 quadrants of each panel ; Wilcoxon signed -rank tests, p > 0.001 ). The loss of sex -biased 441 isoforms in asexual individuals may reflect both relaxed sexual selection and increased drift . 442 However, the absence of increased splicing noise in asexuals and no expression changes for 443 opposite-sex-biased isoforms run counter to predictions from simple drift-barrier models and 444 suggest that sexual selection maintains splicing complexity in sexual Timema species. 445 Concluding Remarks 446 Alternative splicing is a pervasive and structured component of gene regulation in stick 447 insects. It contributes substantially to tissue specificity, sexual dimorphism, and evolutionary 448 diversification. By leveraging long-read transcriptomics across multiple species, tissues, sexes, 449 and reproductive modes, we provide a comprehensive isoform -resolved view of splicing 450 regulation in this system. Our results reveal several key insights. First, stick insects exhibit high 451 levels of isoform diversity, particularly in gonadal tissue, matching or exceeding those 452 reported in other invertebrates. This complexity is not explained by relaxed selection or 453 stochastic splicing noise, but instead is highest in sexual species and in biologically relevant 454 contexts, suggesting a functional role for alternative splicing in reproductive biology. Second, 455 when considering the number of testable genes, sex-specific alternative splicing is as 456 prevalent as sex-biased gene expression in gonads, indicating that splicing is not a secondary 457 or minor contributor to sexual dimorphism, but a core regulatory mechanism. Despite rapid 458 lineage-specific turnover, we identify a conserved core of differentially spliced genes, pointing 459 to deep evolutionary constraints on key components of sexual differentiation. Third, 460 differential splicing and differential expression act on largely distinct gene sets and are subject 461 to different selective regimes , whereby splicing preferentially targets more broadly 462 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 17 expressed, slowly evolving genes under strong purifying selection, consistent with a role in 463 resolving sexual conflict without altering overall gene dosage. Finally, by comparing sexual 464 and asexual lineages, we show that sex-specific splicing patterns are eroded following the loss 465 of sexual reproduction. Contrary to expectations under a simple drift -based model, asexual 466 Timema do not exhibit increased splicing noise , but instead show reduced isoform diversity 467 and a systematic loss of some of the sex-specialized isoforms. This pattern suggests that 468 sexual selection and sex -specific optima actively maintain splicing complexity in sexual 469 species. Taken together, our study establishes alternative splicing as a central, evolutionarily 470 dynamic, yet a selectively constrained layer of gene regulation shaping sexual dimorphism 471 and reproductive evolution. It also highlights the power of long -read transcriptomics to 472 uncover regulatory variation, and positions stick insects as a valuable model for dissecting the 473 evolutionary forces governing transcriptome complexity. 474

475 This work was supported by the European Molecular Biology Organization (EMBO) (grant 476 number 167-2022 to I.D.). We thank the Lausanne Genomic Technologies Facilities for PacBio 477 Iso-seq library preparation and sequencing. We would also like to thank members of the 478 Schwander group and Linley M. Sherin for helpful comments and discussions on this work. 479 Author Contributions 480 I.D. and T.S. designed the research. M.L. performed molecular work. I.D., V.M. analy zed the 481 data with input from T.S.. I.D. and T.S. wrote the paper with input from all authors. 482 Data availability 483 Iso-seq reads are available under NCBI BioProject PRJNA1149183. Scripts for data processing 484 and analysis are available at https://github.com/idarolti/Darolti_etal_2026_AltSpli. 485 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 18

486 Axelsson, E., Hultin-Rosenberg, L., Brandstrom, M., Zwahlen, M., Clayton, D. F., & Ellegren, H. 487 (2008). Natural selection in avian protein -coding genes expressed in brain. Mol Ecol, 488 17(12), 3008-3017. https://doi.org/10.1111/j.1365-294X.2008.03795.x 489 Barbosa-Morais, N. L., Irimia, M., Pan, Q., Xiong, H. Y., Gueroussov, S., Lee, L. J., Slobodeniuc, 490 V., Kutter, C., Watt, S., Colak, R., Kim, T., Misquitta -Ali, C. M., Wilson, M. D., Kim, P. 491 M., Odom, D. T., Frey, B. J., & Blencowe, B. J. (2012). The evolutionary landscape of 492 alternative splicing in vertebrate species. Science, 338(6114), 1587 -1593. 493 https://doi.org/10.1126/science.1230612 494 Bargiela, A., Llamusi, B., Cerro -Herreros, E., & Artero, R. (2014). Two enhancers control 495 transcription of Drosophila muscleblind in the embryonic somatic musculature and in 496 the central nervous system. PLoS One , 9(3), e93125. 497 https://doi.org/10.1371/journal.pone.0093125 498 Bast, J., Parker, D. J., Dumas, Z., Jalvingh, K. M., Tran Van, P., Jaron, K. S., Figuet, E., Brandt, 499 A., Galtier, N., & Schwander, T. (2018). Consequences of Asexuality in Natural 500 Populations: Insights from Stick Insects. Mol Biol Evol , 35(7), 1668 -1677. 501 https://doi.org/10.1093/molbev/msy058 502 Benitiere, F., Necsulea, A., & Duret, L. (2024). Random genetic drift sets an upper limit on 503 mRNA splicing accuracy in metazoans. Elife, 13. https://doi.org/10.7554/eLife.93629 504 Blekhman, R., Marioni, J. C., Zumbo, P., Stephens, M., & Gilad, Y. (2010). Sex -specific and 505 lineage-specific alternative splicing in primates. Genome Res , 20(2), 180 -189. 506 https://doi.org/10.1101/gr.099226.109 507 Bonduriansky, R., & Chenoweth, S. F. (2009). Intralocus sexual conflict. Trends Ecol Evol, 24(5), 508 280-288. https://doi.org/10.1016/j.tree.2008.12.005 509 Bray, N. L., Pimentel, H., Melsted, P., & Pachter, L. (2016). Near-optimal probabilistic RNA-seq 510 quantification. Nat Biotechnol, 34(5), 525-527. https://doi.org/10.1038/nbt.3519 511 Chen, L., Bush, S. J., Tovar -Corona, J. M., Castillo -Morales, A., & Urrutia, A. O. (2014). 512 Correcting for differential transcript coverage reveals a strong relationship between 513 alternative splicing and organism complexity. Mol Biol Evol , 31(6), 1402 -1413. 514 https://doi.org/10.1093/molbev/msu083 515 Chen, W., Zhou, W., Li, Q., & Mao, X. (2023). Sex differences in gene expression and 516 alternative splicing in the Chinese horseshoe bat. PeerJ, 11, e15231. 517 https://doi.org/10.7717/peerj.15231 518 Connallon, T., & Knowles, L. L. (2005). Intergenomic conflict revealed by patterns of sex -519 biased gene expression. Trends Genet , 21(9), 495 -499. 520 https://doi.org/10.1016/j.tig.2005.07.006 521 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 19 Dapper, A. L., & Wade, M. J. (2020). Relaxed Selection and the Rapid Evolution of 522 Reproductive Genes. Trends Genet , 36(9), 640 -649. 523 https://doi.org/10.1016/j.tig.2020.06.014 524 Darolti, I., & Mank, J. E. (2023). Sex-biased gene expression at single-cell resolution: cause and 525 consequence of sexual dimorphism. Evol Lett , 7(3), 148 -156. 526 https://doi.org/10.1093/evlett/qrad013 527 Deduyer, I., Montembault, E., Claverie, M. C., Touly, N., McCusker, D., & Royou, A. (2025). 528 Dual modes of contractility tailor cytokinesis to the requirements of diverse tissues. 529 iScience, 28(11), 113694. https://doi.org/10.1016/j.isci.2025.113694 530 Djordjevic, J., Dumas, Z., Robinson -Rechavi, M., Schwander, T., & Parker, D. J. (2022). 531 Dynamics of sex -biased gene expression during development in the stick insect 532 Timema californicum. Heredity (Edinb) , 129(2), 113 -122. 533 https://doi.org/10.1038/s41437-022-00536-y 534 Djordjevic, J., Tran Van, P., Toubiana, W., Labedan, M., Dumas, Z., Aury, J. M., Cruaud, C., 535 Istace, B., Labadie, K., Noel, B., Parker, D. J., & Schwander, T. (2025). Dynamics of X 536 chromosome hyper -expression and inactivation in male tissues during stick insect 537 development. PLoS Genet , 21(3), e1011615. 538 https://doi.org/10.1371/journal.pgen.1011615 539 Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Batut, P., Chaisson, M., 540 & Gingeras, T. R. (2013). STAR: ultrafast universal RNA -seq aligner. Bioinformatics, 541 29(1), 15-21. https://doi.org/10.1093/bioinformatics/bts635 542 Emlen, D. J., Warren, I. A., Johns, A., Dworkin, I., & Lavine, L. C. (2012). A mechanism of 543 extreme growth and reliable signaling in sexually selected ornaments and weapons. 544 Science, 337(6096), 860-864. https://doi.org/10.1126/science.1224286 545 Emms, D. M., & Kelly, S. (2019). OrthoFinder: phylogenetic orthology inference for 546 comparative genomics. Genome Biol , 20(1), 238. https://doi.org/10.1186/s13059-547 019-1832-y 548 Fear, J. M., Arbeitman, M. N., Salomon, M. P., Dalton, J. E., Tower, J., Nuzhdin, S. V., & 549 McIntyre, L. M. (2015). The Wright stuff: reimagining path analysis reveals novel 550 components of the sex determination hierarchy in Drosophila melanogaster. BMC Syst 551 Biol, 9, 53. https://doi.org/10.1186/s12918-015-0200-0 552 Galouzis, C. C., & Prud'homme, B. (2021). Transvection regulates the sex-biased expression of 553 a fly X -linked gene. Science, 371(6527), 396 -400. 554 https://doi.org/10.1126/science.abc2745 555 Gibilisco, L., Zhou, Q., Mahajan, S., & Bachtrog, D. (2016). Alternative Splicing within and 556 between Drosophila Species, Sexes, Tissues, and Developmental Stages. PLoS Genet, 557 12(12), e1006464. https://doi.org/10.1371/journal.pgen.1006464 558 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 20 Grantham, M. E., & Brisson, J. A. (2018). Extensive Differential Splicing Underlies 559 Phenotypically Plastic Aphid Morphs. Mol Biol Evol , 35(8), 1934 -1946. 560 https://doi.org/10.1093/molbev/msy095 561 Grath, S., & Parsch, J. (2016). Sex -Biased Gene Expression. Annu Rev Genet , 50, 29 -44. 562 https://doi.org/10.1146/annurev-genet-120215-035429 563 Grosso, A. R., Gomes, A. Q., Barbosa -Morais, N. L., Caldeira, S., Thorne, N. P., Grech, G., von 564 Lindern, M., & Carmo -Fonseca, M. (2008). Tissue -specific splicing factor gene 565 expression signatures. Nucleic Acids Res , 36(15), 4823 -4832. 566 https://doi.org/10.1093/nar/gkn463 567 Han, X., Wang, R., Zhou, Y., Fei, L., Sun, H., Lai, S., Saadatpour, A., Zhou, Z., Chen, H., Ye, F., 568 Huang, D., Xu, Y., Huang, W., Jiang, M., Jiang, X., Mao, J., Chen, Y., Lu, C., Xie, J.,…Guo, 569 G. (2018). Mapping the Mouse Cell Atlas by Microwell -Seq. Cell, 173(5), 1307. 570 https://doi.org/10.1016/j.cell.2018.05.012 571 Han, X., Zhou, Z., Fei, L., Sun, H., Wang, R., Chen, Y., Chen, H., Wang, J., Tang, H., Ge, W., Zhou, 572 Y., Ye, F., Jiang, M., Wu, J., Xiao, Y., Jia, X., Zhang, T., Ma, X., Zhang, Q.,…Guo, G. (2020). 573 Construction of a human cell landscape at single -cell level. Nature, 581(7808), 303-574 309. https://doi.org/10.1038/s41586-020-2157-4 575 Harrison, P. W., Wright, A. E., Zimmer, F., Dean, R., Montgomery, S. H., Pointer, M. A., & Mank, 576 J. E. (2015). Sexual selection drives evolution and rapid turnover of male gene 577 expression. Proc Natl Acad Sci U S A , 112(14), 4393 -4398. 578 https://doi.org/10.1073/pnas.1501339112 579 Hu, X. F., Jin, M. J., Gong, Z. X., Lin, Z. L., Zhang, L. Z., Zeng, Z. J., & Wang, Z. L. (2024). Full-580 Length Transcriptome Profile of Apis cerana Revealed by Nanopore Sequencing. Int J 581 Mol Sci, 25(19). https://doi.org/10.3390/ijms251910833 582 Huylmans, A. K., Macon, A., Hontoria, F., & Vicoso, B. (2021). Transitions to asexuality and 583 evolution of gene expression in Artemia brine shrimp. Proc Biol Sci , 288(1959), 584 20211720. https://doi.org/10.1098/rspb.2021.1720 585 Immonen, E., Sayadi, A., Bayram, H., & Arnqvist, G. (2017). Mating Changes Sexually 586 Dimorphic Gene Expression in the Seed Beetle Callosobruchus maculatus. Genome 587 Biol Evol, 9(3), 677-699. https://doi.org/10.1093/gbe/evx029 588 Ingleby, F. C., Flis, I., & Morrow, E. H. (2014). Sex -biased gene expression and sexual conflict 589 throughout development. Cold Spring Harb Perspect Biol , 7(1), a017632. 590 https://doi.org/10.1101/cshperspect.a017632 591 Jacobs, A., & Elmer, K. R. (2021). Alternative splicing and gene expression play contrasting 592 roles in the parallel phenotypic evolution of a salmonid fish. Mol Ecol, 30(20), 4955-593 4969. https://doi.org/10.1111/mec.15817 594 Kalsotra, A., & Cooper, T. A. (2011). Functional consequences of developmentally regulated 595 alternative splicing. Nat Rev Genet, 12(10), 715-729. https://doi.org/10.1038/nrg3052 596 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 21 Kopp, A., Duncan, I., Godt, D., & Carroll, S. B. (2000). Genetic control and evolution of sexually 597 dimorphic characters in Drosophila. Nature, 408(6812), 553 -559. 598 https://doi.org/10.1038/35046017 599 Kraus, M. E., & Lis, J. T. (1994). The concentration of B52, an essential splicing factor and 600 regulator of splice site choice in vitro, is critical for Drosophila development. Mol Cell 601 Biol, 14(8), 5360-5370. https://doi.org/10.1128/mcb.14.8.5360-5370.1994 602 Lande, R. (1980). Sexual Dimorphism, Sexual Selection, and Adaptation in Polygenic 603 Characters. Evolution, 34(2), 292 -305. https://doi.org/10.1111/j.1558-604 5646.1980.tb04817.x 605 Law, J. H., & Crespi, B. J. (2002). The evolution of geographic parthenogenesis in Timema 606 walking-sticks. Mol Ecol , 11(8), 1471 -1489. https://doi.org/10.1046/j.1365-607 294x.2002.01547.x 608 Li, H., Janssens, J., De Waegeneer, M., Kolluru, S. S., Davie, K., Gardeux, V., Saelens, W., David, 609 F. P. A., Brbic, M., Spanier, K., Leskovec, J., McLaughlin, C. N., Xie, Q., Jones, R. C., 610 Brueckner, K., Shim, J., Tattikota, S. G., Schnorrer, F., Rust, K.,…Zinzen, R. P. (2022). Fly 611 Cell Atlas: A single -nucleus transcriptomic atlas of the adult fruit fly. Science, 612 375(6584), eabk2432. https://doi.org/10.1126/science.abk2432 613 Li, Y., Wang, G., Tian, J., Liu, H., Yang, H., Yi, Y., Wang, J., Shi, X., Jiang, F., Yao, B., & Zhang, Z. 614 (2012). Transcriptome analysis of the silkworm (Bombyx mori) by high -throughput 615 RNA sequencing. PLoS One , 7(8), e43713. 616 https://doi.org/10.1371/journal.pone.0043713 617 Lichilin, N., El Taher, A., & Bohne, A. (2021). Sex -biased gene expression and recent sex 618 chromosome turnover. Philos Trans R Soc Lond B Biol Sci , 376(1833), 20200107. 619 https://doi.org/10.1098/rstb.2020.0107 620 Liu, L. J., Zheng, H. Y., Jiang, F., Guo, W., & Zhou, S. T. (2014). Comparative transcriptional 621 analysis of asexual and sexual morphs reveals possible mechanisms in reproductive 622 polyphenism of the cotton aphid. PLoS One , 9(6), e99506. 623 https://doi.org/10.1371/journal.pone.0099506 624 Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and 625 dispersion for RNA -seq data with DESeq2. Genome Biol , 15(12), 550. 626 https://doi.org/10.1186/s13059-014-0550-8 627 Love, M. I., Soneson, C., & Patro, R. (2018). Swimming downstream: statistical analysis of 628 differential transcript usage following Salmon quantification. F1000Res, 7, 952. 629 https://doi.org/10.12688/f1000research.15398.3 630 Magnusson, K., Mendes, A. M., Windbichler, N., Papathanos, P. A., Nolan, T., Dottorini, T., 631 Rizzi, E., Christophides, G. K., & Crisanti, A. (2011). Transcription regulation of sex -632 biased genes during ontogeny in the malaria vector Anopheles gambiae. PLoS One, 633 6(6), e21572. https://doi.org/10.1371/journal.pone.0021572 634 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 22 Mank, J. E. (2017). The transcriptional architecture of phenotypic dimorphism. Nat Ecol Evol, 635 1(1), 6. https://doi.org/10.1038/s41559-016-0006 636 Mank, J. E., Hultin -Rosenberg, L., Zwahlen, M., & Ellegren, H. (2008). Pleiotropic constraint 637 hampers the resolution of sexual antagonism in vertebrate gene expression. Am Nat, 638 171(1), 35-43. https://doi.org/10.1086/523954 639 Mazin, P. V., Khaitovich, P., Cardoso-Moreira, M., & Kaessmann, H. (2021). Alternative splicing 640 during mammalian organ development. Nat Genet , 53(6), 925 -934. 641 https://doi.org/10.1038/s41588-021-00851-w 642 McIntyre, L. M., Bono, L. M., Genissel, A., Westerman, R., Junk, D., Telonis -Scott, M., 643 Harshman, L., Wayne, M. L., Kopp, A., & Nuzhdin, S. V. (2006). Sex-specific expression 644 of alternative transcripts in Drosophila. Genome Biol , 7(8), R79. 645 https://doi.org/10.1186/gb-2006-7-8-R79 646 Meisel, R. P. (2011). Towards a more nuanced understanding of the relationship between sex-647 biased gene expression and rates of protein-coding sequence evolution. Mol Biol Evol, 648 28(6), 1893-1900. https://doi.org/10.1093/molbev/msr010 649 Merkin, J., Russell, C., Chen, P., & Burge, C. B. (2012). Evolutionary dynamics of gene and 650 isoform regulation in Mammalian tissues. Science, 338(6114), 1593 -1599. 651 https://doi.org/10.1126/science.1228186 652 Naftaly, A. S., Pau, S., & White, M. A. (2021). Long -read RNA sequencing reveals widespread 653 sex-specific alternative splicing in threespine stickleback fish. Genome Res , 31(8), 654 1486-1497. https://doi.org/10.1101/gr.274282.120 655 Nilsen, T. W., & Graveley, B. R. (2010). Expansion of the eukaryotic proteome by alternative 656 splicing. Nature, 463(7280), 457-463. https://doi.org/10.1038/nature08909 657 Nowicka, M., & Robinson, M. D. (2016). DRIMSeq: a Dirichlet -multinomial framework for 658 multivariate count outcomes in genomics. F1000Res, 5, 1356. 659 https://doi.org/10.12688/f1000research.8900.2 660 Papa, F., Windbichler, N., Waterhouse, R. M., Cagnetti, A., D'Amato, R., Persampieri, T., 661 Lawniczak, M. K. N., Nolan, T., & Papathanos, P. A. (2017). Rapid evolution of female-662 biased genes among four species of Anopheles malaria mosquitoes. Genome Res , 663 27(9), 1536-1548. https://doi.org/10.1101/gr.217216.116 664 Pardo-Palacios, F. J., Arzalluz -Luque, A., Kondratova, L., Salguero, P., Mestre -Tomas, J., 665 Amorin, R., Estevan-Morio, E., Liu, T., Nanni, A., McIntyre, L., Tseng, E., & Conesa, A. 666 (2024). SQANTI3: curation of long -read transcriptomes for accurate identification of 667 known and novel isoforms. Nat Methods , 21(5), 793 -797. 668 https://doi.org/10.1038/s41592-024-02229-2 669 Parker, D. J., Bast, J., Jalvingh, K., Dumas, Z., Robinson-Rechavi, M., & Schwander, T. (2019a). 670 Repeated Evolution of Asexuality Involves Convergent Gene Expression Changes. Mol 671 Biol Evol, 36(2), 350-364. https://doi.org/10.1093/molbev/msy217 672 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 23 Parker, D. J., Bast, J., Jalvingh, K., Dumas, Z., Robinson-Rechavi, M., & Schwander, T. (2019b). 673 Sex-biased gene expression is repeatedly masculinized in asexual females. Nat 674 Commun, 10(1), 4638. https://doi.org/10.1038/s41467-019-12659-8 675 Parker, D. J., Jaron, K. S., Dumas, Z., Robinson -Rechavi, M., & Schwander, T. (2022). X 676 chromosomes show relaxed selection and complete somatic dosage compensation 677 across Timema stick insect species. J Evol Biol , 35(12), 1734 -1750. 678 https://doi.org/10.1111/jeb.14075 679 Parker, D. J., Dumas Z., Lencero, R. G., Aury, J. M., Labedan, M., Tran Van, P., Istace, B., Cruaud, 680 C., Labadie, K., Noel, B., Freitas, S., Djordjevic, J., & Schwander, T. (2025). Dosage 681 compensation and meiotic sex chromosome inactivation are maintained in the 682 absence of selection. bioRxiv, https://doi.org/10.1101/2025.07.27.667025 683 Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., & Kingsford, C. (2017). Salmon provides fast 684 and bias-aware quantification of transcript expression. Nat Methods, 14(4), 417-419. 685 https://doi.org/10.1038/nmeth.4197 686 Perry, J. C., Harrison, P. W., & Mank, J. E. (2014). The ontogeny and evolution of sex -biased 687 gene expression in Drosophila melanogaster. Mol Biol Evol , 31(5), 1206 -1219. 688 https://doi.org/10.1093/molbev/msu072 689 Price, J., Harrison, M. C., Hammond, R. L., Adams, S., Gutierrez -Marcos, J. F., & Mallon, E. B. 690 (2018). Alternative splicing associated with phenotypic plasticity in the bumble bee 691 Bombus terrestris. Mol Ecol, 27(4), 1036-1043. https://doi.org/10.1111/mec.14495 692 Proschel, M., Zhang, Z., & Parsch, J. (2006). Widespread adaptive evolution of Drosophila 693 genes with sex -biased expression. Genetics, 174(2), 893 -900. 694 https://doi.org/10.1534/genetics.106.058008 695 Qiu, C., Zhang, Y., Fan, Y. J., Pang, T. L., Su, Y., Zhan, S., & Xu, Y. Z. (2019). HITS -CLIP reveals 696 sex-differential RNA binding and alterative splicing regulation of SRm160 in 697 Drosophila. J Mol Cell Biol, 11(2), 170-181. https://doi.org/10.1093/jmcb/mjy029 698 Ranwez, V., Douzery, E. J. P., Cambon, C., Chantret, N., & Delsuc, F. (2018). MACSE v2: Toolkit 699 for the Alignment of Coding Sequences Accounting for Frameshifts and Stop Codons. 700 Mol Biol Evol, 35(10), 2582-2584. https://doi.org/10.1093/molbev/msy159 701 Rogers, T. F., Palmer, D. H., & Wright, A. E. (2021). Sex-Specific Selection Drives the Evolution 702 of Alternative Splicing in Birds. Mol Biol Evol , 38(2), 519 -530. 703 https://doi.org/10.1093/molbev/msaa242 704 Saito, K., Ishizu, H., Komai, M., Kotani, H., Kawamura, Y., Nishida, K. M., Siomi, H., & Siomi, M. 705 C. (2010). Roles for the Yb body components Armitage and Yb in primary piRNA 706 biogenesis in Drosophila. Genes Dev , 24(22), 2493 -2498. 707 https://doi.org/10.1101/gad.1989510 708 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 24 Sartain, C. V., Cui, J., Meisel, R. P., & Wolfner, M. F. (2011). The poly(A) polymerase GLD2 is 709 required for spermatogenesis in Drosophila melanogaster. Development, 138(8), 710 1619-1629. https://doi.org/10.1242/dev.059618 711 Sawyer, S. A., Parsch, J., Zhang, Z., & Hartl, D. L. (2007). Prevalence of positive selection among 712 nearly neutral amino acid replacements in Drosophila. Proc Natl Acad Sci U S A , 713 104(16), 6504-6510. https://doi.org/10.1073/pnas.0701572104 714 Schutt, C., & Nothiger, R. (2000). Structure, function and evolution of sex -determining 715 systems in Dipteran insects. Development, 127(4), 667 -677. 716 https://doi.org/10.1242/dev.127.4.667 717 Schwander, T., Crespi, B. J., Gries, R., & Gries, G. (2013). Neutral and selection -driven decay 718 of sexual traits in asexual stick insects. Proc Biol Sci , 280(1764), 20130823. 719 https://doi.org/10.1098/rspb.2013.0823 720 Schwander, T., Henry, L., & Crespi, B. J. (2011). Molecular evidence for ancient asexuality in 721 timema stick insects. Curr Biol , 21(13), 1129 -1134. 722 https://doi.org/10.1016/j.cub.2011.05.026 723 Simpson, S. J., Sword, G. A., & Lo, N. (2011). Polyphenism in insects. Curr Biol, 21(18), R738-724 749. https://doi.org/10.1016/j.cub.2011.06.006 725 Singh, A., & Agrawal, A. F. (2023). Two Forms of Sexual Dimorphism in Gene Expression in 726 Drosophila melanogaster: Their Coincidence and Evolutionary Genetics. Mol Biol Evol, 727 40(5). https://doi.org/10.1093/molbev/msad091 728 Soneson, C., Love, M. I., & Robinson, M. D. (2015). Differential analyses for RNA -seq: 729 transcript-level estimates improve gene -level inferences. F1000Res, 4, 1521. 730 https://doi.org/10.12688/f1000research.7563.2 731 Soumillon, M., Necsulea, A., Weier, M., Brawand, D., Zhang, X., Gu, H., Barthes, P., Kokkinaki, 732 M., Nef, S., Gnirke, A., Dym, M., de Massy, B., Mikkelsen, T. S., & Kaessmann, H. (2013). 733 Cellular source and mechanisms of high transcriptome complexity in the mammalian 734 testis. Cell Rep, 3(6), 2179-2190. https://doi.org/10.1016/j.celrep.2013.05.031 735 Tardaguila, M., de la Fuente, L., Marti, C., Pereira, C., Pardo-Palacios, F. J., Del Risco, H., Ferrell, 736 M., Mellado, M., Macchietto, M., Verheggen, K., Edelmann, M., Ezkurdia, I., Vazquez, 737 J., Tress, M., Mortazavi, A., Martens, L., Rodriguez -Navarro, S., Moreno-Manzano, V., 738 & Conesa, A. (2018). SQANTI: extensive characterization of long -read transcript 739 sequences for quality control in full -length transcriptome identification and 740 quantification. Genome Res, 28(3), 396-411. https://doi.org/10.1101/gr.222976.117 741 Tong, L., Chen, X., Wang, W., Xiao, Y., Yu, J., Lu, H., & Cui, F. (2022). Alternative Splicing 742 Landscape of Small Brown Planthopper and Different Response of JNK2 Isoforms to 743 Rice Stripe Virus Infection. J Virol , 96(2), e0171521. 744 https://doi.org/10.1128/JVI.01715-21 745 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 25 Toubiana, W., Armisen, D., Dechaud, C., Arbore, R., & Khila, A. (2021). Impact of male trait 746 exaggeration on sex -biased gene expression and genome architecture in a water 747 strider. BMC Biol, 19(1), 89. https://doi.org/10.1186/s12915-021-01021-4 748 Trost, N., Mbengue, N., & Kaessmann, H. (2023). The molecular evolution of mammalian 749 spermatogenesis. Cells Dev , 175, 203865. 750 https://doi.org/10.1016/j.cdev.2023.203865 751 Van den Berge, K., Soneson, C., Robinson, M. D., & Clement, L. (2017). stageR: a general stage-752 wise method for controlling the gene -level false discovery rate in differential 753 expression and differential transcript usage. Genome Biol , 18(1), 151. 754 https://doi.org/10.1186/s13059-017-1277-0 755 Verta, J. P., & Jacobs, A. (2022). The role of alternative splicing in adaptation and evolution. 756 Trends Ecol Evol, 37(4), 299-308. https://doi.org/10.1016/j.tree.2021.11.010 757 Wan, Y., & Larson, D. R. (2018). Splicing heterogeneity: separating signal from noise. Genome 758 Biol, 19(1), 86. https://doi.org/10.1186/s13059-018-1467-4 759 Wang, E. T., Sandberg, R., Luo, S., Khrebtukova, I., Zhang, L., Mayr, C., Kingsmore, S. F., 760 Schroth, G. P., & Burge, C. B. (2008). Alternative isoform regulation in human tissue 761 transcriptomes. Nature, 456(7221), 470-476. https://doi.org/10.1038/nature07509 762 Wang, X., Xu, X., Lu, X., Zhang, Y., & Pan, W. (2015). Transcriptome Bioinformatical Analysis 763 of Vertebrate Stages of Schistosoma japonicum Reveals Alternative Splicing Events. 764 PLoS One, 10(9), e0138470. https://doi.org/10.1371/journal.pone.0138470 765 Wang, Z., & Burge, C. B. (2008). Splicing regulation: from a parts list of regulatory elements to 766 an integrated splicing code. RNA, 14(5), 802-813. https://doi.org/10.1261/rna.876308 767 Whittle, C. A., & Extavour, C. G. (2019). Selection shapes turnover and magnitude of sex -768 biased expression in Drosophila gonads. BMC Evol Biol , 19(1), 60. 769 https://doi.org/10.1186/s12862-019-1377-4 770 Whyte, W. L., Irick, H., Arbel, T., Yasuda, G., French, R. L., Falk, D. R., & Hawley, R. S. (1993). 771 The genetic analysis of achiasmate segregation in Drosophila melanogaster. III. The 772 wild-type product of the Axs gene is required for the meiotic segregation of 773 achiasmate homologs. Genetics, 134(3), 825 -835. 774 https://doi.org/10.1093/genetics/134.3.825 775 Yanai, I., Benjamin, H., Shmoish, M., Chalifa-Caspi, V., Shklar, M., Ophir, R., Bar-Even, A., Horn-776 Saban, S., Safran, M., Domany, E., Lancet, D., & Shmueli, O. (2005). Genome -wide 777 midrange transcription profiles reveal expression level relationships in human tissue 778 specification. Bioinformatics, 21(5), 650-659. 779 https://doi.org/10.1093/bioinformatics/bti042 780 Yang, Z. (2007). PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol, 24(8), 781 1586-1591. https://doi.org/10.1093/molbev/msm088 782 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 26 Yeo, G. W., Van Nostrand, E. L., & Liang, T. Y. (2007). Discovery and analysis of evolutionarily 783 conserved intronic splicing regulatory elements. PLoS Genet , 3(5), e85. 784 https://doi.org/10.1371/journal.pgen.0030085 785 Zheng, S. Y., Pan, L. X., Cheng, F. P., Jin, M. J., & Wang, Z. L. (2023). A Global Survey of the Full-786 Length Transcriptome of Apis mellifera by Single-Molecule Long-Read Sequencing. Int 787 J Mol Sci, 24(6). https://doi.org/10.3390/ijms24065827 788 789 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 27 Figures 790 791 792 Figure 1. Transcriptome complexity across tissues, between sexes, and between 793 reproductive modes. (A) Average n umber of isoforms per gene for each femur, gut, and 794 gonad sample across all species. Only genes expressed in all three tissues are included , and 795 genes with only one isoform in all tissues have been excluded for ease of comparison. (B) The 796 ratio between males and females in the average number of isoforms per gen e. For each 797 tissue, only genes expressed in both sexes are included , and genes with only one isoform in 798 both males and females have been excluded. (C) The ratio between sexual and 799 parthenogenetic females, and separately males, in the average number of isoforms per gene. 800 For each tissue, only genes expressed in both sexual and parthenogenetic individuals are 801 included, and genes with only one isoform have been excluded. In all panels , results from 802 pairwise Wilcoxon rank -sum tests are represented by lowercase letters, where boxplots 803 labelled with di stinct letters are significantly different ( p < 0.05), while boxplots sharing a 804 letter are non-significant. 805 806 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 28 807 Figure 2. Differential gene expression and differential transcript usage between males and 808 females in Timema sexual species. Number of genes with a significant difference in 809 expression (DE) and in transcript usage patterns ( DS) between males and females in femur, 810 gut, and gonad tissues of T. poppense (Tps), T. californicum (Tcm), T. cristinae (Tce), T. podura 811 (Tpa) and T. bartmani (Tbi). The numbers next to the bars represent the percentage of DE and 812 DS genes out of the total number of genes that could be tested in each respective analysis. 813 814 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 29 815 Figure 3. Intersection of orthologous gonad DE genes (A) and DS genes (B) across the five 816 sexual species. The number of DE and DS genes shared across all species is shown in orange. 817 818 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 30 819 Figure 4. Overlap, tissue specificity, and rates of coding-sequence evolution for DE and DS 820 genes in the gonad tissue of sexual Timema species. (A) Venn diagrams showing o verlap 821 between DS and DS genes . (B) Tissue specificity (tau) estimates based on expression 822 information from seven tissues , where values closer to 1 represent higher tissue -specific 823 expression, while values closer to 0 represent broader expression across tissues. In all 824 comparisons, DS genes have a significantly lower tau value compared to DE genes as 825 estimated by Wilcoxon rank -sum tests. (C) Estimates of the ratio between nonsynonymous 826 (dN) and synonymous substitutions ( dS). (D) Estimates of nonsynonymous substitution rates. 827 In both (C) and (D) 95% confidence intervals are shown and significance values are based on 828 1,000 replicates permutation tests. 829 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 31 830 Figure 5. Changes in gonadal sex differences in splicing following transition to asexuality. 831 (A) Schematic on the potential patterns of isoform abundance in asexual females for genes 832 identified as differentially spliced between males and females in the sexual species. Asexual 833 females may maintain the isoform abundance patterns found in sexual females (left), may 834 increase abundance of sexual male -biased isoforms and decrease expression for sexual 835 female-biased isoforms (middle), or may show no correlation in splicing patterns between 836 either sexual females or males (right), indicative of random changes in splicing patterns . 837 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 32 Hexagonal binning plot s showing the relationship between the sex difference in isoform 838 proportions for DS genes identified in the sexual species and the isoform abundance for the 839 same genes in (B) asexual females and (C) asexual males for two Timema species pairs. The 840 plot areas are partitioned into bins, where the darker the color shading of each bin, the higher 841 the number of unique genes. 842 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint 33 Tables 843 Table 1. Key genes with conserved differential splicing pattern across the Timema 844 phylogeny. The gene names and their functions are based on the Drosophila ortholog 845 annotations. 846 Gene name Gene description Relevant functions and pathways SRm160 serine/arginine repeat protein splicing coactivator involved in somatic sex determination (Qiu et al., 2019) B52 serine/arginine-rich splicing factor splicing factor; expressed in ovaries and developing eggs (Fear et al., 2015; Kraus & Lis, 1994) MBL muscleblind major regulator of alternative splicing; involved in tissue differentiation including muscle and neuronal differentiation (Bargiela et al., 2014) ARMI armitage involved in the localization of Piwi into Yb bodies (Saito et al., 2010) AXS aberrant X segregation chromosome segregation during meiosis (Whyte et al., 1993) PBL pebble role in cytokinesis (Deduyer et al., 2025) GLD2 germ line development 2 mRNA translation control; regulator of late spermatogenesis (Sartain et al., 2011) 847 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted February 22, 2026. ; https://doi.org/10.64898/2026.02.21.707240doi: bioRxiv preprint