Nepenthes × ventrata Transcriptome Profiling Reveals a Similarity Between the Evolutionary Origins of Carnivorous Traps and Floral Organs

The emergence of the carnivory syndrome and traps in plants is one of the most intriguing questions in evolutionary biology. In the present study, we addressed it by comparative transcriptomics analysis of leaves and leaf-derived pitcher traps from a predatory plant Nepenthes ventricosa × Nepenthes alata. Pitchers were collected at three stages of development and a total of 12 transcriptomes were sequenced and assembled de novo. In comparison with leaves, pitchers at all developmental stages were found to be highly enriched with upregulated genes involved in stress response, specification of shoot apical meristem, biosynthesis of sucrose, wax/cutin, anthocyanins, and alkaloids, genes encoding digestive enzymes (proteases and oligosaccharide hydrolases), and flowering-related MADS-box genes. At the same time, photosynthesis-related genes in pitchers were transcriptionally downregulated. As the MADS-box genes are thought to be associated with the origin of flower organs from leaves, we suggest that Nepenthes species could have employed a similar pathway involving highly conserved MADS-domain transcription factors to develop a novel structure, pitcher-like trap, for capture and digestion of animal prey during the evolutionary transition to carnivory. The data obtained should clarify the molecular mechanisms of trap initiation and development and may contribute to solving the problem of its emergence in plants.


INTRODUCTION
Predator plants attract arthropods that serve both as pollinators for sexual reproduction and as prey captured and digested in traps (active or passive); traps represent modified leaves developed to support plant survival in the nutrient-deficient environment (Jürgens et al., 2012;Pavlovič and Saganová, 2015;Bartlett, 2017). It is considered that carnivory evolved independently nine times; in particular, passive pitfall traps have at least six autonomous origins in different orders of flowering plants (Givnish, 2015). To date, the list of green predators contains more than 650 species representing 19 genera, 12 families, and five orders of the flowering plants, both monocots and eudicots (Pavlovič and Saganová, 2015). The most diverse lineage of carnivorous plants, Caryophyllales, includes one of the most famous genera of pitcher plants, Nepenthes, which comprises 160-180 species as well as numerous natural and cultivated hybrids, inhabiting humid and sunny lowlands or tropical mountains deficient in nitrogenous nutrients soils (Murphy et al., 2020). Because of the large number of Nepenthes species residing in diverse habitats, the genus has become a model for studying plant adaptive radiation and speciation (Thorogood et al., 2018).
Nepenthes species develop typical photosynthetic leaves whose tips form a tendril as an extension of the midrib, to which a pitcher-like pitfall trap is attached (Owen and Lennon, 1999). The trap has a leaf-like lid (operculum), which initially seals the growing pitcher until it is ripe and ready to capture the prey; it also prevents fluid dilution during rain in mature traps (Bauer et al., 2008). Once open, the pitcher attracts the prey with rim (peristome)-located nectaries forming a slippery zone covered with a thick wax layer to catch and prevent the escape of the prey which then falls into a lower digestive zone containing multicellular glands that secrete an acidic fluid saturated with hydrolytic enzymes (chitinases and proteases) and antifungal factors (Gaume and Forterre, 2007;Scholz et al., 2010;Pavlovic and Mithöfer, 2019). The main pitcher characteristics attracting the prey are color patterns (UV fluorescence/visible wavebands), sweet fragrance of the secreted extrafloral nectar, and high CO 2 levels (Kurup et al., 2013;Baby et al., 2017).
Most studies on Nepenthes focus on the morphophysiological and functional characteristics of the pitcher, including biochemical composition of the secreted fluid in the digestive zone and transcriptomic and proteomic changes in response to prey and during early pitcher opening (Gorb et al., 2004;Hamada, 2008, 2012;Rottloff et al., 2016;Wan Zakaria et al., 2016Baby et al., 2017;Zulkapli et al., 2017;Goh et al., 2020); however, the data regarding the evolutionary transition from the non-carnivorous to carnivorous status in Nepenthes species are scarce. Comparative transcriptomics of carnivorous and non-carnivorous leaves in Cephalotus follicularis identified genetic signatures associated with prey attraction, capture, digestion, and nutrient absorption underlying the developmental switch between carnivorous and non-carnivorous leaf state (Fukushima et al., 2017). Recent transcriptomic analysis of genes controlling pitcher development in Nepenthes khasiana has revealed that ASYMMETRIC LEAVES 1 (AS1) and ERECTA (ER) contribute to the formation of the tendril from the midrib, whereas REVOLUTA (REV) may be positively associated with the initiation and maturation of the pitcher (Dkhar and Pareek, 2019). The most intriguing question is identification of molecular pathways that directed the emergence of the carnivory syndrome, as there are examples of both convergent (across unrelated lineages) and divergent (within genera) evolution of the trap (Thorogood et al., 2018). The pitcher-like trap has diverse morphology adapted to trapping and digestion of various preys. This is especially characteristic for the species of the Nepenthes genus, which is suggested to have originated from a Drosera-like predecessor (Givnish, 2015;Murphy et al., 2020). A comparative analysis of the genome of three carnivorous Droseraceae species showed that the evolution of predation in this family could be facilitated by a whole-genome duplication (WGD) in their last common ancestor, followed by diversification of the multiplied genes, including recruitment of root-specific genes to developing traps, massive loss of genes involved in non-carnivorous nutrition, and expansion of gene families associated with the attracting, catching, digesting, and utilizing prey (Palfalvi et al., 2020). However, the molecular mechanism employed by carnivorous plants to develop such a novel structure is still unclear.
In this study, we aimed to identify transcriptional signatures of the transition from the leaf to the mature pitfall trap by performing comparative transcriptomics of Nepenthes × ventrata leaves and pitchers at three developmental stages. The results revealed activation of genes associated with defense response, shoot apical meristem (SAM) organization, anthocyanin biosynthesis, and flowering, including MADS-domain transcription factors (TFs), in pitchers. Considering that MADS-box genes are thought to have been involved in the evolution of plant flower organs, we suggest that the emergence of the pitcher structure was also due to highly conserved MADS-domain TFs, which usually determine the identity of floral meristems and organs in extant plants.

