Subgenome dominance shapes novel gene evolution in the decaploid pitcher plant Nepenthes gracilis

Subgenome dominance after whole-genome duplication generates distinction in gene number and expression at the level of chromosome sets, but it remains unclear how this process may be involved in evolutionary novelty. Here, we generated a chromosome-scale genome assembly of the Asian pitcher plant Nepenthes gracilis to analyze how its novel traits (dioecy and carnivorous pitcher leaves) are linked to genomic evolution. We found a decaploid karyotype with five complete sets of syntenic chromosomes (2n = 10x = 80) yet with a clear indication of subgenome dominance and highly diploidized gene contents. The male-linked and pericentromerically located region on the putative sex chromosome was identified in a recessive subgenome and was found to harbor three transcription factors involved in flower and pollen development, including a likely neofunctionalized LEAFY duplicate. Transcriptomic and syntenic analyses of carnivory-related genes suggested that the paleopolyploidization events seeded genes that subsequently formed tandem clusters in recessive subgenomes with specific expression in the digestive zone of the pitcher, where specialized cells digest prey and absorb derived nutrients. Novel gene evolution in recessive subgenomes is likely to be prevalent because duplicates were enriched with Nepenthes-specific genes with tissue-specific expression, including those expressed in trapping pitchers. Thus, subgenome dominance likely contributed to evolutionary novelty by allowing recessive subgenomes experiencing relaxed purifying selection to serve as a preferred host of novel tissue-specific duplicates. Our results provide insight into how polyploids, which may frequently be evolutionary dead-ends, have given rise to novel traits in exceptionally thriving high-ploidy lineages.


Introduction 67
Novel phenotypes are often the result of the emergence of new genes that have undergone 68 functional divergence after gene duplication. Polyploidy, which involves the duplication of entire 69 genomes, has been argued to have great potential to drive adaptive evolution (Soltis et al., 2009). 70 However, polyploid lineages are sometimes considered evolutionary dead-ends, with only short-term 71 adaptive potential and high long-term extinction rates ( Van de Peer et al., 2017). There are, however, 72 some exceptions to this trend, such as flowering plants (angiosperms) (AMBORELLA GENOME 73 PROJECT et al., 2013;Chanderbali et al., 2022) and vertebrates (Dehal and Boore, 2005), which have 74 thrived following multiple rounds of polyploidization early in their evolution. 75 After polyploidization, especially following a hybridization event (Edger et al., 2017;Li et al., 76 2021), the coexistence of more than one haploid chromosome set often engenders subgenome 77 dominance, in which one chromosome set is preferentially expressed over the other (Cheng et al., 2018). 78 This is one of the early changes that occur after polyploidization, and the differential retention of 79 duplicate genes on dominant versus recessive subgenomes can render versatile adaptive prospects 80 between these subgenomes over time. Although such a genome-wide effect may constrain the long-81 term evolutionary paths of a large number of genes, it remains unclear how subgenome dominance may 82 be involved in the evolution of novel gene functions. To study the long-term effects of subgenome 83 dominance, it is necessary to examine genomes including relatively ancient polyploidization events, yet 84 with well-preserved subgenome structures. 85 The Old World pitcher plant Nepenthes offers a unique opportunity to study the relationship 86 between polyploidization-related genomic features and the evolution of novel genes. This genus of 87 carnivorous plants includes at least 169 species of perennial vines (Cross et al., 2020) and is known for 88 its highly modified leaves (pitchers) adapted for trapping and obtaining nutrients from arthropod prey 89 (Ellison and Adamec, 2018). In addition to carnivory, Nepenthes possess other atypical traits, such as 90 dioecy. Less than 1% of all angiosperms are carnivorous (Cross et al., 2020), and approximately 6% of 91 angiosperms are dioecious, possessing distinct male and female individuals (Renner and Ricklefs, 92 1995). No carnivorous plant species outside of Nepenthes are known to be dioecious. Thus, no other 93 plant lineages possess this unique combination of traits. Nepenthes has undergone ancient whole-94 genome duplications (WGDs), as indicated by transcriptomic signatures (Walker et al., 2017;Yang et 95 al., 2018;Palfalvi et al., 2020), and exhibits highly stable karyotypes among many species (Heubl and 96 Wistuba, 1997). 97 In this study, we present a chromosome-scale genome assembly of Nepenthes gracilis and 98 demonstrate that the Nepenthes genome is decaploid, bearing five subgenomes with a clear signature of 99 subgenome dominance. Our analysis suggests that recessive subgenomes played a crucial role in 100 permitting the adaptive evolution of tissue-specific genes associated with the unique biology of 101 Nepenthes, and we discuss how gene divergence occurs under the influence of WGDs, subgenome 102 dominance, and subsequent small-scale duplications (SSDs). 103 104 105 ( Supplementary Fig. 5), indicates that Nepenthes has a complex polyploid background, having 153 undergone at least two sequential WGDs following the gamma hexaploidization event that occurred in 154 the common ancestor of all extant core eudicots (Jaillon et al., 2007). This history has also included the 155 addition of a fifth subgenome at some point during the evolution of the Nepenthes lineage, possibly 156 after splitting from its closest relatives ( Supplementary Fig. 4). Given that the chromosome number is 157 stable in Nepenthes, including in the earliest diverged species N. pervillei (Heubl and Wistuba, 1997), 158 decaploidy must have been established before the diversification of extant species (6.4-18.2 MYA 159 (Scharmann et al., 2021)). This suggests the possibility that these ancient WGDs may have played a 160 role in the evolution of novel traits in Nepenthes.

162
Biased gene fractionation and subgenome dominance. Despite having a complete decaploid 163 karyotype (n = 5x = 40 and 2n = 10x = 80), the gene content in Nepenthes is highly fractionated, as is 164 evident in a syntenic comparison with grape (Vitis vinifera), which has maintained its ploidy level since 165 the gamma hexaploidization event (i.e., 1x in Vitis versus 5x in Nepenthes) (Jaillon et al., 2007) 166 ( Fig. 2b). In total, 94.5% (5,002/5,293) of 100-gene genomic windows in the Nepenthes genome 167 retained fewer than 1.5 syntelog copies on average (compared to the ploidy-based expectation of 5.0 168 copies). While fractionation is advanced overall, there is significant heterogeneity among syntenic 169 chromosome groups (Fig. 2b). For example, the syntenic group corresponding to grape chromosome 6 170 has an average of 1.4 syntelog copies, while the group corresponding to grape chromosome 16 has an 171 average of only 1.1 copies. Homeologous Nepenthes chromosome groups that are highly fractionated 172 tend to show a clear distinction between one dominant, gene-rich chromosome and four recessive, gene-173 poor chromosomes (i.e., fractionation bias (Joyce et al., 2017)). This distinction is less obvious among 174 groups that are less fractionated ( Fig. 2c and Supplementary Fig. 5). In fact, the degree of fractionation 175 bias is significantly correlated with the degree of diploidization (Fig. 2d). The dominant subgenome 176 carries 56% of the 1x equivalent syntelog set (average of 5,293 sliding 100-gene windows), and every 177 chromosome in the recessive subgenomes has at least 46 single-copy syntelogs (i.e., not detected in the 178 other chromosomes) ( Supplementary Fig. 6). Moreover, dominant chromosomes tend to be larger in 179 assembled scaffold sizes (Fig. 2e). These results suggest that the establishment of dominant and 180 recessive subgenomes played a crucial role in efficient gene fractionation and in the differentiation of 181 chromosomal sizes. 182 In line with subgenome fractionation dominance, genes in the dominant subgenome tend to 183 show higher expression levels compared with corresponding syntelogs in recessive subgenomes 184 (Fig. 2f). Thus, the dominant/recessive distinction is clear in both gene retention and gene expression. 185 It is noteworthy that, despite the clear chromosome-level expression dominance, 47% (972/2,055) of 186 syntelog pairs showed higher expression levels in copies on the recessive subgenomes (Fig. 2g), 187 showing a significant contribution to gene activity by these subgenomes as well. 188 The odd-numbered subgenomes of Nepenthes with clear subgenome dominance strongly 189 suggest an allopolyploid history underlying this dominance pattern (Alger and Edger, 2020). Our 190 attempts to use subgenome-specific k-mers identified by SubPhaser, which has successfully phased the 191 subgenomes of many well-known allopolyploids (Jia et al., 2022), did not result in complete subgenome 192 phasing in Nepenthes ( Supplementary Fig. 7). Disentangling exact polyploidization history is 193 sometimes difficult, especially in high-ploidy lineages (e.g., decaploid cotton (Wang et al., 2016)). 194 Chromosome-scale syntenic comparisons between Nepenthes and related lineages, once they are 195 available, will be necessary to discern among alternative scenarios.

197
Sex chromosome evolution. Individual Nepenthes plants exclusively produce male or female flowers 198 (dioecy), in contrast to the hermaphroditic flowers of their closest relatives and most other angiosperms 199 (Fig. 1b). An XY-type sex determination system has independently evolved in this genus (Scharmann 200 et al., 2019), but the sex chromosome has so far remained unidentified. To locate the Y-linked 201 chromosomal region (YLR) in our genome assembly of a male N. gracilis, we mapped sequencing reads 202 of double digest restriction-site associated DNA (ddRAD-seq) from wild, sex-identified N. gracilis 203 individuals (Scharmann et al., 2019). The result showed a clear, sharply delimited signal of male-204 specific sequences in a 1 Mb region on chromosome 20 ( Fig. 3a; Supplementary Fig. 8). Chromosome 205 20 is a member of one recessive subgenome (Fig. 2a), suggesting that the sex chromosome evolved 206 from an ancestral chromosome that experienced relaxed purifying selection. Interestingly, the YLR has 207 no homologous region on the corresponding X-chromosome in the female genome assembly 208 ( Supplementary Fig. 9), implying that it is fully hemizygous, similar to the independently evolved YLR 209 of Asparagus (Harkess et al., 2020).

210
A major unresolved problem in sex chromosome research is the evolution of suppressed 211 recombination. The lack of recombination throughout the Nepenthes YLR is implicit in the population-212 resequencing approach we used for the delimitation of the YLR. Although hemizygosity would in itself 213 imply a lack of recombination, an alternative explanation is that the Nepenthes YLR is recombination-214 depauperate due to adjacency to the satellite-rich putative centromere ( Fig. 3b and Supplementary Fig.  215 10). Pericentromeric regions, due to their heterochromatin and low recombination rates, might have 216 favored the formation of Y-linked regions and other supergenes, as suggested in other plant and animal 217 systems (Schwander et al., 2014;Li et al., 2015;Horiuchi et al., 2022;Potente et al., 2022;Akagi et al., 218 2023). 219 220 Male-specific genes in the Y-linked chromosomal region. Notably, the male-specific region, which 221 spans approximately 1 Mb in the 16.7-Mb chromosome, contains DYSFUNCTIONAL TAPETUM 1 222 (DYT1), which is the only fully male-linked gene known to date in Nepenthes (Scharmann et al., 2019) 223 (Fig. 3b). DYT1 encodes a bHLH transcription factor that, in Arabidopsis, functions in the cell 224 maturation of tapetum, the nutritive cell layer that aids microsporogenesis in the developing anther 225 (Zhang et al., 2006). An ortholog of MALE MEIOCYTE DEATH 1 (MMD1), which encodes a PHD-226 finger transcription factor whose loss causes male meiotic defects (Yang et al., 2003), is another 227 transcription factor gene located in the male-specific region. While the characterized functions of DYT1 228 and MMD1 suggest their involvement in microspore development, which begins later in anther 229 development, female flowers in Nepenthes not only lack microspores but also, almost entirely, staminal 230 structures (Subramanyam and Narayana, 1971). We, therefore, hypothesized that another gene upstream 231 of DYT1 and MMD1 locates to the male-specific region and determines floral sex. Consistent with this 232 idea, we found a male-specific copy (LFY-Y) of the LEAFY (LFY) gene, which in hermaphroditic 233 angiosperms encodes a plant-specific transcription factor that assigns the floral fate of meristems 234 (Moyroud et al., 2010). LFY is one of the earliest expressed genes in flower development where it acts 235 as a master regulator. Loss of LFY gene function converts lateral floral organs into leaf-like identity 236 (Moyroud et al., 2010). 237 In Nepenthes, the three male-specific transcription factor genes (DYT1, MMD1, and LFY-Y) are 238 embedded in a tandem cluster of long interspersed nuclear elements (LINEs), a group of 239 retrotransposons ( Fig. 3b and Supplementary Fig. 8). Transposable elements have not only accumulated 240 in these intergenic regions but also in the introns of LFY-Y itself ( Supplementary Fig. 11). It seems 241 likely that this transposable element accumulation is associated with the suppressed recombination 242 discussed above.