Plant Materials and Growth Conditions
Nepenthes × ventrata Hort. ex Fleming (Nepenthaceae Dumort.) is a natural hybrid between two Nepenthes clade 2 species: N. ventricosa Blanco (Insignes clade) and N. alata Blanco (Graciiflora clade) endemic to the northern forests of the Philippines (Murphy et al., 2020). Plants were obtained from a nursery garden and grown in a greenhouse under controlled conditions (16/8-h light-dark cycle with light intensity of 150-200 µmol m −2 s −1 , 20/24 • C night/day temperature cycle, and 75-90% humidity). Samples of mature leaves at the middle section and pitchers at three developmental stages: primordial pitcher (or early pitcher) (∼1 cm), unopened young pitcher (3-4 cm), and open mature pitcher (or late pitcher) (7-11 cm; the first day of opening; there are no insects inside) were collected from the upper part of plant; three biological replicates were analyzed.

RNA Isolation, Library Preparation, and Transcriptome Sequencing
Total RNA was extracted from 300 mg of each of the 12 tissue samples using a modified CTAB-based technique (Filyushin et al., 2019)  Before assembly, adapter sequences were removed from the reads using cutadapt v1.17 (Martin, 2011), low quality reads were trimmed using sickle v1.33 1 , and paired overlapping reads were merged using FLASH v1.2.11. De novo transcriptome assembly was carried out using Trinity v2.6.5 (Haas et al., 2013) with default parameters, using all reads combined together.
Coding regions in the assembled transcripts were predicted by TransDecoder v 5.1.0 2 within the Trinity software and the putative transcripts were screened for homology to known protein-coding plant genes using Diamond v.0.9 homology search in NCBI-NR database, and the completeness of the resulting assembly was assessed using BUSCO v.5.0.0 (Seppey et al., 2019). Potential contamination was assessed by checking the taxon of the best homolog for each predicted protein using a diamond search in the Uniref90 database which contains most of publicly available protein sequences clustered with 90% identity.
Relative transcription levels of protein-coding genes were calculated by mapping clean reads on the assembled transcripts using Trinity scripts with RSEM (Li and Dewey, 2011).
Differential expression was analyzed using Trinity scripts with the edgeR method; genes were considered differentially expressed if changes in the expression levels computed by EdgeR were not less than two-fold [| log 2 (FC)| ≥ 1], and false discovery rate (FDR) was ≤0.05.
Phylogenetic analysis was performed using the Fast Minimum Evolution method in NCBI 13 and the Maximum Likelihood method based on the JTT matrix-based model in MEGA7.0 (Kumar et al., 2016).

Validation of RNA-Seq Data
Real-time quantitative PCR results for 13 selected genes were validated with RT-qPCR. First-strand cDNA was synthesized with the Reverse Transcription System (Promega, Madison, WI, United States) using an oligo-dT and quantified by fluorimetry. RT-qPCR was performed with SYBR Green and ROX RT-PCR mixture (Syntol, Moscow, Russia), 3 ng cDNA, and 10 µM gene-specific primers in a CFX96 Real-Time PCR Detection System (Bio-Rad Laboratories, United States) at the following cycling conditions: initial denaturation at 95 • C for 5 min, 40 cycles of denaturation at 95 • C for 15 s and annealing/synthesis at 62 • C for 50 s. As homologs of the two house-keeping genes, Ubiquitin (UBQ) and Elongation Factor (ELF) previously suggested for N. khasiana gene expression normalization (Dkhar and Pareek, 2019) did not show uniform expression, we normalized the data using two other reference genes, UBQ-ligase Praja-2 and Actin 7. Experiments were carried out in three biological and three technical replicates and statistical analysis was performed using GraphPad Prism version 7.02 (GraphPad, San Diego, CA, United States 14 ). P ≤ 0.05 was considered to indicate significant difference.

Nucleotide Sequence Accession Numbers
The raw sequences obtained in this study have been deposited in the NCBI Sequence Read Archive under the accession numbers SRX5724495-SRX5724506.

Nepenthes × ventrata Transcriptomes
Tissue samples of mature leaves and pitchers at three stages of development-early pitcher, young pitcher, and late pitcher (Figure 1) were analyzed by RNA-seq.
The reads of the four samples were combined and assembled. Overall, 245,999 transcripts (with size from 201 to 14,825 bp, total length 234.2 Mbp, and N50 size 1,602 bp) were generated. The completeness of the assembly was verified for 1,614 marker genes of the Embryophyta taxon using BUSCO, which revealed Arrows indicates the part of the organ used for transcriptomics. dz, digestive zone; sz, slippery zone; pe, peristome (rim); op, operculum (lid); te, tendril; le, leaf; MP, late pitcher; YP, young pitcher. Dashed lines indicate organ parts taken for analysis. Scale bar = 1 cm. that 93.9 and 4.0% of the marker genes were complete and partial, respectively.
Annotation of the assembled transcripts using diamond search (cut-off E-value 10 −5 ) against NCBI non-redundant database of plant protein-coding genes revealed that among 85,558 proteincoding sequences, 62,037 had significant homology to nonhypothetical plant proteins (Supplementary Table 2). The 85,558 sequences were used for further analysis.
We analyzed the taxonomy of the best homolog for each putative protein; only 135 proteins showed better homology with bacterial proteins. The total expression of these bacterial transcripts ranged from 0.10 to 0.24% depending on the samples; seven and one transcripts were identified as belonging to fungi and insects, respectively. Of the 53,400 proteins that matched in the Uniref90 database, 53,065 had plant proteins as the best homologs, which indicated an extremely low level of contamination.
Analysis with the Trinity software identified many groups of transcripts that may belong to the same gene or its isoforms, alternative splice variants, and paralogs, which could be explained by multiple paleopolyploidy events (at least seven genome duplications, according to Walker et al., 2017) in the evolutionary history of non-core Caryophyllales, including Nepenthes, as well as by the hybrid origin of N. × ventrata.

Gene Ontology Annotation
Gene ontology classification of the transcripts according to the function of the translated products revealed that 4,827 of them belonged to Biological Process (BP), 2,752-to Molecular Function (MF), and 1,020-to Cellular Component (CC) categories representing 97 sub-aspects (Supplementary Tables 4, 5). The BP sub-aspects with the maximum number of transcripts were linked to stress response, cellular component organization, and biosynthetic processes (Supplementary Table 4). Most transcripts belonged to BP terms associated with general regulation of plant development, especially leaf development, senescence, and morphogenesis of leaves, and with photosynthesis, including carotenoid metabolism (Supplementary Table 5). The pitcher primordium, a redcolored organ with waxy inner walls, defines the efficiency of insect catching and processing; therefore, transcripts were classified in terms related to SAM development, biosynthesis of wax and anthocyanins, and plant-insect interactions. Of a particular interest are transcripts relevant to flowering initiation and formation of inflorescences and flowers, including reproductive organs, nectaries, pollen, ovules, and seeds (Supplementary Table 5).
The transcripts of the most enriched MF category were associated with utilization of ATP energy and the uptake and transfer of metal ions to specific locations, which is one of the basic requirements in photosynthetic plants (Yruela, 2013). The CC-classified transcripts were found to encode mainly nuclear, membrane, cytoplasmic, and plasma membrane proteins (Supplementary Table 5).
To obtain preliminary information about activation or suppression of various functions in leaves or pitchers, differentially expressed genes (DEGs) were identified based on the relative expression change of >2 times between the transcriptomes (Supplementary Table 6) and classified in GO terms (Supplementary Table 7). The most enriched BP/MF/CC terms in early and young pitchers vs. leaves were related to embryo development and active developmental and growth processes, whereas in late pitcher vs. leaves they were linked to plant-type cell wall biogenesis and organization, and in leaves vs. early/young/mature pitcher-to photosynthesis (Supplementary Table 7 and Supplementary Figure 1). Late pitcher was similar to leaves in the upregulation of genes involved in oxidation-reduction and downregulation of those related to DNA replication. Late pitcher differed from all other transcriptomes by the downregulation of defense response genes. Figure 2 shows transcriptional suppression of genes associated with photosynthesis and leaf developmental processes (early/young/mature pitcher), and induction of those involved in cell division, differentiation, and development of the growing (early/young pitcher) pitcher.