243
LFY is well known to be a highly conserved gene in angiosperms, where it is maintained as a 244 single copy in most species (Moyroud et al., 2009). The Nepenthes genome harbors what is likely to be 245 the principal LFY copy on an autosome (LFY-A on chromosome 3), which is not homeologous to the 246 putative sex chromosome. Gene phylogeny (Fig. 3c), chromosomal syntenic groups (Fig. 2a), and the 247 presence of introns ( Supplementary Fig. 11) suggest that the two LFY copies emerged by a lineage-248 specific SSD rather than via WGD in Nepenthes after the split of the other carnivorous lineages in 249 Caryophyllales, consistent with the emergence of dioecy in Nepenthaceae alone. This situation contrasts 250 with DYT1 and MMD1, both of which are maintained as single-copy genes in the decaploid genome of 251 Nepenthes. As expected, DYT1, MMD1, and LFY-Y were expressed in the developing buds of male 252 flowers but not in female buds (Fig. 3d,e). Furthermore, the male-specific expression of LFY-Y was 253 maintained until flower maturation ( Supplementary Fig. 12a). By contrast, we found no significant 254 difference in LFY-A expression between male and female buds. Although this expression analysis has 255 been performed by heterologous mapping to the N. gracilis genome using RNA-seq reads from related 256 species (see Methods), the read mapping rates were comfortably high (75.9-83.9% in 31 samples), and 257 we further confirmed by transcriptome assembly that the transcripts of DYT1, MMD1, and LFY-Y were 258 detected only in male samples . In the YLR, only these three genes could 259 be functionally annotated with sequence similarity against the UniProt database (Supplementary Table  260 5), and they showed the highest male/female transcript abundance ratios in developing buds among all 261 physically linked gene models in the region ( Supplementary Fig. 12e). In summary, the three 262 transcription factor genes in YLR could have been involved in establishing dioecy in Nepenthes.

264
Neofunctionalization of the male-specific LFY-Y gene copy. To examine how two LFY copies 265 acquired distinct expression patterns, we characterized cis-regulatory motifs in the promoters by 266 mapping their motif sequences using the JASPAR database (Castro-Mondragon et al., 2022). With the 267 false-discovery rate cutoff of 0.05, we detected 23 types of putative transcription factor-binding motifs 268 in the LFY-Y promoters (Supplementary Table 7), of which three motifs were specifically found in LFY-269 Y among promoter sequences from other species: i.e., motifs bound by ARF18, TCX3, or ZHD6 270 (Fig. 3f). In Arabidopsis, ARF18 is highly expressed in developing carpels (Fig. 3g) and encodes a 271 transcriptional repressor (Liu et al., 2015), whereas TCX3 and ZHD6 are expressed mainly in anthers 272 (Fig. 3g), suggesting roles for these regulators in modulating LFY-Y expression in both male and female 273 organs.

274
These LFY-Y-specific motifs (ARF18, TCX3, and ZHD6) were detected 1.0-1.5 kb upstream 275 from the start codon. This region may have been newly acquired by LFY-Y during sex chromosome 276 evolution. Within this promoter region, we also detected paired PISTILLATA-binding (PI-binding) 277 motifs only in Nepenthes LFY-Y, and interestingly, in LFY promoters of the comparator species spinach 278 (Spinacia oleracea). LFY in fact also upregulates PI expression in Arabidopsis (Honma and Goto, 279 2000). PI is a B-class MADS-box protein that specifies petal and stamen identity in Arabidopsis and 280 other plants (Theißen et al., 2016). Spinach evolved dioecy independently from Nepenthes ( Fig. 1b), 281 and the suppression of spinach PI converts male flowers into phenotypically normal female flowers 282 (Sather et al., 2010), suggesting a potentially convergent mechanism underlying the evolution of dioecy 283 among different Caryophyllales lineages through regulatory interactions between LFY, PI and other 284 transcription factors. Taken together, these results suggest that LFY-Y was neofunctionalized in terms 285 of expression pattern, likely as a consequence of the change in cis-regulatory elements. We detected 286 five amino acid substitutions specifically found in LFY-Y in its C-terminus DNA-binding domain 287 ( Supplementary Fig. 13), suggesting potential changes in protein properties in addition to the expression 288 patterns.

290
Pitcher-tissue-specific prey capture responses. Besides dioecy, the trap leaf organization also serves 291 as an evolutionary innovation in Nepenthes. A single Nepenthes leaf consists of segmented parts, 292 including a carnivorous pitcher trap, which is further elaborated by tissue differentiations along the 293 proximodistal axis, representing one of the most complex leaf shapes known among angiosperms 294 (Tsukaya, 2014). In Caryophyllales, plant carnivory evolved before the origin of Nepenthaceae 295 (Fig. 1b), presumably with the flypaper-type trapping mechanism and in a conventional leaf 296 organization (Albert et al., 1992;Heubl et al., 2006;Renner and Specht, 2011), with the pitcher 297 organization representing an evolutionary novelty that emerged in Nepenthes.

298
To examine whether gene expression patterns parallel unique pitcher tissue differentiation, we 299 performed RNA-seq experiments in six dissected leaf parts: from proximal to distal, flat part, tendril, 300 digestive zone, waxy zone, peristome, and lid ( Fig. 4a-b). For a cross-species comparison, we also 301 generated the pitcher tissue-specific transcriptomes of Cephalotus, which evolved pitcher leaves 302 independently from Nepenthes in an altogether different angiosperm order (Fig. 1b), thus representing 303 convergent evolution (Albert et al., 1992;Freund et al., 2022). We further integrated additional RNA-304 seq data from Dionaea (Bemm et al., 2016;Iosip et al., 2020;Procko et al., 2021) and Arabidopsis 305 (Klepikova et al., 2016) to form a 4-species expression level dataset normalized with trimmed mean of 306 M-values (TMM) (Robinson and Oshlack, 2010) for a total of 3,572 single-copy genes. Dimensionality 307 reduction analysis with different methods consistently showed that, in gene expression profiles, the 308 pitcher tissues of Nepenthes were distinct from conventional photosynthetic organs 309 ( Supplementary Fig. 14a). This contrasts with the pitcher leaves of Cephalotus, whose tissues showed 310 expression profiles largely overlapping with the photosynthetic organs, suggesting that pitcher tissues 311 in Nepenthes traps are much more strongly differentiated from ancestral photosynthetic leaf structures.