Kyoto Encyclopedia of Genes and Genomes Annotation
Kyoto Encyclopedia of Genes and Genomes classification of transcripts in terms of KEGG PATHWAY (KP) and KEGG BRITE (KB; protein families) revealed 26,235 significant matches in the KP database, which were assigned to seven main classes (Supplementary Figure 2a), 45 sub-categories, and 401 pathways (Supplementary Table 8 Table 8). The presence of multi-transcript KP sub-categories "Environmental adaptation" and "Immune system" (Supplementary Figure 2c), as well as pathways "Plant hormone signal transduction, " "Plant-pathogen interaction, " and "MAPK signaling pathwayplant" (Supplementary Table 8) may reflect defense mechanisms and adaptation to nutrient deficiency in the carnivorous plant.
It should be mentioned that the plant data used by KEGG do not include characteristics of the carnivorous syndrome. Digestive fluid in Nepenthes pitcher contains proteases (such as nepenthesin), chitinases, glucanases, phosphatases, proteases, and pathogenesis-related proteins (Fukushima et al., 2017;Ravee et al., 2018;Pavlovic and Mithöfer, 2019). In view of this, we briefly reviewed KP "Digestive system" and found that all our genes of interest were absent there and scattered throughout different other pathways such as "Amino sugar and nucleotide sugar metabolism, " "MAPK signaling pathway-plant, " "Plant hormone signal transduction, " "Plant-pathogen interaction, " etc. (Supplementary Table 8). Venn diagrams from KEGG gene ontology were used to visualize which tissues have similar gene expression and which are the most differentiated. Comparison of the number of activated and suppressed transcripts in KEGG categories showed the greatest similarity between early and young pitchers (Figures 3A,B).
Comparative analysis of KEGG terms enriched (by ≥2 times upregulated transcripts, Supplementary Table 6) in each transcriptome revealed that in leaves (vs. early pitcher), they were mostly associated with aging and primary and secondary metabolism, whereas in early and young pitcher (vs. leaves)-with genetic information processing, cell growth and death, development and regeneration, immune system, and  Table 9). Late pitcher was similar to leaves in the downregulation of genetic information processing, cell growth and development, immune system, and signal transduction, and to early/young pitcher-in the suppression of genes assigned to primary and secondary metabolism, and activation of transport and catabolism, immune system, and signaling (Supplementary Table 9 and Figure 3). Compared to the other transcriptomes, late pitcher showed transcriptional induction in the categories of "Digestive system" and "Biosynthesis of other secondary metabolites, " and reduction in the "Endocrine system" and "Drug resistance" (Supplementary Table 9 FIGURE 3 | Comparison of differential expression of leaf and pitcher transcripts in Kyoto encyclopedia of genes and genomes (KEGG) terms. Number of transcripts upregulated (A) and downregulated (B) in early pitcher (PP), young pitcher (YP), and late pitcher (MP) vs. leaves (L) (Supplementary Table 8). Heatmap (C) based on relative expression of transcripts associated with KEGG terms (Supplementary Table 9) in leaves vs. early pitcher, and in PP, YP, and MP vs. leaves. and Figure 3); the latter is possibly correlated with the downregulation of "Defense response" in GO enrichment analysis (Figure 2).

Plant Transcription Factor Database Annotation
Nepenthes transcripts encoding TFs were additionally characterized using plant-specific PlantTFDB. In total, 3,298 TF-encoding transcripts were classified into 48 families, among which the most represented (≥20 members) were MYB/MYBrelated, bHLH, AP2/ERF, C2H2, C3H, bZIP, NAC, WRKY, G2-like, HD-ZIP, B3/ARF, GRAS, Dof, and MADS-domain protein families (Supplementary Table 10). Volcano plots were used to visualize the differences in expression levels between the transcripts, including TF-coding transcripts (Figure 4). Analysis of TF differential expression revealed similarity of the regulatory machinery between leaves and late pitcher and between early and young pitchers (Figure 4 and Supplementary Table 10).

MapMan Annotation
Next, we used MapMan to facilitate comparison of upregulated transcripts (Supplementary Table 6) in early, young, and late pitchers vs. leaves, late and young pitchers vs. early pitcher, and late pitcher vs. young pitchers. Secondary metabolism patterns were similar between early pitcher and young pitchers and between leaves and late pitcher. Compared to leaves and late pitcher, early pitcher and young pitchers were found to have several secondary metabolism pathways upregulated, including biosynthesis of phenylpropanoids, anthocyanins, and alkaloids ( Figure 5 and Supplementary Figure 3). Compared to leaves, photosynthesis was sharply reduced in pitchers at the early developmental stages (early and young pitchers) ( Figure 5) but not at the late stage (mature pitcher) (Supplementary Figure 4). Genes related to biotic and abiotic stress responses were upregulated in pitchers at all developmental stages compared to leaves ( Figure 5 and Supplementary Figure 5).
Overall, these results indicate a loss of leaf identity, i.e., photosynthetic activity, in the pitcher, and the upregulation of defense-related and metabolic processes during its development.