312
This view is consistent with the contrasting patterns of the tissue-specific photosynthetic activity in the 313 two pitcher plant lineages, with Nepenthes pitchers showing almost no photosynthetic assimilation 314 (Pavlovič et al., 2007;Pavlovič, 2011). 315 Next, we characterized tissue-specific prey capture responses in Nepenthes. Pitchers were 316 treated with an insect homogenate to mimic prey capture, and tissues were harvested after 24 h with the 317 expected onset of jasmonic acid-mediated (JA-mediated) prey responses (Yilamujiang et al., 2016). JA 318 is a known regulator of gene expression in plants and can mediate plant-herbivore and pathogen 319 interactions (Glazebrook, 2005;Erb et al., 2012), as well as prey-capture responses in the 320 Caryophyllales lineage of carnivorous plants including Nepenthes (Pavlovič and Mithöfer, 2019). As 321 anticipated, the number of significantly differentially expressed genes (DEGs, FDR < 0.05) was highest 322 in the digestive zone, where digestive/absorptive glands (Freund et al., 2022) are directly exposed to 323 the insect homogenate (Fig. 4b). A clear response of the digestive zone was further confirmed by 324 dimensionality reduction analysis ( Supplementary Fig. 14b). Significantly enriched Gene Ontology 325 (GO) terms (Supplementary Table 8; Supplementary Table 9; Supplementary Table 10;  326  Supplementary Table 11) among the up-regulated genes clearly indicated the activation of translation 327 machinery upon prey capture (Supplementary Table 8), which would be anticipated for rapid response 328 by the palette of necessary digestive enzymes. Other enriched terms include those presumably 329 associated with prey digestion, nutrient metabolism, and biological interactions (Fig. 4b). The second-330 highest number of DEGs was found in the tendril, which connects the pitcher to the rest of the plant 331 body. Interestingly, with the above threshold, we did not detect DEGs in other trapping leaf tissues, 332 including the flat part, where absorbed nutrients flow through. These results suggest that prey response 333 is well compartmentalized within a Nepenthes trapping leaf.

335
Tandem gene clusters with pitcher-tissue-specific expression. The leaf complexity of Nepenthes may 336 be reflected by considerable sub-or neofunctionalized duplicate gene expression in novel tissues, and 337 amplification of gene copies is often mediated by tandem gene duplications (Durand and Hoberman, 338 2006). We therefore searched for tandem gene clusters with high tissue specificity in gene expression, 339 measured by a metric called τ (Yanai et al., 2005), which outperformed alternative measures in a 340 benchmark study (Kryuchkova-Mostacci and Robinson-Rechavi, 2016). A conspicuous example we 341 detected was gene clusters of SENESCENCE-RELATED GENE 1 (SRG1) orthologs (Callard et al., 342 1996). This senescence marker gene encodes cytoplasm-localized melatonin 3-hydroxylase (M3H) and 343 produces cyclic 3-hydroxymelatonin, whose antioxidant activity is 15-fold higher than that of its 344 precursor melatonin, and promotes growth in Arabidopsis and rice (Lee et al., 2016;Choi and Back, 345 2019; Lee and Back, 2022). This family of genes forms tandem duplicates on separate Nepenthes 346 chromosomes within the same syntenic group (Fig. 4c), with the largest cluster (19 gene models with 347 15 cases with complete protein domain structure ( Supplementary Fig. 15)) on chromosome 20 (the 348 YLR-containing chromosome), which belongs to a recessive subgenome. Phylogenetic analyses 349 suggested that the first copy of this cluster originated as a WGD duplicate whose counterpart in the 350 dominant subgenome (chromosome 6) remains a single-copy gene (Nepgr007022). The tandem genes 351 in the dominant subgenome showed well-conserved microsynteny among eudicots (Supplementary Fig.  352 15), suggesting ancient origins (>131 MYA, Fig. 1b). Interestingly, the tandem array on chromosome 353 20 contains many genes with expression specific to the digestive zone in the pitcher, which forms the 354 interface to the digestive fluid, where peroxidases likely produce ROS to facilitate proteolysis for prey 355 degradation (Chia et al., 2004;Hatano and Hamada, 2012;Fukushima et al., 2017;Wal et al., 2022).

356
Given this characteristic expression, the SRG1 proteins may participate in scavenging cytotoxic reactive 357 oxygen species (ROS) produced during prey digestion and nutrient absorption.

358
To further examine whether other carnivory-related genes similarly formed tandem clusters, we 359 analyzed genes encoding digestive enzymes. Using the list of experimentally confirmed digestive fluid 360 proteins, many of which are digestive enzymes (Fukushima et al., 2017), we characterized their 361 phylogeny, synteny, and expression ( Supplementary Fig. 16). Among the 11 families of genes we 362 analyzed, we detected eight tandem gene clusters (six on recessive subgenomes) that were formed in 363 Nepenthes after its split from close relatives in Caryophyllales. Tandem clusters encoding typical 364 digestive enzymes (aspartic protease, class III peroxidase, glycoside hydrolase family 19 chitinase, 365 purple acid phosphatase, RNase T2, and β-1,3-glucanase) included genes whose transcript abundance 366 was highest in the digestive zone, with increased expression following the feeding treatment. These 367 results suggest that specific tandem cluster evolution is not restricted to SRG1 and may be found in 368 many other genes involved in Nepenthes carnivory.

370
Unequal subgenomic contributions for Nepenthes-specific genes. Since the male-specific region 371 harboring LFY-Y and many tandem gene clusters expressed in the digestive zone (SRG1 and digestive 372 enzyme genes) were found in recessive subgenomes rather than the dominant subgenome, subgenomes 373 may differentially serve as hosts of novel duplicated genes. However, these observations may also be 374 explained by chance, since there is a 70% probability that a novel gene will occur in the recessive 375 subgenome if the occurrence is proportional to the chromosome assembly size (528/753 Mb), and a 376 67.7% probability if it occurs proportionally to the number of genes (23,039/34,010 genes). Therefore, 377 statistical analysis was necessary. To this end, we analyzed lineage-specific genes in the Nepenthes 378 genome. First, DIAMOND BLASTP searches (Buchfink et al., 2021) were conducted against 20 other 379 plant genomes to identify Nepenthes genes that are most similar to another Nepenthes gene rather than 380 to genes from other species and thus likely to have emerged by gene duplication after the split from 381 those other lineages. Identified lineage-specific genes were significantly enriched with those having 382 specific expression in all analyzed tissues (Fig. 5b). Next, the types of duplication (i.e., WGD or SSD) 383 and involved subgenomes (dominant or recessive) were estimated by whether the best-hit gene was on 384 the same chromosome, on homeologous chromosomes, or on other chromosomes (Fig. 5a). Finally, the 385 frequency of each category in tissue-specific genes was compared with the overall average obtained 386 from all expressed genes. Interestingly, the recessive subgenomes tended to host more tissue-specific 387 genes than the dominant subgenome. Such distribution shifts from the null expectation with all 388 expressed genes were statistically significant in genes specifically expressed in six out of the seven 389 analyzed tissues (Fig. 5b), including pitcher tissues, whose differentiation evolved in Nepenthes to form 390 pitfall traps. The difference is likely to have emerged during functional divergence after gene 391 duplications because the frequency of tandem duplications themselves was comparable between 392 dominant and recessive subgenomes ( Supplementary Fig. 17).

393
The increased contribution to the tissue-specific genes by recessive subgenomes was most 394 pronounced among within-chromosome SSD duplicates, which are enriched in tandemly duplicated 395 genes. The rate of protein evolution (dN/dS) was higher in SSD duplicates than in WGD duplicates, 396 with the within-chromosome SSD duplicates showing the highest rates ( Supplementary Fig. 18).

397
Because there was no overall difference in dN/dS between dominant and recessive subgenomes, it is 398 possible that subgenome dominance may have only transiently relaxed purifying selection in the 399 recessive subgenomes, and that its signature is difficult to capture in the extant gene sequences, except 400 for the substantial gene losses via fractionation discussed above (Fig. 2 'response to oomycetes', whereas recessive tandems were enriched, for example, with GOs such as 419 'methyl jasmonate methylesterase activity' and 'salicylic acid glucosyltransferase (glucoside-forming) 420 activity', which could conceivably be associated with prey recognition pathways (Hedrich and 421 Fukushima, 2021).

422
Taken together, these results suggest that recessive subgenomes play an important role in 423 hosting novel tissue-specific genes that evolved through SSDs, possibly in an environment of relaxed 424 purifying selection compared to that in the dominant subgenome, and thereby to the unique biology of 425 Nepenthes, including dioecy ( Fig. 3) and carnivory (Fig. 4).

427
Discussion 428 In this study, we elucidated the decaploid structure of the Nepenthes genome and identified a 429 clear signature of 1:4 subgenome dominance (Fig. 2). We also highlighted how the four recessive 430 subgenomes have contributed to evolutionary novelties in Nepenthes (Fig. 1), especially in relation to 431 dioecy and carnivorous trapping leaves. 432 Our analysis of the organization of the male-specific chromosomal region (Fig. 3) suggested 433 that a series of mutational events may have led to the evolution of dioecy in Nepenthes 434 ( Supplementary Fig. 20a). A likely scenario is an evolutionary path from hermaphroditism via 435 gynodioecy towards full dioecy (Pannell and Jordan, 2022 to explain the evolution of the bisexual angiosperm flower from separate male and female reproductive 454 axes controlled by two distinct LFY copies in gymnosperms.

455
In addition, our analysis of tissue-specific genes in Nepenthes pitcher leaves ( Fig. 4) showed 456 how molecular evolution likely paralleled the increased complexity of tissue organization in Nepenthes. 457 As one of the most prominent examples, SRG1 genes, which may participate in scavenging ROS during 458 prey digestion and nutrient absorption, formed a massive tandem cluster in a recessive subgenome, and 459 many of its members acquired tissue-specific expression in Nepenthes-specific tissues, such as the 460 digestive zone. In contrast to the LFY duplication, the mode of functional divergence (i.e., neo-or 461 subfunctionalization) in the SRG1 cluster is not clear. However, it seems clear from gene expression 462 and phylogenetic relationship data ( Fig. 4c) that successive functional divergence has taken place 463 during cluster formation. These instances of novel genes (i.e., LFY-Y and SRG1) occurred in recessive 464 subgenomes, as do many other Nepenthes-specific genes that have acquired tissue-specific expression 465 (Fig. 5), including digestive enzyme genes ( Supplementary Fig. 16). This genome-wide trend (Fig. 5) 466 suggests that WGDs influenced the fates of subsequently produced SSD duplicates through recessivity 467 in a system with strongly divergent subgenome dominance patterns.

468
Our findings of novel gene duplicates suggest that the higher occurrence of tissue-specific genes 469 after SSDs in recessive subgenomes may have contributed to adaptive potential in Nepenthes, and thus, 470 to the maintenance of recessive subgenomes. Czech Republic). The plants were grown on a mixture of sphagnum moss and perlite in the greenhouse.

486
Since it was difficult to obtain fresh flower samples from N. gracilis, we obtained developing 487 inflorescences from cultivars and other species (Supplementary Table 4 High-molecular-weight genomic DNA isolation. Young leaf tissues of male and female Nepenthes 499 gracilis individuals from the wild were gathered, cleaned, and flash-frozen in liquid nitrogen, and then 500 stored at -80°C prior to extraction. About 10 g of flash-frozen tissue was used for high-molecular-weight 501 (HMW) genomic DNA isolation. The first step followed the BioNano NIBuffer nuclei isolation 502 protocol, in which frozen leaf tissue was homogenized in liquid nitrogen, with a subsequent nuclei lysis 503 step using IBTB buffer with spermine and spermidine added and filtered just before use. IBTB buffer 504 consists of Isolation Buffer (IB; 15 mM Tris, 10 mM EDTA, 130 mM KCI, 20 mM NaCl, 8% (w/v) 505 PVP-10, pH9.4) with 0.1% Triton X-100, and 7.5% (V/V) β-mercaptoethanol mixed in and chilled on 506 ice. The mixture of homogenized leaf tissue and IBTB buffer was strained to remove undissolved plant 507 tissue. 1% Triton X-100 was added to lyse the nuclei before centrifugation at 2000 ×g for 10 min to 508 pellet the nuclei. Once the nuclei pellet was obtained, we proceeded with CTAB 509 (cetyltrimethylammonium bromide) DNA extraction with modifications for the Oxford Nanopore 510 sequencer as described previously (Michael et al., 2018). The quality and concentration of HMW 511 genomic DNA was checked using Thermo Scientific™ NanoDrop™ Spectrophotometer, as well as on 512 agarose gel electrophoresis following standard protocols. The genomic DNA was further purified with 513 a Qiagen® Genomic-tip 500/G according to the protocol provided by the developer. 514 515 Genome sequencing. HMW DNA for a male and a female plant was used to generate a sequencing 516 library for use with Oxford Nanopore Technology (ONT) (SQK-LSK109; Oxford, England). The 517 resulting libraries were run on a PromethION sequencer running for 48 h. All bases were called on the 518 PromethION using Guppy v2.0, and the resulting fastq files were used for genome assembly. 519 520 Genome assembly and polishing. The fastq sequencing data was filtered for >10-kb reads using seqtk 521 v1.2 (https://github.com/lh3/seqtk). The filtered reads were assembled using WTDBG2 v2.2 (wtdbg2 -522 t64 -p19 -AS2 -e2 -L10000) (https://github.com/ruanjue/wtdbg2). A consensus genome assembly was 523 generated by mapping reads >10 kb to the assembly with minimap2 v0.2 524 (https://github.com/lh3/minimap2), and then running Racon v1.3.1 (https://github.com/isovic/racon); 525 the consensus process was repeated three times. The contig assembly was further polished using paired-526 end 2×150 Illumina sequence by first aligning the reads to the consensus assembly using minimap2 527 followed by running the assembly tool Pilon v1.18 (https://github.com/broadinstitute/pilon) three times 528 using the paired-end 2×150 Illumina sequence data. Purge Haplotigs v1.0.0 529 (https://bitbucket.org/mroachawri/purge_haplotigs/src/master/) was applied to both the female and 530 male Nanopore assemblies separately. The raw Illumina reads were aligned to the genome assemblies 531 using bwa mem v0.7.17 (https://github.com/lh3/bwa), and input files were prepared using SAMtools 532 v1.3 (https://github.com/samtools/samtools). Purge Haplotigs was then run using the prepared bam files 533 and genome assembly. 534 535 Hi-C scaffolding and syntenic path assembly. While Illumina-polished ONT-based wtdbg2 536 assemblies were generated independently for male and female specimens, we further improved the 537 contiguity of the male assembly using HiRise scaffolding of Chicago and Dovetail Hi-C libraries by  Table 1). 539 Heterologous Hi-C scaffolding of the female genome assembly was performed using the female wtdbg2 540 assembly and the male Hi-C sequencing data. A list of Hi-C contact positions for the female was 541 generated with Juicer v1.6 (https://github.com/aidenlab/juicer). This file was then used as input for 3D-542 DNA v170123 (https://github.com/aidenlab/3d-dna) to order and orient fragments of the genome 543 assembly. Because of the lack of detectable synteny to any regions in the male assembly, an unnaturally 544 large scaffold of the female Hi-C assembly was collapsed to its original contigs. We also employed 545 another approach to scaffolding the female ONT assembly. With the male Hi-C assembly as a reference, 546 we applied the syntenic path assembly to the female genome using the SynMap function of CoGe 547 (Haug-Baltzell et al., 2017). Assembly statistics are available in Supplementary Table 2. The scaffold 548 numbering in the final assembly corresponds to the chromosome numbers discussed in this paper. 549 550 Feeding experiment. Dried mealworms (Batch No. L400518, MultiFit Tiernahrungs GmbH) were 551 ground into a fine powder with a mortar and a pestle to prepare a 100 mg/ml homogenate in water. The 552 homogenate was subsequently centrifuged at 2,000 ×g for 30 sec, and the supernatant (mealworm 553 extract) was obtained. Upon the feeding experiment, we measured the total volume of digestive fluid 554 and added the mealworm extract adjusted to 10% of the total volume to ensure a comparable 555 concentration among pitchers with different sizes. The pitchers were then sealed using parafilm and left 556 in the greenhouse for 24 h. The same procedure was applied to control plants fed with water alone. 557 558 RNA extraction and transcriptome sequencing. The leaves of N. gracilis were dissected into the lid, 559 peristome, waxy zone, digestive zone, tendril, and flat part ( Supplementary Fig. 21). Digestive fluid 560 was discarded. After washing with water and drying with a paper towel, the tissues were immediately 561 frozen in liquid nitrogen. Tissues from five leaves (from multiple individuals) were pooled for one 562 replicate, and a total of three biological replicates were prepared. Root samples were collected from 563 untreated plants, and each replicate was derived from a single individual. Frozen samples were 564 homogenized in liquid nitrogen using mortar and pestle. RNA extraction was performed with PureLink 565 Plant RNA Reagent solution (Thermo Fisher) according to the manufacturer's instructions. The RNA 566 pellet was resuspended in RNase-free water at 4°C. After centrifuging at 14,000 rpm for 10 min at 4°C, 567 the solution was transferred to a new tube and purified using an RNeasy Mini kit (Qiagen). While the 568 above extraction method was used for Nepenthes, Ancistrocladus, Triphyophyllum, and Cephalotus 569 ( Supplementary Fig. 23) as well as leaves and roots of Drosophyllum, a modified CTAB protocol was 570 employed for Drosophyllum flowers. The frozen Drosophyllum flower samples were homogenized with 571 mortar and pestle. The ground sample was transferred to a pre-cooled 2-ml tube, and 0.75 ml of 572 preheated 2×CTAB buffer (65°C) was added. The tube was shaken vigorously by using a vortex mixer 573 and incubated in a thermoshaker at 1,400 rpm for 20 min at 60°C. An aliquot of 0.75 ml 574 chloroform:isoamyl alcohol (25:1) was added and mixed vigorously. The homogenate was centrifuged 575 at 15,000 rpm for 15 min at room temperature. The supernatant was then transferred to a new 2-ml tube. 576 For RNA precipitation, 1.5 ml of 100% ethanol was mixed with 60 µl of 3 M aqueous sodium acetate 577 and added to the supernatant, which was then shaken at 500 rpm for 60 min at room temperature. The 578 homogenate was centrifuged at 15,000 rpm for 20 min at room temperature. The supernatant was 579 discarded carefully so as not to disturb the pellet. The pellet was then washed with 1 ml of 75% ethanol. 580 The supernatant was discarded, and the pellet was air-dried for 15 min in a thermoshaker at 37°C. The 581 pellet was resuspended in 100 µl of RNase-free water at 4°C by shaking at 500 rpm for 15 min. The 582 samples were again centrifuged at 15,000 rpm for 5 min at 4°C. The solution was transferred to a new 583 tube, and the RNA was purified using the RNeasy Mini kit (Qiagen). The quality of RNA was examined 584 using Nanodrop (Thermo Fisher) and the Qubit IQ assay kit (Invitrogen). Total RNA was sent to 585 Novogene UK, where paired-end mRNA sequencing libraries were prepared using the Novogene NGS 586 RNA Library Prep Set (PT042). Briefly, after the poly-A enrichment and fragmentation, the first strand 587 cDNA was synthesized using random hexamer primers followed by the second strand cDNA synthesis, 588 end repair, A-tailing, adapter ligation, size selection for 250-300-bp insert size, amplification, and 589 purification. Libraries were paired-end sequenced for 150 bps with an Illumina Novaseq 6000 platform. 590 591 Transcriptome assembly. Transcriptome assemblies for A. abbreviatus, D. lusitanicum, and T. 592 peltatum were generated with the RNA-seq reads (Supplementary Table 4). The reads were 593 preprocessed with fastp v0.20.0 (https://github.com/OpenGene/fastp) and assembled using Trinity (see 594 Supplementary Table 6 Table 3). For further identification and annotation of repetitive elements, we used 604 EDTA v2.1.0 (https://github.com/oushujun/EDTA). In order to identify the overall repetitiveness of 605 genomes, we performed de novo repeat discovery with RepeatExplorer2 (Novák et al., 2020). We used 606 a repeat library obtained from the RepeatExplorer2 analysis of Illumina paired-end reads. All clusters 607 representing at least 0.005% of the genomes were manually checked, and the automated annotation was 608 corrected if needed. Contigs from the annotated clusters were used to build a repeat library. 609 Transposable element protein domains (Neumann et al., 2019) found in the assembled genomes were 610 annotated using the DANTE tool available from the RepeatExplorer2 Galaxy portal (https://galaxy-611 elixir.cerit-sc.cz/). Further repeat density distribution plots shown in Supplementary Fig. 10  Analysis of the male-specific region. The male-specific region was delimited on the male genome 669 assembly using re-sequencing data (ddRAD-seq, (Peterson et al., 2012)) of 11 male and ten female 670 individuals sampled from wild populations, including data from previous work (Scharmann et al.,671 2019). Male versus female read coverage was evaluated by mapping the read data to the genome with 672 bwa mem v0.7.17-r1188, and alignments were filtered with SAMtools v1.12 against non-primary and 673 supplementary alignments, and in the case of paired-end reads, enforcing the "properly-paired" status. 674 Read depth per sample was counted in genomic windows using bedtools v2.30.0 675 (https://github.com/arq5x/bedtools2), and compared between the sexes by taking the log2 of the ratio of 676 male and female normalized depth sums given by 677 where ( is the sum of male read counts and , is the sum of female read counts. Window-specific null 679 distributions for this statistic were obtained by 1,000 permutations of the sex assignment of the 21 680 individuals. Male-specific sequences were further evaluated by counting k-mers (. = 16) in the 681 ddRAD-seq data of males and females using KMC v3.1.0 (https://github.com/refresh-bio/KMC), 682 requiring k-mers to occur at least two times. Observed k-mers were classified as male-specific if they 683 were present in at least nine out of 11 male samples and absent in all ten female samples, using 684 kmc_tools. All possible alignments of such male-specific k-mers in the reference genome were found 685 by bwa mem, allowing at most one mismatch. The abundance of such k-mer alignments was counted in 686 genomic windows using bedtools. Window-specific null distributions for the abundance of male-687 specific k-mers were generated by 1,000 permutations of the sexes and repetition of the above 688 kmc_tools and alignment procedure. 689 690 Analysis of differential expression and GO enrichment.  Albert VA, Oppenheimer DG, Lindqvist C (2002)      Digestion & metabolism structural constituent of ribosome (134 genes) polysaccharide catabolic process (10 genes) chitin catabolic process (10 genes) serine-type carboxypeptidase activity (8 genes) isoleucine biosynthetic process (6 genes) sulfate assimilation (4 genes) threonine catabolic process (3 genes) hydro-lyase activity (3 genes) Interactions with insects & microbes response to oomycetes (5 genes) negative regulation of defense response to insect (4 genes) jasmonic acid hydrolase (4 genes) salicylic acid biosynthetic process (3 genes) cellular response to toxic substance (3 genes) Chromosome 38 (recessive) N e p g r 0 2 1 2 5 6 N e p g r 0 2 1 2 5 7 N e p g r 0 2 1 2 5 8 N e p g r 0 2 1 2 6 1 N e p g r 0 2 1 2 6 3 N e p g r 0 2 1 2 6 5 N e p g r 0 2 1 2 6 6 N e p g r 0 2 1 2 6 8 N e p g r 0 2 1 2 6 9 N e p g r 0 2 1 2 7 0 N e p g r 0 2 1 2 7 4 13.2M 13.4M 13.6M Chromosome 20 (recessive) Figure 5. Differential subgenomic contributions to novel gene evolution. in tissue-specific genes were compared with that of expressed genes using χ 2 tests controlled for false 1217 discovery rates (FDR  Schäferhoff et al., 2010;Yang et al., 2018;Yao et al., 2019) or to the ADD clade (Meimberg et al., 1248(Meimberg et al., 2000Cameron et al., 2002;Renner and Specht, 2011;Walker et al., 2017;Walker et al., 2018). This 1249 may be explained by the possibility that Nepenthaceae emerged by the allopolyploidization between 1250 two lineages: one close to Droseraceae and another close to the ADD clade. To test this possibility, we 1251 analyzed the 17 eudicot species (Supplementary Table 6) using Grampa v1.3.1 1252 (https://github.com/gwct/grampa), which accounts for individual gene tree topology to distinguish auto-1253 and allopolyploidy events (Thomas et al., 2017). To obtain a high-confidence set of gene trees, including 1254 lineage-specific gene duplication events that provide the signal of WGDs, we extracted BUSCO genes 1255 with the eudicot dataset of OrthoDB v10 (eudicots_odb10), including those labeled as duplicated. Using 1256 the translated sequences aligned with MAFFT and TrimAl, the ML trees were reconstructed using IQ-1257 TREE with the LG+R4 substitution model and were rooted by NOTUNG (Chen et al., 2000). The 1258 analysis of the 2,326 gene trees using Grampa with N. gracilis assigned as H1 yielded the best multi-1259 labeled tree which infers an allopolyploidization between the ancestral lineages of Droseraceae and the 1260 ADD clade, with 21 mapped duplication events that were judged as supporting the sister relationship 1261 between N. gracilis gene(s) and Droseraceae gene(s) ( Supplementary Table 17). However, manual 1262 inspection of the allopolyploidization-supporting 24 gene trees that contain multiple N. gracilis genes 1263 revealed prevalent N. gracilis-specific duplications and paralogy or misplacement of N. gracilis genes. 1264 As a result, no analyzed gene trees actually showed the topology expected from the allopolyploidization 1265 scenario: i.e., ((N. gracilis gene 1, ADD genes), (N. gracilis gene 1, Droseraceae genes)). 1266 Another aspect of the data that could potentially indicate a history of allopolyploidization is the 1267 topologies of phylogenies for genes that have returned to single-copy after the event. Here, we took 1268 advantage of the advanced diploidization which has occurred since polyploidization events took place 1269 in Nepenthes: the set of BUSCO genes (eudicots_odb10) is largely single-copy in Nepenthes 1270 ( Supplementary Fig. 3d). If Nepenthes involved allopolyploid hybridization between the ADD and 1271 Droseraceae lineages followed by re-diploidization, we would expect that single-copy genes place 1272 Nepenthes either as sister to the Droseraceae or sister to the ADD families (and their relative proportions 1273 possibly manifesting as subgenome dominance, see main text). Strong gene tree incongruence of this 1274 kind is indeed observed in the data (quartet support only ~0.45 in the ASTRAL species tree, 1275 Supplementary Fig. 3e). However, gene tree incongruence will not necessarily derive from 1276 hybridization but can also be caused by incomplete lineage sorting (ILS): i.e., rapid and nearly 1277 simultaneous speciation events among the three lineages (Nepenthaceae, Droseraceae, and ADD 1278 families). To distinguish between these alternatives, we applied Species Networks applying Quartets 1279 (SNaQ) (Solís-Lemus and Ané, 2016), an ILS-aware (coalescent) method to infer phylogenetic 1280 networks on protein sequences of single-copy genes in Nepenthes. Since the inference of phylogenetic 1281 networks is computationally intensive, the number of taxa on a gene tree must be kept relatively small. 1282 We, therefore, pruned the available gene trees to subsets of nine taxa, and in the case that taxa were 1283 represented by more than one gene copy, these taxa were discarded entirely (treated as missing data; 1284 consistent with the method used for species tree inference in the main text; Supplementary Fig. 24). We 1285 evaluated ten taxon subsets including the same seven ingroup taxa (Nepenthes, Triphyophyllum, 1286 Drosophyllum, Ancistrocladus, Dionaea, Aldrovanda, and Drosera), and two randomly chosen 1287 outgroup taxa. Starting from a bifurcating topology inferred by ASTRAL v5.7.8 (Zhang et al., 2018), 1288 phylogenetic networks with zero and then one hybrid edge were inferred and optimized using the SNaQ 1289 approach with PhyloNetworks v0.15.3 (https://github.com/crsl4/PhyloNetworks.jl) with 24 parallel 1290 attempts. While network optimization with one hybrid edge by SNaQ was successful in all ten taxon 1291 sets, Nepenthes did not descend from a hybridization event in any of them. Instead, it was usually the 1292 Droseraceae lineage which was constructed as hybridized with a ghost lineage directly descending from 1293 the ingroup stem (seven taxon sets), or the ADD clade as descending from such an event (two taxon 1294 sets), or the ADD lineage as hybridized with the Nepenthes lineage (one taxon set). Although these 1295 analyses do not support the hypothesis that Nepenthes is an allopolyploid hybrid involving ancient 1296 Droseraceae and/or the ancient ADD clade, it must be noted that there is as yet no satisfying explanation 1297 for the observed gene tree incongruence. In particular, neither the bifurcating species trees nor the 1298 optimized phylogenetic networks with one hybrid edge do fit the observed gene trees well 1299 (QuartetNetworkGoodnessFit v0.4.0 (Cai and Ané, 2021), all empirical P values < 10 -5 ). 1300 To be thorough, we analyzed gene location on subgenomes. If allopolyploidy occurred between 1301 the two distant lineages, genes on particular subgenomes should show the signal of particular 1302 phylogenetic affinity. A total of 33.4% (129/386) of the gene trees with the allopolyploidization signal 1303 (including those containing only one N. gracilis gene) are located on the dominant subgenome, which 1304 is comparable and not significantly different from the expectation with all genes (32. 3%, 10,971/22,621, 1305 418 gene models are on unanchored scaffolds and are thus excluded; χ 2 test, P = 0.79, χ 2 = 0.069). 1306 Therefore, our results do not support the possibility of the allopolyploidization events between the two 1307 distant lineages (i.e., ancestral Droseraceae and ancestral ADD clade). Still, there are phylogenetic 1308 signals that support the affinity of N. gracilis genes with Droseraceae, and this may be explained by 1309 ancient hybridization events before polyploidization. In this scenario, to explain the subgenomic 1310 frequencies of gene trees supporting alternative topologies, enough time would have been needed for 1311 allelic recombination to occur after hybridization and before polyploidization. Homeologous exchange 1312 ( Bird et al., 2018), if it did occur, would also contribute to gene shuffling between subgenomes. 1313 1314 • tau: Yanai's tau for the tissue specificity of gene expression obtained from TMM-FPKM values. 1362 • is_enough_tau: Whether the tau value is greater or equal to 0.67. 1363 •  Supplementary Fig. 2. The scaffold size distribution of the Nepenthes gracilis male genome 1413 assembly. Scaffolds smaller than 1 Mb are merged to the right. 1414 1415   scaffold2  scaffold3  scaffold4  scaffold5  scaffold6  scaffold7  scaffold8  scaffold9  scaffold10  scaffold11  scaffold1  scaffold12  scaffold13  scaffold14  scaffold15  scaffold16  scaffold17  scaffold18  scaffold19  scaffold20  scaffold21  scaffold22  scaffold23  scaffold24  scaffold25  scaffold26  scaffold27  scaffold28  scaffold29  scaffold30  scaffold31  scaffold32  scaffold33  scaffold34  scaffold35  scaffold36  scaffold37  scaffold38  scaffold39  scaffold40 22  24  23  38  35  36  37  29  20  15  5  19  18  14  11  40  27  17  12  8  7  4  3  6   1  9  13  16  25  26  32  30  31  2  10  38  39  33  28  35  5  19  18  11  17  15  36  37  29  20  14  21  24  23  22  34  40  27  12  8  7  3  6  4   2  10  1  16  9  13  32  21  30  34  24  25  26  31  38  23  22  39  40  35  33  28  36  29  37  5  19  18  20  17  11  15  14  7  6  12  8  3  4  27 5 chromosomes assigned to 5 subgenomes in specific k-mer detection (chromosomes 2, 1, 8, 11, and 12) 1, 8, 11, and 12) 10 chromosomes assigned to 2 subgenome groups (dominant/recessive) in specific k-mer detection (chr. 2 and 3 vs. 1, 8, 11, 12, 17, 23, 24, and 40) 20 chromosomes assigned to 2 subgenome groups (dominant/recessive) in specific k-mer detection (chr. 2, 3, 4, and 5 vs. 16 recessive chromosomes) Supplementary Fig. 8 outputs that are available for download from the following CoGe links: male Hi-C assembly versus 1478 female Hi-C assembly (https://genomevolution.org/r/1ir4m) and male Hi-C assembly versus female 1479 syntenic path assembly (https://genomevolution.org/r/1nqd6).    e   21051  21052  21053  21054  21055  21056  21057  21058  21059  21060  21061  21062  21063  21064  21065  21066  21067  21068  21069  21070  21071  21072  21073  21074  21075  21076  21077  21078  21079  21080  21081  21082  21083  21084  21085  21086  21087  21088  21089  21090  21091  21092  21093  21094  21095  21096  21097  21098  21099  21100  21101  21102  21103  21104  21105  21106  21107  21108  21109  21110  21111  21112  21113  21114  21115  21116  21117 Nepenthes gracilis gene ID      F  FF  W  T  TW WF  T  TF  W  FF  WF TFTF   F  LF  T  T  WF  F F  L  F  WF L  T  LF  TF  DF  PF  WF  L  LF  TF TF  R   Chr. 22 Tandemness index test statistic (Brunner and Munzel, 2000). Gland zone Supplementary Fig. 24. Generating single-copy gene alignments. All detected BUSCO genes 1616 marked as single-copy (S) or fragmented (F) were extracted, while those marked as duplicated (D) or 1617 missing (M) were treated as missing data. 1618