Differentially Expressed Genes in Pitcher Initiation, Growth, and Maturation
Next, we selected DEGs upregulated >4 times and annotated them, using Arabidopsis thaliana UniProt data for the pairwise comparisons of early, young or late pitchers with leaves, and leaves with early pitcher (Supplementary Table 11). The results revealed significant upregulation of transcripts related to flowering initiation and flower development, biosynthesis of secondary metabolites such as anthocyanins, terpenes, jasmonate, and wax, as well as stomatal development and movement, water transport, immunity, and stress response ( Table 1).
In accordance with pitcher initiation on the tendril apical meristem, early pitcher differed from leaves in the activation of processes associated with the initiation and organization of embryonic SAM, and determination of organ abaxial-adaxial polarity ( Table 1 and Supplementary Table 11).
Initiation of the pitcher primordium was accompanied by transcriptional activation of genes associated with the redox state, abiotic stress, respiration, and signaling as well as TFs and R-genes ( Figure 5 and Supplementary Figures 3-5). Early pitcher was similar to young pitcher but both differed from late pitcher in the expression pattern of genes regulating the redox state and induction of signaling, respiration, and TF activity (especially ethylene-responsive factors) (Supplementary Figure 5 and  Furthermore, initiated and developing pitchers showed upregulation of stress-response factors encoding pathogenesisrelated proteins (PR-1C, Fra a 1.03, PRHA, STH-2, STH-21, PTI5, PRHA, SNI1, and thaumatin-like), LRR receptor-like serine/threonine protein kinases, digestive enzymes (peroxidases, glucanases, alfa-and beta-amylases, phospholipases, chitinases, and aspartic proteinases nepenthesin-1 and -2), ABC transporters, dehydration-responsive elements, heat shock proteins, and jasmonate signaling-related molecules (Supplementary Tables 2, 6, 11).
Among the TFs upregulated in early pitcher vs. leaves by more than two times there were multiple homologs of many floweringrelated MADS-domain TFs (Table 1 and Supplementary  Tables 10, 11). We performed phylogenetic analysis of the found 45 proteins in comparison with the known MADS-TFs of a model plant A. thaliana and a representative of Caryophyllales, B. vulgaris; in addition, each of the 45 proteins was separately analyzed with the most similar sequences available from the NCBI GenBank. As a result, the NveMADS1-45 genes were classified as homologs of AGL24, AGL42, SVP, SOC1, AP1/FUL, SEP1/SEP2, SEP3, SEP4, AG, STK, AP3, AGL6, AGL62, AGL65, AGL104, AGL80, and AGL82 (Supplementary Table 12 and Supplementary Figure 6).
Heatmap shows activation of most MADS-box genes during pitcher initiation and development ( Figure 6A). In early
Regarding other important pigments, carotenoids, it should be mentioned that very few genes of the carotenoid pathway, associated with fruit-specific sequestration of pigment in chromoplasts and abscisic acid synthesis, rather than photosynthesis and photoprotection, were found upregulated: chromoplast-specific carotenoid-associated protein CHRC (in leaves and late pitcher vs. early pitcher), 9-cis-epoxycarotenoid dioxygenase NCED1 (in leaves vs. late pitcher and in early pitcher vs. leaves), and carotenoid-cleavage dioxygenase NCED4 (in late pitcher vs. leaves and early pitcher) (Supplementary Table 11).

Validation of RNA-Seq Data by Real-Time Quantitative PCR
Thirteen DEGs were selected for RT-qPCR validation of the transcriptomic data with specific primers (Supplementary Table 16): 10 genes of the anthocyanin biosynthetic pathway (regulatory MYB113 and bHLH001 and structural CHS1, CHS2, CHI, F3H, F3'5'H, DFR, ANS, and UFGT) and three TF genes associated with flowering time and floral organ identity (MYB17 and MADS-box genes AG and SEP3). The results showed significant upregulation of all analyzed genes in early pitcher compared to leaves, with the exception of MYB113 activated in late pitcher (Supplementary Figure 7), which was consistent with the transcriptome data (Figure 6).

DISCUSSION
The Nepenthes pitcher, derived from the tip of the leaf tendril, is thought to emerge from a spontaneously occurred epiascidiate leaf folded into a tubular structure (outer adaxial side inward) with fused margins, which provided selective advantage through improved storage of water and nutrients (Givnish, 2015). To survive in nitrogen-deficient conditions, Nepenthes (like other known carnivorous species with similar trap structures) turned the tubular leaf into a trap-pitcher with prey attraction and digestion abilities (Givnish, 2015).

Stress-Response Proteins Trigger Trap Initiation
Comparative genomic and transcriptomic studies of true leaves and leaf traps in different carnivorous species, including Nepenthes, have shown that signaling pathways involved in prey catching and digestion are similar to those providing protection against pathogens in non-predatory plants (Renner and Specht, 2013;Fukushima et al., 2017). For example, jasmonates associated with response to stress, microbial pathogens, and pests in non-carnivorous plants, in carnivores can perceive signals from captured prey and participate in the development of digestive cavities loaded with hydrolytic enzymes (Pavlovic and Mithöfer, 2019), including pathogenesisrelated proteins associated with stress response, which acquired digestive properties (Pavlovic and Mithöfer, 2019). These data suggest that carnivorous plants evolved using pest defense mechanisms (Fukushima et al., 2017;Pavlovic and Mithöfer, 2019). The results of this study revealed that genes related to stress response and pathogen/pest attacks were upregulated in the early pitcher transcriptome (Supplementary  Tables 6, 11, Figure 5, and Supplementary Figures 3-5), which is consistent with the hypothesis of defense-to-carnivory pitcher evolution.

Pitcher Primordium Organization
In response to defense-related signals, the central vein of the leaf lengthens, producing a tendril with an apical meristem at the tip, which forms the pitcher in the process directed by phyllotaxis signals that promote cell division activity and prevent premature cell differentiation. We found that transcription of the WOX genes was upregulated in pitchers compared to leaves ( Table 1, Supplementary Table 11, and Figure 4), which may be related to lateral organ formation (WOX1 and WOX3), embryo patterning (WOX8 and WOX9), floral transition (WOX13), and determination of vascular stem cell niche (WOX4) (van der Graaff et al., 2009;Zhou et al., 2015).
Consistent with cell division patterns in adaxial tissues shown in carnivorous Sarracenia purpurea (Fukushima et al., 2015), in N. × ventrata early pitcher the genes regulating the abaxialadaxial polarity and initiation/organization of embryonic SAM were activated ( Table 1, Supplementary Table 11, and Figure 4), including abaxial-adaxial polarity genes shown to be upregulated in pitcher-bearing shoots of carnivorous Cephalotus (Fukushima et al., 2017), as well as AS1 and REV genes activated at the site of N. khasiana trap initiation and formation (Dkhar and Pareek, 2019). KAN and PHB genes upregulated in the N. × ventrata pitcher may determine its initial asymmetric growth and YABBY genes may specify the fate of abaxial cells, thus contributing to the proper abaxial/adaxial organization of the trap as shown in Arabidopsis (Eshed et al., 2004). Also similar to Arabidopsis (Inagaki et al., 2006;Rast and Simon, 2012;Shpak, 2013;Amanda et al., 2016), the longitudinal growth of N. × ventrata trap primordium may depend on the mutual activity of the upregulated ER, ERL1, and ERL2 genes, whereas cell division and differentiation (including epidermal cells) may rely on TEBICHI and DEK1, and symmetry-on AS2. GRF genes are known to be involved in the control of cell expansion in the leaf (Dkhar and Pareek, 2014) and, hence, may be important for N. × ventrata trap outgrowth at the tip of the leaf.
Folding of the tubular trap may require SRL2, which regulate leaf rolling through abaxial cell differentiation (Liu et al., 2016), as well as bHLH30 involved in upward leaf curling  and FIL and YABBY2 related to the twisted leaf phenotype (Bowman, 2000).
Overall, the activation and cooperative activity of these genes in N. × ventrata early pitcher may contribute to the epiascidiate leaf folding.

Prey Capture
The pitcher of N. alata (an N. × ventrata parent) catches the prey with the help of anisotropic slippery, highly wet peristome (with water and/or nectar films), waxy inner walls, and release of attractants and sugar-rich nectar (Bohn and Federle, 2004;Bauer et al., 2008). The captured prey is then digested in the lower zone fluid enriched with hydrolytic enzymes (Bauer et al., 2008). Consistent with these data, in N. × ventrata pitcher transcriptomes we observed activation of genes involved in sucrose, wax/cutin, and alkaloid biosynthesis, as well as those encoding proteases and oligosaccharide hydrolases (Supplementary Table 11 and Figure 5).
To attract insects with well-developed CO 2 receptors, closed Nepenthes pitchers accumulate CO 2 emitted after opening, which decreases pitcher photosynthetic capacity, while promoting its growth and respiration and increasing carbohydrate synthesis, cuticular wax density, and humidity levels (Baby et al., 2017). Moreover, CO 2 dissolved in the digestive fluid maintains the optimum pH for the activity of hydrolytic enzymes and nutrient absorbance (Baby et al., 2017). Because of the high CO 2 concentrations, "modified stomata" evolved on the inner side wax layers of Nepenthes pitchers have altered morphology (Baby et al., 2017). The transcriptome data obtained in this study confirmed the role of CO 2 in the development and functioning of the N. × ventrata pitcher, as evidenced by the upregulation of genes related to reversible CO 2 hydration and fixation, CO 2 and water transport during transpiration, and stomatal and guard cell development and function (Supplementary Table 11).
Traps produced by Nepenthes plants have low or no photosynthetic activity and, consequently, lower content of chlorophylls and carotenoids but higher content of anthocyanins (Pavlovič and Saganová, 2015). Accordingly, our comparative transcriptome analysis revealed downregulation of photosynthesis-related genes and upregulation of anthocyanin pathway genes in N. × ventrata pitchers at all stages compared to leaves (Supplementary Table 11 and Figures 5, 6B). Another recent study of Nepenthes pitchers (Goh et al., 2020) also showed the downregulation of genes associated with photosynthesis and increased biosynthesis of secondary metabolites in a pitcher depleted of fluid proteins.

Anthocyanin Pathway in Prey Attraction
Carnivorous plants, including Nepenthes spp., use insects both as pollinators and as food, participating in pollinatorprey conflict (Jürgens et al., 2012). Nepenthes pitchers, which do not or weakly reflect UV radiation visible to insects (Schaefer and Ruxton, 2008), acquire the same red color in pitchers and flower tepals by increasing anthocyanin production (Scharmann et al., 2019).
Red coloration could have initially developed in carnivorous plants as an adaptive trait, since anthocyanin accumulation is often associated with stress responses or nutritional deficiency in plants (Schaefer and Ruxton, 2008). At the same time, it increased prey capture efficiency of the traps by providing attractive visual signals. Capture rates are positively correlated with levels of red pigmentation, probably because insects can detect differences in red light intensity compared to the green background (Schaefer and Ruxton, 2008).
We found that in N. × ventrata red-burgundy pitchers, most of the anthocyanin pathway genes identified in Arabidopsis were expressed, except for F3'H, whose transcripts were not detected (Supplementary Tables 6, 11, Figures 5, 6B, and Supplementary Figure 6). F3'H encodes flavonoid 3'-hydroxylase catalyzing conversion of dihydrokaempferol to dihydroquercetin and its absence may suggest that N. × ventrata lacks the branch of cyanidin synthesis. Products of the other two branches are blue delphinidins and orange-red pelargonidins (Tanaka and Ohmiya, 2008) and, considering N. × ventrata pitcher color, we suggest that biosynthesis of pelargonidin may be the most pronounced branch of the anthocyanin pathway in N. × ventrata pitchers ( Figure 7C).
In some families of core Caryophyllales, flowers and fruits produce betalains: red-violet betacyanins and yellow betaxanthins, which are structurally and biochemically unrelated to anthocyanins and exclude their synthesis (Timoneda et al., 2019). Nepenthaceae are noncore Caryophyllales that produce anthocyanins (Timoneda et al., 2019); nevertheless, the gene encoding a main enzyme of the betalain pathway, 4,5-DOPAextradiol-dioxygenase (DODA), was found to be activated in N. × ventrata, late pitcher (Supplementary Table 6). DODA catalyzes the conversion of 3,4-dihydroxy-L-phenyalanine (L-DOPA) to 4,5-seco-DOPA, which spontaneously forms betalamic acid, the precursor of all betalain compounds. Spontaneous conjugation of betalamic acid with cyclo-DOPA produced by cyclization of oxidized L-DOPA generates redviolet betacyanins (Polturak and Aharoni, 2018). Along with DODA activation, transcription of anthocyanin biosynthesis genes was decreased in late pitcher compared to young pitcher (Figure 6B), which confirms mutual exclusion of anthocyanins and betalains (Timoneda et al., 2019) and indicate a possibility for betalain biosynthesis at later developmental stages of noncore Caryophyllales.

MADS-Domain TFs as Possible Regulators of Pitcher Origin and Development
As part of a long-lasting discussion about the origin of the carnivorous syndrome in plants, Charles Darwin, already in the 19th century, pointed out the similarities between traps and reproductive parts of the flower (stamens and pistils), but since that time, there has been no experimental evidence to support this possibility (Pennisi, 2002). In this study, the MADS-domain TFs ( Table 1 and Supplementary Tables 2, 6, 11, 12), which play decisive roles in plant reproductive development (Smaczniak et al., 2012;Theißen et al., 2016) (Supplementary Tables 14, 15), were revealed as the most intriguing set of transcripts upregulated in N. × ventrata pitchers. In higher plants, MADS-domain TFs regulate flowering initiation (AGL24, AGL42, SOC1, and AP1) and repression (SVP) (Dorca-Fornell et al., 2011;Hu et al., 2014;Valentim et al., 2015). We found that AGL24, two of the four AP1/FUL homologous genes, and six of the eight SOC1/AGL42 homologous genes were upregulated, whereas three of seven SVP homologous genes were downregulated in early pitcher vs. leaf ( Figure 6A). We also observed activation of genes from other gene families (FY, FCA, FLK, FD, and CONSTANS-like) known to control flowering time through regulation of MADS-box gene transcription (Table 1 and Supplementary Tables 6, 11). In N. × ventrata vegetative SAM, AGL24, AGL42, and SOC1 encode proteins which may induce the expression of inflorescence meristem identity gene LEAFY (LFY) that in turn may activate floral meristem identity gene AP1 as shown for Arabidopsis (Smaczniak et al., 2012). We found all these genes, with the exception of LFY, both in the leaves and pitchers (Table 1 and  Supplementary Table 12), suggesting that the expression of these genes is independent of the flowering process and that the single leaf-tendril-trap unit may represent a novel structure (adapted for plant nutrition), as has been shown for grapevine tendrils adapted for plant climbing (Calonje et al., 2004). Transcripts of all four genes homologous to AP1/FUL were found in all types of transcriptomes (Supplementary Table 12), which supports the idea that the leaf, tendril, and trap represent a single structure. In this case, LFY could be transcribed in the meristem prior to the initiation of leaf-tendril-trap differentiation.
Duplication, diversification, and neo-functionalization of the MADS-box genes are thought to underlie the origin of flower organs through leaf modification (Melzer et al., 2010;Theißen et al., 2016). As a result, highly conserved MADSdomain TFs in various combinations represent a mechanism controlling flower meristem differentiation (Theißen et al., 2016). The overwhelming variety of existing flower forms indicates the ongoing evolution of the flower MADS-box set through selective changes in protein-protein interactions, downstream target genes, and regulatory patterns (Bartlett, 2017).
Besides traps, dioecious Nepenthes species form female and male flowers on separate plants and, therefore, must have a complete set of highly conserved MADS-box genes associated with flowering (Subramanyam and Narayana, 1971). Indeed, the N. khasiana male inflorescence transcriptome (Scharmann et al., 2019) contained the similar set of the MADS-box genes as N. × ventrata leaves and pitchers (Supplementary Tables 2, 12).
It is worth noting the mimicry of pitchers as flowers in the patterns of scent emission (Di Giusto et al., 2010) and the similarity of the pitcher with a petal in color and stomatal patterning and that of the pistil with a tubular structure, extrafloral nectaries in the peristome, and lower glands at the base inside the trap, which produce digestive fluid (Gaume and Forterre, 2007;Bauer et al., 2008). In many plant species, the surface of the carpel stigma is often uneven, tuberous, and wet because of the wax cover and sugar-rich sticky exudate, which contribute to more efficient pollen adhesion (Lau et al., 2017). The peristome also has an anisotropic surface coated with a film of water and/or attractive and slippery sugar-rich nectar for prey capture (Bohn and Federle, 2004).
It is postulated that the carpels, as well as other floral organs, emerged as a result of leaf modifications, which is confirmed by the inter-conversion between leaves and flower organs when the A, B, and C classes or E-class of MADS-box genes are overexpressed or inactivated (Scutt et al., 2006).
Besides flower specification, flowering-related MADS-box genes are involved in plant responses to various abiotic stresses and wounding (e.g., from insect feeding) (Castelán-Muñoz et al., 2019), which is consistent with carnivory origin through plant defense mechanisms (Pavlovic and Mithöfer, 2019). Moreover, some studies indicate that MADS-box TFs may positively affect biosynthesis of anthocyanins (Jaakola et al., 2010;Zhao et al., 2019) used by traps to attract insects. The presence of highly conserved MADS-box genes associated with flowering and defense in all angiosperms may be the reason why a pitfall trap originated independently six times in diverse plant lineages.
Thus, non-carnivorous Nepenthes predecessor plants could use the MADS-box gene set specifying floral meristems and organs for adaptation to stressful conditions and nutrient deficiency. As a result, they develop a new structure, leaf-tendriltrap, based on the vegetative leaf, using the existing pathways for the formation of floral organs (tepal and pistil) (Figure 7). A recently shown WGD event in the last common ancestor of Droseraceae carnivorous plants (Palfalvi et al., 2020) led to the MADS-box genes duplication, and it can be speculated that leaves have acquired the ability to express such MADS-box paralogous genes, which, among other mechanisms, resulted in the leaftendril-trap development.
It has been shown that MADS-box genes involved in the initiation of flowering and floral organogenesis can also play an important role in the development of leaves (Burko et al., 2013;Gregis et al., 2013) and roots (Gan et al., 2005;Alvarez-Buylla et al., 2019). Therefore, the leaf-tendril-trap structure can be considered a specialized modified organ with its own unique regulatory pathway that evolved through co-option of genes from the networks controlling SAM identity and organization, leaf and root development, and flower morphogenesis.
The results obtained may clarify the genetic patterns of pitcher trap initiation and development and help answer the question of the origin of the plant carnivory syndrome. For example, given that a protocol for Nepenthes mirabilis in vitro regeneration and Agrobacterium-mediated transformation has been developed (Miguel et al., 2020), it is possible to obtain transgenic Nepenthes plants in which individual MADS-box genes are either overexpressed or silenced through various mechanisms such as CRISPR/Cas9-based genome editing. It would also be useful to compare the flower-and trapspecific paralogs of the MADS-box genes, including their regulatory regions. Evaluation of the morphology, genomics, and proteomics of traps from such transgenic plants could shed light on the trap-specific functions of the analyzed MADS-box genes and their roles in trap evolution. Also, modern methods for analysis of protein-protein interactions and transcription factor target genes should make possible to compare the functional activity of MADS-box genes in the trap with that of MADSbox genes known to be involved in flowering initiation and flower development.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: NCBI, accession number: PRJNA487526 (https://www.ncbi.nlm. nih.gov/bioproject/PRJNA487526).

AUTHOR CONTRIBUTIONS
NR and EK: conceptualization. MF, EK, and EG: plant material. NR, AB, and AM: methodology. AB: software. MF and MS: validation. AB and AS: formal analysis. NR: data curation. AS: Writing-original draft preparation. AS, EK, and NR: writingreview and editing. NR: supervision. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the Russian Science Foundation  and by the Ministry of Science and Higher Education of the Russian Federation.

ACKNOWLEDGMENTS
We are sincerely grateful to the late Konstantin Skryabin, who initiated this project and provided the opportunity to explore unusual unique plants at our Institute. We would like to thank Marina Chuenkova for English language editing. This work was performed using the experimental climate control facility in the Institute of Bioengineering (Research Center of Biotechnology, Russian Academy of Sciences).