Integrated Omics Analyses Identify Key Pathways Involved in Petiole Rigidity Formation in Sacred Lotus

Sacred lotus (Nelumbo nucifera Gaertn.) is a relic aquatic plant with two types of leaves, which have distinct rigidity of petioles. Here we assess the difference from anatomic structure to the expression of genes and proteins in two petioles types, and identify key pathways involved in petiole rigidity formation in sacred lotus. Anatomically, great variation between the petioles of floating and vertical leaves were observed. The number of collenchyma cells and thickness of xylem vessel cell wall was higher in the initial vertical leaves’ petiole (IVP) compared to the initial floating leaves’ petiole (IFP). Among quantified transcripts and proteins, 1021 and 401 transcripts presented 2-fold expression increment (named DEGs, genes differentially expressed between IFP and IVP) in IFP and IVP, 421 and 483 proteins exhibited 1.5-fold expression increment (named DEPs, proteins differentially expressed between IFP and IVP) in IFP and IVP, respectively. Gene function and pathway enrichment analysis displayed that DEGs and DEPs were significantly enriched in cell wall biosynthesis and lignin biosynthesis. In consistent with genes and proteins expressions in lignin biosynthesis, the contents of lignin monomers precursors were significantly different in IFP and IVP. These results enable us to understand lotus petioles rigidity formation better and provide valuable candidate genes information on further investigation.


Introduction
Plant architecture is the embodiment of space utilization. There are numerous factors contribute to plant architecture including stem height, leaf and branching mode [1]. In crops, several factors such as stem rigidity, branching pattern, and plant height which strongly influences crop yields and efficiency of harvesting have been extensively studied [2,3]. The sacred lotus is a relic aquatic plant with two leaf types, namely the floating leaf and the vertical leaf. Floating leaf petiole exhibits flexibility and floating tendency despite it being longer than the depth of the water, while the vertical leaf petiole is rigid and erect ( Figure 1). Moreover, there is no reported interconvertibility between these two kinds of leaves during the entire life span of lotus, meaning that lotus is not a heterophyllous plant. Stem or petiole rigidity is species-specific and plastic [4] controlled by environmental factors including light, temperature, water and nutrients, among others.
Cell wall makes plant cell different from animal cell. Cell wall contributes to plant architecture, organ development, transport of water and nutriment and defense. In order to accomplish variety function mentioned above, cell in plant was specialized into collenchyma, sclerenchyma and xylem [5][6][7]. Generally, cell wall can be divided into middle lamella, primary wall and secondary wall if the cell wall has secondary thickening. The main components in secondary cell wall are polysaccharide and lignin [8]. Polysaccharide mainly includes cellulose, hemicellulose and pectin. Lignin is compounded of three primary monolignols including coniferyl alcohol, sinapyl alcohol and p-coumaryl alcohol in varieties of proportions [9]. The components of cell wall are diversified depending on the plant taxa, age, tissue and cell type [8]. Studies have verified that many genes are associated with cell wall biosynthesis and secondary cell wall thickening in Arabidopsis [10], Populus [6] and other species [11]. Typically, cellulose synthase [12,13], glycosyltransferases and glycosyl hydrolases [8,14] play a crucial role in polysaccharide biosynthesis, and enzymes involved in phenylalanine deamination, hydroxylation, methylation and redox reactions in phenylpropanoid pathway play a key role in lignin biosynthesis [15]. Transcription factor families are engaged more in the upstream of cell wall biosynthesis and modification, including NAC domain (VND6/7:vascular-related NAC domain 6/7, SND1: secondary wall-associated NAC domain 1 and NST1: NAC secondary wall thickening promoting factor) and HD-ZIP homeobox, which were reported that they play sufficient role in influencing secondary cell wall synthesis in xylem vessel [16][17][18]. Members in MYB family are also informed as regulators involved in secondary cell wall formation [19], while members in E2F family were affirmed as a fundamental upstream transcriptional regulator of VND6/7 and other secondary cell wall biosynthesis genes [10].
Despite the increasing understanding on cell wall biosynthesis in the past two decades, a lot still remain elusive regarding related gene regulation. The findings in heterophyllous plants recommend that light and temperature are involved in transformation of heterophyllous leaves [20], and internal factors including phytohormones, such as gibberellin (GA), abscisic acid (ABA) and ethylene play a dominant role in the induction of heterophylly in plants [21]. Besides, some studies focused on Arabidopsis petioles indicating light quality, phytohormones and crosstalk between light and phytohormones mediate petiole elongation [22], and it is ascertained that numerous genes are involved in cell wall formation and secondary cell wall thickening in Arabidopsis [23], rice [24] and Populus [25]. Since there is little research done, it is unclear whether lotus has the same mechanism of gene network regulation as typical heterophyllous plants. Furthermore, lotus is the only aquatic plant which has a combination of floating and vertical leaves, and it is a good example to probe the mechanism that controls differentiation of these two leaves. These findings will facilitate us to understand better the mechanism of cell wall formation.
High-throughput profiling techniques used in transcriptome and proteome are effective in exploring complex biologic processes at the overall level [26]. The expression levels of mRNAs and proteins can be accurately measured using RNA-sequencing, iTRAQ (isobaric tags for relative and absolute quantitation) and TMT (tandem mass tag) technologies. As we known, the regulation of genes at transcriptional, translational and post-translational levels are different, integrated analysis of these data are necessary [27]. In Arabidopsis, transcriptome and proteome technologies were used to identified genes involved in cell wall biosynthesis. It also indicates no clear correlation was found between the omics data demonstrating complementary in different methods [28,29].
In this study, the IFP (initial floating leaves' petiole) and IVP (initial vertical leaves' petiole) were used to as the research material for analysis through RNA sequencing and protein labeling quantification with tandem mass tags (TMT) technology to explore the mechanism that controls differentiation of these two types of leaves in lotus. By analyzing DEGs (genes differentially expressed between IFP and IVP), DEPs (proteins differentially expressed between IFP and IVP) and the abundance of metabolite in the IFP and IVP, we found that genes highly expressed in IVP in several central biologic processes, such as cell wall biosynthesis, organization, assembly and lignin biosynthesis were associated with lotus petiole rigidity. The integrated analysis of transcriptomic and proteomic data provides a clue to understand the formation of different petioles in lotus.

Phenotypic Evaluation of Lotus Petioles at IFP and IVP Stages
Lotus, as an herbaceous perennial plant, sprouts as temperature rises toward spring. At the beginning, several initial leaves (three to five leaves) grow and float on, but not rise above the water surface in spite of long enough petioles. Then the succeeding leaves stand out of the water (Figure 1a,b, Figure S1). Generally speaking, initial leaf floats on the water surface all the time. Vertical leaf rises above the water surface once the petiole length is larger than water depth (Figure 1c, Figure S1). Moreover, both the folded and unfolded vertical leaves can rise above the water surface (Figure 1d, Figure S1). We have not observed any interconvertibility between IFP and IVP under green house and field conditions. Numbers of lateral branches grow continuously during its growing season, and each branch develops the floating leaf first, followed by the vertical leaf.

Phenotypic Evaluation of Lotus Petioles at IFP and IVP Stages
Lotus, as an herbaceous perennial plant, sprouts as temperature rises toward spring. At the beginning, several initial leaves (three to five leaves) grow and float on, but not rise above the water surface in spite of long enough petioles. Then the succeeding leaves stand out of the water ( Figure  1a,b, Figure S1). Generally speaking, initial leaf floats on the water surface all the time. Vertical leaf rises above the water surface once the petiole length is larger than water depth (Figure 1c, Figure S1). Moreover, both the folded and unfolded vertical leaves can rise above the water surface (Figure 1d, Figure S1). We have not observed any interconvertibility between IFP and IVP under green house and field conditions. Numbers of lateral branches grow continuously during its growing season, and each branch develops the floating leaf first, followed by the vertical leaf. and arrow 4 indicates mature vertical leaves' petiole (MVP), respectively; (e) measurements of breaking resistances in two type petioles were performed using a prostrate tester; (g,h) show status of petioles before and after putting pressure on it; (f) shows the length of petioles which were cultured in pot with 20 cm depth water. Data are means ± SD from 9 independent petioles (IFP and IVP). Asterisks indicate significant changes according to Student's t-test (* p < 0.05).
For the purpose of quantifying the phenotype difference between the floating and vertical leaves, the breaking resistance of two type petioles was measured. The breaking resistance of IVP was significantly higher than that in IFP ( Figure 1e). Furthermore, we observed two types of petioles in lotus with low water depth. The minimum length of IFP is over water depth at that condition, and IVP was greater than IFP. Even though the petiole length of floating leaf was greater than water depth as it matures, it still could not stand erect on water (Figure 1f). This result suggests the breaking indicates mature vertical leaves' petiole (MVP), respectively; (e) measurements of breaking resistances in two type petioles were performed using a prostrate tester; (g,h) show status of petioles before and after putting pressure on it; (f) shows the length of petioles which were cultured in pot with 20 cm depth water. Data are means ± SD from 9 independent petioles (IFP and IVP). Asterisks indicate significant changes according to Student's t-test (* p < 0.05).
For the purpose of quantifying the phenotype difference between the floating and vertical leaves, the breaking resistance of two type petioles was measured. The breaking resistance of IVP was significantly higher than that in IFP ( Figure 1e). Furthermore, we observed two types of petioles in lotus with low water depth. The minimum length of IFP is over water depth at that condition, and IVP was greater than IFP. Even though the petiole length of floating leaf was greater than water depth as it matures, it still could not stand erect on water (Figure 1f). This result suggests the breaking resistance, but not the petiole length being the main factor affecting the leaf emergence from water surface.

Anatomic Structure Analysis of Petioles at IFP and IVP Stages
To explore the structure difference of petioles in floating and vertical leaf, we studied the sections of IFP and IVP. There were four primarily air cavities (the numbers of air cavities will increase as the petioles grow) in both transection of IFP and IVP. Numbers of small air cavities were scattered around the chief air cavities (Figure 2a,b), and some of them could develop into the chief air cavities along with petioles growth. The vascular bundles, which increase cell wall resistance and rigidity, were scattered around the big air cavities as well. These data showed the number of vascular bundles in IVP is more than those in IFP (Figure 2a,b,e). Furthermore, the epidermal cell had high level of cutinization in IVP than that in IFP. The outer cortex cells in IVP have specialized into collenchyma by uneven thickening of cell wall to strength the rigidity of petioles. In IFP, by contrast, the cortex cells had no cell wall thickening. (Figure 2c,d). These vascular bundles and collenchyma could contribute to cell wall resistance and rigidity thus supporting petioles emergence from water in IVP. resistance, but not the petiole length being the main factor affecting the leaf emergence from water surface.

Anatomic Structure Analysis of Petioles at IFP and IVP Stages
To explore the structure difference of petioles in floating and vertical leaf, we studied the sections of IFP and IVP. There were four primarily air cavities (the numbers of air cavities will increase as the petioles grow) in both transection of IFP and IVP. Numbers of small air cavities were scattered around the chief air cavities (Figure 2a,b), and some of them could develop into the chief air cavities along with petioles growth. The vascular bundles, which increase cell wall resistance and rigidity, were scattered around the big air cavities as well. These data showed the number of vascular bundles in IVP is more than those in IFP (Figure 2a,b,e). Furthermore, the epidermal cell had high level of cutinization in IVP than that in IFP. The outer cortex cells in IVP have specialized into collenchyma by uneven thickening of cell wall to strength the rigidity of petioles. In IFP, by contrast, the cortex cells had no cell wall thickening. (Figure 2c,d). These vascular bundles and collenchyma could contribute to cell wall resistance and rigidity thus supporting petioles emergence from water in IVP. (e) statistics of vascular bundles number in IFP and IVP, data are means ± SD from 10 independent transverse sections (IFP and IVP); (f) determination of crude cellulose and (g) lignin. Data are means ± SD from 3 independent biologic repeats. Asterisks and double asterisks indicate significant changes compared to the control as assessed according to Student's t-test (* p < 0.05 and ** p < 0.01).
To further study the structure character of petiole, thickness of xylem vessels cell wall was measured (Figure 3a-d). It showed that xylem vessels cell wall in IFP was 0.4 μm compared with IVP (0.8 μm) (Figure 3g). This result suggests that xylem vessel in IVP may provide more strength to support IVP to remain upright. The main components of cell wall are polysaccharide and lignin, and the content of cellulose and lignin in cell wall were measured. This result noted that the cellulose (e) statistics of vascular bundles number in IFP and IVP, data are means ± SD from 10 independent transverse sections (IFP and IVP); (f) determination of crude cellulose and (g) lignin. Data are means ± SD from 3 independent biologic repeats. Asterisks and double asterisks indicate significant changes compared to the control as assessed according to Student's t-test (* p < 0.05 and ** p < 0.01).
To further study the structure character of petiole, thickness of xylem vessels cell wall was measured (Figure 3a-d). It showed that xylem vessels cell wall in IFP was 0.4 µm compared with IVP (0.8 µm) (Figure 3g). This result suggests that xylem vessel in IVP may provide more strength to support IVP to remain upright. The main components of cell wall are polysaccharide and lignin, and the content of cellulose and lignin in cell wall were measured. This result noted that the cellulose contents were 18.7% and 24.8% in IFP and IVP, and the lignin contents were 13.3% and 21.9% in IFP and IVP, respectively (Figure 2f,g). Both cellulose and lignin contents were significantly higher in IVP than that in IFP. This result is consistent with morphologic results which displayed more vascular bundles, collenchyma and cell wall thickening in IVP compared to IFP.  (Figure 2f,g). Both cellulose and lignin contents were significantly higher in IVP than that in IFP. This result is consistent with morphologic results which displayed more vascular bundles, collenchyma and cell wall thickening in IVP compared to IFP.

Overview of Transcriptomic Analysis and Transcriptomic Data Validation
To explore the different genes involved in cell wall formation in IFP and IVP, comparative transcriptomic experiment was done. RNA-seq was performed for IFP and IVP groups on the Illumina platform and each group included three biologic replicates. Generally, 87%-90% of the clean reads were mapped to the genome of the sacred lotus (https://bioinformatics.psb.ugent.be/plaza/versions/plaza_v4_dicots/). The expression estimates were calculated as fragments per kilobase per million reads (FPKM). Pearson's correlation coefficient was used to evaluate repeatability of individual sample ( Figure S2a). A total of 20,109 genes or transcripts (about 75% predicted genes from genome) were obtained from the six samples ( Figure  4a). Differentially expressed genes (DEGs) were screened based on an absolute fold change value of |log2 ratio| ≥ 1 and false discovery rate (FDR) ≤ 0.01. A total of 1422 DEGs were found, including 401 upregulated genes and 1021 downregulated genes (IVP vs. IFP) (Figure 4b, Figure S2b, Table S1).
The expression levels of 14 genes involved in the cell wall biosynthesis and 8 genes with different expression patterns in IFP and IVP were determined by quantitative reverse transcription polymerase chain reaction (qRT-PCR) to validate the transcriptomic data. A total of 19 genes shared the same expression pattern between transcriptomic and qRT-PCR results ( Figure 4c, Table S2). The transcriptomic data and qRT-PCR results of these genes were highly correlated (r = 0.8339; Figure  4d). These results further proved that the transcriptomic data were reliable.

Overview of Transcriptomic Analysis and Transcriptomic Data Validation
To explore the different genes involved in cell wall formation in IFP and IVP, comparative transcriptomic experiment was done. RNA-seq was performed for IFP and IVP groups on the Illumina platform and each group included three biologic replicates. Generally, 87-90% of the clean reads were mapped to the genome of the sacred lotus (https://bioinformatics.psb.ugent.be/plaza/versions/ plaza_v4_dicots/). The expression estimates were calculated as fragments per kilobase per million reads (FPKM). Pearson's correlation coefficient was used to evaluate repeatability of individual sample ( Figure S2a). A total of 20,109 genes or transcripts (about 75% predicted genes from genome) were obtained from the six samples ( Figure 4a). Differentially expressed genes (DEGs) were screened based on an absolute fold change value of |log2 ratio| ≥ 1 and false discovery rate (FDR) ≤ 0.01. A total of 1422 DEGs were found, including 401 upregulated genes and 1021 downregulated genes (IVP vs. IFP) (Figure 4b, Figure S2b, Table S1).
The expression levels of 14 genes involved in the cell wall biosynthesis and 8 genes with different expression patterns in IFP and IVP were determined by quantitative reverse transcription polymerase chain reaction (qRT-PCR) to validate the transcriptomic data. A total of 19 genes shared the same expression pattern between transcriptomic and qRT-PCR results ( Figure 4c, Table S2). The transcriptomic data and qRT-PCR results of these genes were highly correlated (r = 0.8339; Figure 4d). These results further proved that the transcriptomic data were reliable.

Overview of Quantitative Proteomics Analysis and Data Validation Using PRM
To scrutinize the different proteins involved in cell wall formation of IFP and IVP, quantitative and comparative proteomics was done. A total of 25,529 peptides belonging to 5808 proteins (about 22% predicted genes from genome) were identified in the petiole samples. Among these identified proteins, 4855 (about 18% predicted genes from genome) of them were quantified (Figure 4a). Pearson's correlation coefficient was used to evaluate repeatability of individual sample ( Figure S3a). After filtering with an expression fold change of >1.5 and p-value < 0.05, 904 out of the quantified proteins were defined as differentially expressed (DEPs) between IFP and IVP. Among these DEPs, the abundance of 421 proteins in IFP was 1.5 times higher than in IVP, and the abundance of 483 proteins in IVP was 1.5 times higher than in IFP (Table S3).
To validate the reliability of tandem mass tag-labeling (TMT-labeling) protein quantification in this study, parallel reaction monitoring (PRM) method was used. The results showed that 11 out of 15 selected proteins were successfully quantified using PRM methods, and the changes in all 11 proteins abundances between IFP and IVP were consistent with TMT-labeling protein quantification ( Figure S4a, Table S4). The PRM data and proteomics results of these proteins were highly correlated (r = 0.6301, Figure S4b). This result confirms that our protein quantification result was reliable.

Integrated Analysis of Transcriptome and Proteome Data
To analyze the correlation of DEGs and DEPs identified in IFP and IVP, data from transcriptome and proteome were compared. Globally, a total of 20,109 transcripts (approximate 75.36% genome)

Overview of Quantitative Proteomics Analysis and Data Validation Using PRM
To scrutinize the different proteins involved in cell wall formation of IFP and IVP, quantitative and comparative proteomics was done. A total of 25,529 peptides belonging to 5808 proteins (about 22% predicted genes from genome) were identified in the petiole samples. Among these identified proteins, 4855 (about 18% predicted genes from genome) of them were quantified ( Figure 4a). Pearson's correlation coefficient was used to evaluate repeatability of individual sample ( Figure S3a). After filtering with an expression fold change of >1.5 and p-value < 0.05, 904 out of the quantified proteins were defined as differentially expressed (DEPs) between IFP and IVP. Among these DEPs, the abundance of 421 proteins in IFP was 1.5 times higher than in IVP, and the abundance of 483 proteins in IVP was 1.5 times higher than in IFP (Table S3).
To validate the reliability of tandem mass tag-labeling (TMT-labeling) protein quantification in this study, parallel reaction monitoring (PRM) method was used. The results showed that 11 out of 15 selected proteins were successfully quantified using PRM methods, and the changes in all 11 proteins abundances between IFP and IVP were consistent with TMT-labeling protein quantification ( Figure S4a, Table S4). The PRM data and proteomics results of these proteins were highly correlated (r = 0.6301, Figure S4b). This result confirms that our protein quantification result was reliable.

Integrated Analysis of Transcriptome and Proteome Data
To analyze the correlation of DEGs and DEPs identified in IFP and IVP, data from transcriptome and proteome were compared. Globally, a total of 20,109 transcripts (approximate 75.36% genome) were identified in lotus petioles, and transcripts were detected for 93.12% of the proteins (Figure 4a). To explore the relationship between protein abundance and their corresponding gene expression, we conducted a correlation analysis using Spearman's rank correlation coefficient method between the transcriptome and quantitative proteome data. After analyzing, the expression levels of all the transcripts and their corresponding quantified proteins from IVP vs. IFP showed weak correlation (r = 0.2765, Figure 5a). The weak correlation was also observed in various multi-omics researches in Arabidopsis, rice and eggplant [26,30,31]. At the same time, a stronger correlation was discerned between the DEGs and their corresponding DEPs (r = 0.6148, Figure 5d). The expression ratio of proteins and their corresponding mRNAs with the same or opposite trend (both upregulated or both downregulated) were also plotted, and higher positive or negative correlation was indicated (Figure 5b,c,e,f). were identified in lotus petioles, and transcripts were detected for 93.12% of the proteins (Figure 4a). To explore the relationship between protein abundance and their corresponding gene expression, we conducted a correlation analysis using Spearman's rank correlation coefficient method between the transcriptome and quantitative proteome data. After analyzing, the expression levels of all the transcripts and their corresponding quantified proteins from IVP vs. IFP showed weak correlation (r = 0.2765, Figure 5a). The weak correlation was also observed in various multi-omics researches in Arabidopsis, rice and eggplant [26,30,31]. At the same time, a stronger correlation was discerned between the DEGs and their corresponding DEPs (r = 0.6148, Figure 5d). The expression ratio of proteins and their corresponding mRNAs with the same or opposite trend (both upregulated or both downregulated) were also plotted, and higher positive or negative correlation was indicated ( Figure  5b,c,e,f). In the present study, a total of 1422 DEGs and 904 DEPs were found. Both the DEGs and DEPs were classified into various functional groups ( Figure 4b). Among these DEGs and DEPs, 135 of them had quantitative information both in transcript and protein level. For ease of description, the 135 genes both identified in transcriptional and translational levels were named as core genes. MapMan analysis showed the 135 core genes were classified in 22 functional groups. Except for unknown function genes, 72% of core genes were sorted into five larges groups including miscellaneous, stress, signaling, protein and cell wall ( Figure S5a). Among the 135 core genes, 126 of them showed same trend in transcriptional and translational level, while only 9 genes showed the opposite trend at the two levels ( Figure S5b, Tables S1 and S3). This suggests that these core genes play constant role at both RNA and protein levels in petioles development. Besides, a total of 67 transcription factors were identified in DEGs and 3 transcription factors were identified in DEPs. Only AP2/ERF and B3 domain-containing transcription factor RAV1-like both identified in DEGs and DEPs (Table S5). Combining MapMan and GO annotation results, 11 non-redundant genes out of 135 core genes were predicted as being involved in cell wall biosynthesis, degradation and assembly (Table 1). GO and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis of the 135 core genes showed that majority genes were enriched in phenylpropanoid biosynthesis and chitin catabolic process ( Figure S5c,d). This suggested that these genes may play an important and sustained role in lotus petioles formation. However, the few core genes identified both in In the present study, a total of 1422 DEGs and 904 DEPs were found. Both the DEGs and DEPs were classified into various functional groups (Figure 4b). Among these DEGs and DEPs, 135 of them had quantitative information both in transcript and protein level. For ease of description, the 135 genes both identified in transcriptional and translational levels were named as core genes. MapMan analysis showed the 135 core genes were classified in 22 functional groups. Except for unknown function genes, 72% of core genes were sorted into five larges groups including miscellaneous, stress, signaling, protein and cell wall ( Figure S5a). Among the 135 core genes, 126 of them showed same trend in transcriptional and translational level, while only 9 genes showed the opposite trend at the two levels ( Figure S5b, Tables S1 and S3). This suggests that these core genes play constant role at both RNA and protein levels in petioles development. Besides, a total of 67 transcription factors were identified in DEGs and 3 transcription factors were identified in DEPs. Only AP2/ERF and B3 domain-containing transcription factor RAV1-like both identified in DEGs and DEPs (Table S5). Combining MapMan and GO annotation results, 11 non-redundant genes out of 135 core genes were predicted as being involved in cell wall biosynthesis, degradation and assembly (Table 1). GO and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis of the 135 core genes showed that majority genes were enriched in phenylpropanoid biosynthesis and chitin catabolic process ( Figure S5c,d).
This suggested that these genes may play an important and sustained role in lotus petioles formation. However, the few core genes identified both in transcriptional and translational levels may be due to the fact that the genes involved in lotus petioles formation do not express synchronously. As reported by Casas-Vila et al. in Drosophila melanogaster, they found unusual behavior of RNA and protein during embryogenesis [32]. This indicates a strong post-translational regulation and the necessity of joint analysis of the transcriptome and proteome in biology. Superscript a indicates 11 core genes in DEPs (differentially expressed proteins) and DEGs (differentially expressed genes). Null means proteins or transcripts were not found.

Functional Enrichment of Quantified DEGs and DEPs
Excepting the 135 core genes, the function of other DEGs and DEPs was analyzed by MapMan software, GO and KEGG annotation tools. The DEPs were classified into 4 clusters according to the fold changes. An abundance of 164 proteins in IFP was 2-fold higher than IVP (Q1), while the abundance of 257 proteins in IFP was 1.5-2-fold higher than IVP (Q2). In addition, the abundance of 346 proteins in IVP was 1.5-2-fold higher than IFP (Q3), while the abundance of 137 proteins in IVP was 2-fold higher than IFP (Q4) ( Figure S3b). Functional enrichment including GO enrichment, protein domain enrichment, and KEGG pathway enrichment of the DEPs in four clusters was analyzed. Among these 164 proteins in cluster Q1, 13 non-redundant proteins were enriched in cell wall formation and degradation. Among these 137 proteins in cluster Q4, 7 non-redundant proteins were enriched in cell wall biogenesis, cell wall organization and cell wall assembly. In cluster Q4, a total of 9 non-redundant proteins were enriched in polysaccharide including cellulose and hemicellulose biosynthetic process. Nonetheless, in cluster Q2 and Q3, no protein was enriched in cell wall formation and degradation related biologic processes (Figure 6a). This protein functional enrichment analysis coincides with the anatomic results and it indicates that cell wall biosynthesis in IFP and IVP may be different.
Generally, the GO enrichment analysis of these DEGs were shown in ( Figure S6a,b). The DEGs were analyzed in upregulated and downregulated gene clusters. It showed that a total of 219, 141 and 87 upregulated genes in IVP were enriched in anatomic structure development, anatomic structure morphogenesis and cell growth. The principal point is that a total of 92 upregulated genes in IVP were enriched in cell wall organization or biogenesis. However, the upregulated genes in IFP were not enriched in pathways related to cell wall organization or biogenesis ( Figure S2c). In MapMan analysis, a total of 46 non-redundant DEGs and 63 non-redundant DEPs were identified in cell wall biosynthesis and degradation pathway. Moreover, a total of 10 DEGs and 16 DEPs were identified in lignin biosynthesis. More DEPs than DEGs were identified in cell wall biosynthesis and lignin biosynthesis pathway may indicate proteins contribute more to the difference among samples.

Functional Enrichment of Quantified DEGs and DEPs
Excepting the 135 core genes, the function of other DEGs and DEPs was analyzed by MapMan software, GO and KEGG annotation tools. The DEPs were classified into 4 clusters according to the fold changes. An abundance of 164 proteins in IFP was 2-fold higher than IVP (Q1), while the abundance of 257 proteins in IFP was 1.5-2-fold higher than IVP (Q2). In addition, the abundance of 346 proteins in IVP was 1.5-2-fold higher than IFP (Q3), while the abundance of 137 proteins in IVP was 2-fold higher than IFP (Q4) ( Figure S3b). Functional enrichment including GO enrichment, protein domain enrichment, and KEGG pathway enrichment of the DEPs in four clusters was analyzed. Among these 164 proteins in cluster Q1, 13 non-redundant proteins were enriched in cell wall formation and degradation. Among these 137 proteins in cluster Q4, 7 non-redundant proteins were enriched in cell wall biogenesis, cell wall organization and cell wall assembly. In cluster Q4, a total of 9 non-redundant proteins were enriched in polysaccharide including cellulose and hemicellulose biosynthetic process. Nonetheless, in cluster Q2 and Q3, no protein was enriched in cell wall formation and degradation related biologic processes (Figure 6a). This protein functional enrichment analysis coincides with the anatomic results and it indicates that cell wall biosynthesis in IFP and IVP may be different.
Generally, the GO enrichment analysis of these DEGs were shown in ( Figure S6a,b). The DEGs were analyzed in upregulated and downregulated gene clusters. It showed that a total of 219, 141 and 87 upregulated genes in IVP were enriched in anatomic structure development, anatomic structure morphogenesis and cell growth. The principal point is that a total of 92 upregulated genes in IVP were enriched in cell wall organization or biogenesis. However, the upregulated genes in IFP were not enriched in pathways related to cell wall organization or biogenesis ( Figure S2c). In MapMan analysis, a total of 46 non-redundant DEGs and 63 non-redundant DEPs were identified in cell wall biosynthesis and degradation pathway. Moreover, a total of 10 DEGs and 16 DEPs were identified in lignin biosynthesis. More DEPs than DEGs were identified in cell wall biosynthesis and lignin biosynthesis pathway may indicate proteins contribute more to the difference among samples.  Abundance of 164 proteins in IFP was 2-fold higher than IVP (Q1) and abundance of 257 proteins in IFP was 1.5-2-fold higher than IVP (Q2). Also, abundance of 346 proteins in IVP was 1.5-2-fold higher than IFP (Q3) and abundance of 137 proteins in IVP was 2-fold higher than IFP (Q4). Red star indicates cell wall related process; (b) MapMan analysis shows that DEGs and DEPs in different metabolic pathways in lotus petioles. Red squares indicate enhancement and green squares were the opposite. Squares above the blue line represent the transcripts and squares below the blue line represent the proteins.N means not found. Arabic numerals 1 to 7 in red frame represent cell wall proteins, pectin esterases, cellulose synthesis, cell wall degradation, cell wall precursor synthesis, major CHO metabolism and lignin biosynthesis pathways.

Polysaccharide and Lignin Biosynthesis Pathway Were Highly Activated in IVP Compared to IFP
Plant cell walls are composed of polysaccharides, including cellulose, hemicellulose and pectin, lignin, proteins and minerals. Here numerous genes involved in cellulose, hemicellulose, and pectin biosynthesis pathways were identified with a changed expression both at transcriptional and protein level ( Table 1). In the DEGs and DEPs identified in this study, 12 genes (3 transcripts and 9 proteins) involved in cell wall precursor synthesis were identified, and all the 12 genes showed expression increment in IVP. A total of 18 genes (6 transcripts and 12 proteins) involved in cellulose and hemicellulose synthesis were identified, and 16 of them showed expression increment in IVP except two transcripts. Moreover, a total of 26 genes (14 transcripts and 12 proteins) involved in cell wall biosynthesis and modifications were identified, and 15 of them showed expression increment in IVP. Nonetheless, a total of 24 genes (11 transcripts and 13 proteins) involved in cell wall degradation were identified and 13 of them showed expression decline in IVP (Table 1, Figure 6b). It indicated that more genes involved in cell wall precursor synthesis, cellulose and hemicellulose synthesis expressed highly in IVP may play foremost roles in its rigidity formation. Cell wall biosynthesis related genes like cellulose synthases (CESA), UDP-glucose 6-dehydrogenases (UGD), UDP-glucuronic acid decarboxylases (UXS), irregular xylem (IRX), and sucrose synthases (SUSY) showed significant higher expression level in IVP than that in IFP ( Figure S6c).
In the KEGG pathway enrichment analysis, a total of 114 proteins were enriched in 6 different pathways such as starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, glycosphingolipid biosynthesis, phenylpropanoid biosynthesis, flavonoid biosynthesis and biosynthesis of amino acids (Table 2). Remarkably, 14 proteins highly expressed in IVP and 29 proteins highly expressed in IFP were enriched in phenylpropanoid biosynthesis. Furthermore, the abundance of 10 enzymes including phenylalanine ammonia-lyases (PAL, 3 PAL proteins appear in DEPs), cinnamate 4-hydroxylase (C4H, 1 C4H protein appears in DEPs), 4-(hydroxy)cinnamoyl CoA ligase (4CL, 1 4CL protein appears in DEPs), cinnamoyl CoA reductase (CCR, 1 CCR protein appears in DEPs), cinnamyl alcohol dehydrogenase/sinapyl alcohol dehydrogenase (CAD/SAD, 1 CAD protein appears in DEPs), caffeic acid/5-hydroxyferulic acid O-methyltransferase (COMT, 1 COMT protein appears in DEPs), ferulate 5-hydroxylase (F5H, 1 F5H protein appears in DEPs) and hydroxycinnamoyl CoA: shikimate hydroxycinnamoyltransferase/hydroxycinnamoyl CoA: quinate hydroxycinnamoyltransferase (CST/CQT, 1 CST protein appears in DEPs), which are enzymes activated on the upstream of lignin biosynthesis, were higher in IVP than that in IFP. Peroxidase is considered as the last enzyme to catalyze substrate in lignin biosynthesis. In this study, two of peroxidases were highly expressed in IVP compared to IFP and 25 of them were on the contrary. Moreover, the abundance of beta-glucosidase (BGLU, 3 BGLU proteins appear in DEPs) and aldehyde dehydrogenase (REF1, 1 REF protein appears in DEPs), which consume substrates without lignin output in lignin biosynthesis, were highly expressed in IFP (Figure 7). Among 43 identified proteins exhibiting significant changes in lignin biosynthesis pathway between IFP and IVP, 17 of them were selected for qRT-PCR to check their RNA expression. All the seventeen examined genes showed higher mRNA abundance in IVP than that in IFP except gene NNU_03827-RA. Tendencies of transcription level in genes were coincided well to their translation level excluding gene NNU_03827-RA, NNU_19318-RA and NNU_03385-RA (Figure 7). These results indicated that the genes regulating lignin biosynthesis pathway in petioles are regulated similarly at the transcription and translation level, and more genes in lignin biosynthesis pathway were more regulated in IVP compared to IFP.
To evaluate the consequence of genes in lignin biosynthesis pathway being more regulated at transcription and translation level in IVP, metabolites of lignin biosynthesis pathway were quantified in the two petioles types. As shown in (Figure 7, Table S6), a total of thirteen detected metabolites were quantified in IFP and IVP. Among these metabolites, nine of them exhibited significant changes in abundance between IFP and IVP. The changes abundance of metabolites between IFP and IVP was not exactly coinciding with genes expression, but the precursor abundance of syringyl lignin monomers and hydroxyphenyl lignin monomer were significantly higher in IVP compared to IFP. Considering the higher lignin content in IVP than IFP (Figure 2), the abundance of syringyl lignin and hydroxyphenyl lignin in IVP could be higher than in IFP, and the abundance of guaiacyl lignin may not change. were quantified in IFP and IVP. Among these metabolites, nine of them exhibited significant changes in abundance between IFP and IVP. The changes abundance of metabolites between IFP and IVP was not exactly coinciding with genes expression, but the precursor abundance of syringyl lignin monomers and hydroxyphenyl lignin monomer were significantly higher in IVP compared to IFP. Considering the higher lignin content in IVP than IFP (Figure 2), the abundance of syringyl lignin and hydroxyphenyl lignin in IVP could be higher than in IFP, and the abundance of guaiacyl lignin may not change. Red color means abundance of substrate, product or enzymes were higher in IVP than IFP, green color means abundance of substrate, product or enzymes were higher in IFP than IVP and blue color means abundance of substrate and product were similar in two petioles or enzymes which catalyze same action show opposite abundance between IVP and IFP. Color bar shows changes in abundance of enzymes which were quantified in proteomics; (b,c) abundance of selected enzymes in lignin biosynthesis pathway and mRNA expression level of genes encoding selected enzymes in lignin biosynthesis pathway; (d) abundance of metabolites in lignin biosynthesis pathway. Data are means ± SD from 3 independent biologic repeats. Asterisks and double asterisks indicate significant changes compared to control as assessed according to Student's t-test (* p < 0.05 and ** p < 0.01).

Figure 7.
Overview of lignin biosynthesis pathway on transcription level, translation level and metabolite level. (a) Substrate, product and enzymes in whole pathway are written in red, green, blue and black color. Red color means abundance of substrate, product or enzymes were higher in IVP than IFP, green color means abundance of substrate, product or enzymes were higher in IFP than IVP and blue color means abundance of substrate and product were similar in two petioles or enzymes which catalyze same action show opposite abundance between IVP and IFP. Color bar shows changes in abundance of enzymes which were quantified in proteomics; (b,c) abundance of selected enzymes in lignin biosynthesis pathway and mRNA expression level of genes encoding selected enzymes in lignin biosynthesis pathway; (d) abundance of metabolites in lignin biosynthesis pathway. Data are means ± SD from 3 independent biologic repeats. Asterisks and double asterisks indicate significant changes compared to control as assessed according to Student's t-test (* p < 0.05 and ** p < 0.01).

Discussion
The molecular mechanism underlying petiole rigidity formation in lotus is important, but far from being fully clarified. Genome sequencing and analysis promote molecular biology research [33]. High throughput sequencing was commonly used in lotus to analyze rhizome formation, flowering time, SNPs and alternative splicing finding and evolution [34,35]. In the present study, gene expression at RNA and protein levels in two petioles types was studied using RNA sequencing and TMT techniques. Based on the integrative analysis of anatomic, transcriptomic and proteomic data, certain genes involved in polysaccharide and lignin biosynthesis pathway were deduced as potentially playing important role in petioles rigidity formation both in transcription and translation level. However, more work needs to be done to explore key genes involved in lotus petioles rigidity formation in future.

Genes Involved in Lotus Petioles Formation Exhibit No Synchrony at Transcriptional and Translation Level
In this study, approximately 75.36% and 18.19% annotated genes in lotus genome have found their transcripts and protein evidence in petioles (Figure 4a). A total of 93.12% proteins had their respective transcripts evidenced in transcriptome data and this ratio is similar to that in other model plant tissue such as Arabidopsis root and rice grain [26,30]. Undetected transcripts and proteins may be due to the fact that genes have different spatiotemporal expression characters or the limitation of RNA and protein extraction as well as the detection techniques. Moreover, a total of 335 proteins from proteomics data had no relevant transcripts in transcriptomics data. It was speculated that these genes may have a brief transcription time while their proteins have relatively longer life span in petioles. Similar to a previous study, our results also pointed out the un-synchronous in transcriptional and translation of genes in lotus [32]. The poor correlation of total identified transcripts and proteins suggest that a proportion of transcripts may not be translated into proteins and the higher correlation of DEGs and their corresponding DEPs suggest these genes have similar function in petioles development (Figure 5a,d). Inconsistency of genes in mRNA and protein level was similar with other previous studies [30][31][32]. Furthermore, it suggests that the transcriptomics and proteomics analysis are strategically complementary and have equal importance in lotus petioles development analysis [28,29].

Anatomic Evaluation of Lotus Petioles Indicated Cell Wall Thickness and Cell Differentiation Affect Petiole Rigidity
In spite of some anatomic evaluations being presented on Sacred lotus, no research has specifically explored the difference of the two types of petioles. In this present study, we found that IFP and IVP have similar anatomic structure. However, there are some differences between IFP and IVP. Usually, IVP had more collenchyma tissue, more vascular bundles and thicker xylem vessel cell wall (Figure 2), which enable the plant has more strength to keep its architecture [36]. Genes involved in cell wall polysaccharides biosynthesis may play key roles in its anatomic structure formation. In this study, 5 cellulose synthases, which plays important roles in cell wall formation [37,38], showed higher abundance in IVP than that in IFP. Specifically, there were more layers of collenchyma cell in IVP than that in IFP to support petioles stand upright. The vascular bundles and the thicker xylem vessel cell wall are the usual structure in stem or petioles of Arabidopsis [39], rice [40], Populus [41,42], as well as bamboo [43] to provide mechanical support. The more abundance of these structures in IVP indicate that more mechanical strength is generated. Additionally, the different amount of cellulose and lignin and the diverse structure of the cell wall of IVP and IFP may suggest that the cell wall structure could be the main reason for different petioles observed in lotus.

Abundance of Lignin Biosynthesis Related Proteins Significantly Higher in IVP Compared to IFP Supports Contribution of Lignin Biosynthesis in Petiole Rigidity
It has been established that a constant energy supply is a prerequisite for plant growth. Therefore, high carbohydrate metabolism in cell is necessary for growth of petioles. Despite the rapid growth, the floating leaf petiole was phenotypically weak, indicating less deposition of cell wall strengthening compounds and biologic process enrichment gives more clues for cell wall biosynthesis or organization related proteins, demonstrating that they are highly expressed in IVP.
Most importantly, after KEGG pathway enrichment analysis, numerous proteins related to cell wall polysaccharides and lignin biosynthesis were identified from 904 significantly differentially expressed proteins ( Figure 6, Figure S3). Phenylpropanoid and other flavonoids were reported to respond to both abiotic and biotic stimuli [44]. Additionally, flavonoids biosynthesis and lignin biosynthesis were reported have some connections in plant [45,46]. Lignin gets deposited in the matrix of cellulose micro fibril, which strengthens the cell wall. This study showed the significantly enriched phenylpropanoid biosynthesis/lignin biosynthesis pathway and identified major enzymes involved in lignin biosynthesis ( Figure 7, Table 2). Furthermore, 10 abundant enzymes that trigger lignin biosynthesis at the beginning until the formation of lignin monomer were higher in IVP than IFP, which indicates that IVP had more effective lignin deposition compared to IFP.
The transcription level of the selected enzymes genes in lignin biosynthesis pathway were checked in petioles. The results suggested that the genes involved in the regulation of lignin biosynthesis pathway appear at both transcription and translation levels exhibiting similar tendencies in both. Consequently, analyses of metabolites in lignin biosynthesis imply that more lignin was synthesized in IVP compared to IFP (Figures 2 and 7).

Abundance of Cell Wall Polysaccharides Biosynthesis Related Proteins Significantly Higher in IVP Than That in IFP
Cell wall polysaccharides biosynthesis is more complicated than lignin biosynthesis. Biologic processes relevant to cellulose and hemicellulose biosynthesis affect cell wall polysaccharides biosynthesis [47,48]. Pathways such as starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, and biosynthesis of amino acids were significantly enriched in this study (Figure 7, Table 2). Cellulose is the most abundant biopolymer deposited on the plant cell wall. Approximately 40-50% cellulose stored in secondary cell wall of wood tissue, which mainly contributes to the bodily structure [47]. Studies in Arabidopsis, Populus and other plants suggest that various genes involved in cellulose biosynthesis play fundamental roles in cell wall formation and thickening. A total of 10 cellulose synthases [49] mediate cellulose synthesis in both primary [50] and secondary walls in Arabidopsis [51,52]. Previous researches showed that some cellulose synthases in Populus participate in cell wall formation [25,53,54] and other genes including KORRIGAN1 (KOR1), sucrose synthase gene and transcription factors are also involved [42,[55][56][57]. Specifically, the abundance of four CESAs (cellulose synthases) including NNU_09561-RA, NNU_12044-RA, NNU_21632-RA and NNU_01080-RA were significantly expressed higher in IVP compared to IFP, suggesting that more cellulose is deposited in IVP (Figures 2 and 3). Similarly, abundance of three UGD (UDP-glucose 6-dehydrogenase) and one SUSY (sucrose synthase), for instance, NNU_03659-RA, NNU_07386-RA, NNU_04520-RA and NNU_19077-RA were significantly enriched higher in IVP. These proteins have similar functions in cellulose biosynthesis as their homologous proteins in Arabidopsis and Populus [58,59]. Furthermore, hemicellulose biosynthesis related proteins such as HEX (hexokinase), UXS (UDP-glucuronic acid decarboxylase), IRXs (irregular xylem), GATL (galacturonosyltransferase) and UXT (UDP-xylose transporter) have been reported to be involved in cell wall hemicellulose biosynthesis [39,[60][61][62]. In our study, several UXS (NNU_15122-RA, NNU_12302-RA, NNU_10172-RA and NNU_25253-RA), IRX (NNU_20069-RA, NNU_12476-RA, NNU_13626-RA, NNU_13619-RA, NNU_22319-RA, NNU_03759-RA, NNU_18746-RA and NNU_11544-RA) and UXT (NNU_12576-RA, NNU_24078-RA) were identified and most of them showed significantly higher abundance in IVP than in IFP, suggesting the huge contribution of hemicellulose biosynthesis in cell wall formation. The function of these proteins involved in polysaccharides biosynthesis need more experiments in vivo to confirm because the presence of glycoside hydrolases does not mean that the polysaccharides are degraded [63]. Their expression of these genes at mRNA and protein level may have more complicated relation with polysaccharides deposition in cell wall.

Plant Growth and Petiole Collection
The sacred lotus cultivar "Ancient Chinese lotus" was cultivated in Wuhan Botanical Garden, Chinese Academy of Sciences (N30 • 32 44.02", E114 • 24 52.18"). There are two ecotypes of N. nucifera: temperate lotus and tropical lotus (N19 •~N 43 • in China) [64]. "Ancient Chinese lotus" is temperate lotus. Its seed starts to germinate in April, flower blooms from June to August, leaf withers in September and October and dormant buds remains from November to March of the following year with an enlarged rhizome [65]. In middle May, the dormant buds sprout underground. In terms of leaf development, it could be divided into two stages. At floating leaf development stage, the newly developed leaf is the floating leaf, which is folded until emerging above the water surface. As the first floating leaf grows into maturity, the second leaf appears. After about several floating leaves have emergence, the succeeding leaves are vertical leaves. At the vertical leaf development stage, the newly developed vertical leaf is folded, which will grow bigger and unfold when the petiole is long enough to support leaf far above the water surface. Four kinds of petioles samples were collected at different stages. At the floating leaf development stage, initial floating leaves' petiole (IFP) was collected from floating leaf when it just floats on water surface. Mature floating leaves' petiole (MFP) was collected from floating leaf when its shape and size had stopped changing. At the vertical leaf development stage, initial vertical leaves' petiole (IVP) was collected when the folded leaf just emerges out of water surface. Mature vertical leaves' petiole (MVP) was collected when its shape and size stopped changing ( Figure S1). In this study, rhizome of "Ancient Chinese lotus" was planted in concrete container (length, 2 m; width,1 m; depth, 80 cm) with loam soil is about 50 cm in depth on 5 April 2016. The containers were constructed outside under natural light conditions. Fifty grams of fertilizer (SAN AN; STANLEY, Shandong, China) was applied to container every 2 weeks during the growing season. In field, lotus usually was cultured in farmland or pond with different depth of water (20 cm to 100 cm). The two different types of petioles were always observed in different depth of water in filed. In order to exclude the effect of water depth on ability of leaf stand erect on water, lotus was cultivated in low water depth. this means the length of two types of petioles was longer than water depth.

Petiole Microscopic Observation
Petiole segments (the middle section of the petiole) were cut and fixed in FAA (formaldehyde-acetic acid-alcohol solution) solution (5% glacial acetic acid, 5% formaldehyde and 70% ethanol) at room temperature for 24 h. The segments were then rinsed with tap water and dehydrated in ethanol baths ranges of 70% ethanol to 100% ethanol. After dehydration, ethanol was replaced by chloroform progressively. Then finely ground paraffin was added into chloroform and incubated it at 36 • C for 2 h. The small segments were transferred into xylene solution containing 50% and 75% paraffin and incubated at 42 • C and 50 • C for 2 h successively. After that, the small segments were transfer into pure paraffin solution at 58 • C for 1 h and this process was repeated three times. Finally, the sample were embedded using pure paraffin. Samples were cut with a rotary microtome Leica RM2265 (Leica, Benshein, Germany). Sections (10 µm) were stained using fast green and counterstained using safranin solution (Biosharp, Hefei, China) and observed under an optical microscope (Olympus BX53, Tokyo, Japan). Safranin appears red in lignified, suberized or cutinized cell walls. Fast Green (Biosharp, Hefei, China) presents green in cytoplasm and cellulosic cell walls. Petiole segments from three independent plants were collected, and at least 3 sections per plant were employed for microscopy observation.

Breaking Resistance Measurement
Petioles (petiole length was about 30 cm) with same diameter (5 mm) from IFP and IVP were cut into equally long stem segment (20 cm). The breaking resistance of the middle point of petiole was measured using a prostrate tester (DIK 7400, daiki rika kogyo co. ltd., Tokyo, Japan). The distance between tester fulcra was set at 10 cm (Figure 1g,h). The petiole was pressed using prostrate tester until petiole breaking point and the test readings were recorded. Breaking resistance was represented by the manual force added on petiole. Breaking resistance (F) = (test reading × 39.2 N) ÷ 40 [66,67]. Six biologic petioles of each type of petiole were tested in this experiment.

Transmission Electron Microscope
Petiole segment (the middle section of the petiole) was infiltrated in 4% (v/v) glutaraldehyde in phosphate buffered solution at 4 • C overnight (PBS, 33-mM Na 2 HPO 4 , 1.8-mM NaH 2 PO 4 and 140-mM NaCl [pH 7.2]). Then the petiole segment was fixed in 1% (w/v) osmium tetroxide (Biosharp, Hefei, China) and dehydrated through a gradient of ethanol and eventually embedded in Spurr's resin (Biosharp, Hefei, China). Petiole sections (100 nm) were stained using uranyl acetate and lead citrate. Lastly, the petiole sections from three independent petioles from three biologic lotus were visualized using HT7700 transmission electron microscope (HITACHI, Tokyo, Japan) and the cell-wall thickness was measured from ten cells on each section using Hitachi TEM system microscopy software (Gatan, Pleasanton, CA, USA).

Determination of Crude Cellulose and Lignin
Crude cellulose was measured by acid digestion method using cellulose detective kit (Sinobestbio, Shanghai, China). The weighed samples (the middle section of the petiole) were grinded with 80% ethanol and washed with ethanol and acetone. After the collected residue was digested by concentrated sulfuric acid and the supernatant was used to detect cellulose content. Lignin content was measured according to Klason lignin methods [68]. Six biologic petioles of each type of petiole were tested in this experiment.

Total RNA Extraction
Total RNA was extracted from about 0.5 g of samples with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the standard protocol. RNA quality and concentration were evaluated using the 6000 Pico LabChip of the Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) and a NanoDrop spectrophotometer (Thermo Scientific, Waltham, MA, USA), respectively.

cDNA Synthesis and Illumina Sequencing
The cDNA was synthesized by cDNA synthesis kit (Bio-Rad, Hercules, CA, USA) according to the manufacturer's instructions. After purification, fragments were enriched by PCR amplification for cDNA library construction. Then, library products were sequenced with Illumina HiSeq TM 2000 platform (San Diego, CA, USA). Each sample had three biologic replicates from three individual plants.
The raw data of this project was deposited in NCBI (www.ncbi.nlm.nih.gov/bioproject) and can be downloaded with identifier PRJNA642670.

Bioinformatics Analysis
RNA sequencing data processing, gene expression calculation, differentially expressed genes (DEGs) identification, gene ontology (GO) enrichment analysis and pathway enrichment were performed in BMK Cloud (Beijing, China, http://www.biocloud.net/). A false discovery rate (FDR) <1% and an absolute value of log2Ratio >1 was set as the threshold to judge the significance of gene expression difference between two petioles types. Generally, the raw RNA sequencing data were cleaned by removing adapter sequences, reads containing ploy-N and filtered with low-quality sequences (Q < 20) from raw data. All the downstream analyses were based on clean data with high quality. Clean reads were aligned to the reference genome sequence using the algorithm Tophat2 [69].
Tolerance parameters in Tophat2 were set as default with no more than two bases mismatches. Then gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins);Swiss-Prot (A manually annotated and reviewed protein sequence database); KO (KEGG Ortholog database); GO (Gene Ontology). Gene expression levels were estimated by FPKM (fragments per kilobase of transcript per million fragments mapped). Differential expression analysis of two groups was performed using the DESeq R package (1.10.1). DESeq provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting p-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate.

Protein Extraction and Trypsin Digestion
Petioles were ground and transferred to a 5-mL centrifuge tube containing homogenizing buffer (8-M urea, 2-mM ethylenediaminetetraacetic acid (EDTA), 10-mM dl-dithiothreitol (DTT) and 1% protease inhibitor cocktail). The protein was precipitated using cold TCA-acetone solution for 2 h at −20 • C. Then, the remaining precipitate was washed with cold acetone. The protein was re-dissolved in lysis buffer (8-M urea, 100-mM triethylammonium bicarbonate (TEAB), pH 8.0) and the protein concentration was measured by BCA methods [70]. The protein solution was reduced by 10-mM DTT for 30 min at 37 • C and alkylated with 20-mM iodoacetamide (IAM) at 37 • C in darkness for 45 min. Trypsin was added with a 1:50 trypsin-to-protein mass ratio for overnight digestion. All regents used in this section were supplied by Sangon Biotech Company (Shanghai, China). Six biologic petioles of each type of petiole were used in this experiment.

Quantitative Proteomic Analysis by LC-MS/MS
Peptides in each sample were separated into 18 fractions and dried by vacuum centrifugation (Labconco, Kansas, MO, USA). Then peptides were dissolved in 0.1% formic acid (FA) and loaded directly onto the reversed-phase pre-column (Acclaim PepMap 100, Thermo Scientific, Waltham, MA, USA). The peptides were separated by the reverse phase analysis column (Acclaim PepMap RSLC, Thermo Scientific, Waltham, MA, USA) according to the following conditions. The 6% to 22% solvent B (0.1% FA in 98% ACN, acetonitrile) for 22 min, 22% to 36% for 10 min, escalating to 85% in 5 min and held at 85% for the last 3 min. Flow rate was held at 300 nL/min on an EASY-nLC 1000 ultra-performance liquid chromatography system (UPLC)(Thermo Scientific, Waltham, MA, USA). Peptides were detected in the Orbitrap (Thermo Scientific, Waltham, MA, USA) at a resolution of 70,000 and instrument parameters were set as follows: Normalized collision energy (NCE) value was 30; ion fragments resolution was 17,500. The electrospray voltage was 2.0 kV. The m/z scan range is 350 to 1800. Fixed first mass is 100 m/z [70]. A data-dependent analysis method that alternated between primary MS scan followed by 20 MS/MS scans was applied for the top 20 precursor ions above a threshold ion count of 1E4 in the MS survey scan with 30.0 s dynamic exclusion.

Database Searching
The acquired MS/MS data were processed using Proteome Discoverer (version 2.1.0.81, Thermo Scientific, Waltham, MA, USA). Tandem mass spectra were searched against the sacred lotus genome database with mascot search engine (https://bioinformatics.psb.ugent.be/plaza/versions/ plaza_v4_dicots/). The parameters were set as follows: Trypsin/P was chosen as enzyme and the maximum missing cleavages was 2. The number of modification and charge in each peptide were set up to 5. Mass error tolerance was set to 10 ppm and 0.02 Da for precursor and fragment ions, respectively. Fixed modification included carbamidomethylation on Cys and TMT-6-plex on Lys and N-term. Variable modification contained oxidation on Met and FDR (false discovery rate) threshold value 1% was applied in protein, peptide and modification site identification. Minimum peptide length was set at 7.

Protein Annotation, Functional Classification and Enrichment
Proteins were classified by GO annotation method. GO annotation proteome was derived from the UniProt-GOA database (http://www.ebi.ac.uk/GOA/) [71,72]. KEGG database was used in the identification of enriched pathways [71,72]. A two-tailed Fisher's exact test was employed to test the enrichment of DEPs against all identified proteins. The GO or pathway with a corrected p-value < 0.05 was considered significant.

Metabolite Extraction
Freeze-dried petiole was crushed using a mixer mill (MM 400, Retsch, Germany). One hundred milligrams of powder were weighed and stirred in 70% methanol (1.0 mL). After centrifugation, the supernatant containing extracts were absorbed (CNWBOND Carbon-GCB SPE Cartridge, 250 mg, 3 mL, Shanghai, China and filtrated before LC-MS analysis. Six biologic petioles of each type of petiole were tested in this experiment.

Metabolite Analysis
Metabolite profiling was carried out using a widely targeted metabolome method by Wuhan Metware Biotechnology Co., Ltd. (Wuhan, China, http://www.metware.cn/). The sample extracts were analyzed with LC-ESI-MS/MS system. The analytical conditions were as follows. Column (Waters, 1.8 µm, 2.1 mm × 100 mm) and buffer A (water containing 0.04% acetic acid) and buffer B (acetonitrile containing 0.04% acetic acid) were used. The gradient program in HPLC system (SCIEX, Redwood City, CA, USA) was set as follows: At the first 12 min, mobile phase containing buffer A 5% and buffer B 95%. In the next three min, mobile phase was changed into buffer A 95% and buffer B 5%. The flow rate was set as 0.40-mL/min and temperature was set at 40 • C. A total of 5-µL sample was injected in each run. The eluate was alternatively connected to MS (API 4500 Q TRAP LC/MS/MS System, SCIEX, Redwood City, CA, USA).

qRT-PCR
Total RNA was extracted by TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and 2 µg total RNA was used for reverse transcription using Rever Tra Ace-α-First Strand cDNA synthesis Kit (TOYOBO, Osaka, Japan). qRT-PCR was performed on CFX96 Real Time System (BIO-RAD, Hercules, CA, USA) with SYBR green fluorescence. The data were normalized based on relative transcript level of actin (NNU_24864). Three biologic petioles of each type of petiole were tested in this experiment.

Parallel Reaction Monitoring (PRM)
The digested peptides were dissolved in 0.1% formic acid and directly loaded onto a reversed-phase analytical column (15 cm length, 75 µm i.d.). The parameters of UPLC (Thermo Scientific, Waltham, MA, USA) were set as mentioned above. The parameters of MS were set as follows: The electrospray voltage was 2.0 kV. The m/z scan range was 350 to 1000 for full scan and the resolution was 35,000. Peptides were then selected for MS/MS using NCE setting at 27 and the resolution of fragments was 17,500. The maximum IT was set at 20 milliseconds for full MS and auto for MS/MS. The isolation window for MS/MS was set at 2.0 m/z. The acquired MS data were identified and quantified using Skyline software (v.3.6) [73]. MS proteomics data in this study were deposited in the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the PRIDE partner repository [74]. Data are available via ProteomeXchange with identifier PXD009136.

Statistical Analysis
The data were analyzed using the Student's t-test to evaluate the statistical significance to compare the two groups, and one-way ANOVA test to compare multiple groups (p < 0.05).

Conclusions
In lotus growth, the floating leaves and erect rigid leaves develop subsequently. These two petioles types appear like lateral branches growing during its growing season. Due to the fact that floating leaf and vertical leaf experience similar ambient, the differentiation of petioles rigidity in lotus is not affected by environmental changes. To investigate the mechanism controlling the rigidity of lotus petiole, we performed integrated omics analysis in this study. Anatomically, more vascular bundles, collenchyma and thicker cell wall were observed in IVP than that in IFP. Many genes involved in cell wall biosynthesis and lignin biosynthesis exhibited differentially expression at both mRNA and protein levels. Gene function and pathway enrichment analysis displayed that DEGs and DEPs were significantly enriched in cell wall biosynthesis and lignin biosynthesis. Interestingly, the bulk of identified DEPs in lignin biosynthesis were up regulated in IVP, suggesting that the differences in lignin biosynthesis in the two types of lotus petioles may be responsible for the observed different rigidity. Despite our results are still far from underlying the rigidity of lotus construction, these findings afford novel clues for better understanding of the gene network involved in lotus petioles formation. In addition, roles of decisive candidate gene involved in cell wall and lignin biosynthesis in lotus should be accurately interpreted in future studies. The 905 DEPs were classified into four groups. Abundance of 164 proteins in IFP was two-fold higher than IVP (Q1) and abundance of 257 proteins in IFP was 1.5-2-fold higher than IVP (Q2). Also, abundance of 346 proteins in IVP was 1.5-2-fold higher than IFP (Q3), and abundance of 137 proteins in IVP was two-fold higher than IFP (Q4); (c) KEGG pathway enrichment analysis. A total of 105 proteins were enriched in six different pathways such as starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, glycosphingolipid biosynthesis, phenylpropanoid biosynthesis, flavonoid biosynthesis and biosynthesis of amino acids, Figure S4: Validation of proteomics data using PRM. (a) Comparison of proteomics data and PRM validation; (b) Correlation of proteomics data and PRM validation; (c-f) show that the abundance of two peptides (NNU_21632-RA andNNU_03827-RA) belongs to cellulose synthase and phenylalanine ammonia-lyase, Figure  S5: Function analysis of 135 core genes. (a) Gene function analysis of 135 core genes by MapMan bincodes; (b) Cluster analysis of 135 core genes. Asterisks indicate that DEGs and DEPs have the opposite expression of RNA and protein level. Dot indicates that core genes involved in cell wall biosynthesis. Red dot means DEGs and DEPs have the opposite expression of RNA and protein level, and black dots mean DEGs and DEPs have the same expression of RNA and protein level; (c,d) indicate the function and pathway enrichment analysis of 135 core genes, Figure S6: GO enrichment analysis of DEGs and DEPs. (a,b) GO enrichment analysis of DEPs and DEGs; (c) DEPs involved in pathways related to cellulose and hemicellulose biosynthesis. Red color characters mean abundance of proteins were higher in IVP than IFP, green color characters mean abundance of proteins were higher in IFP than IVP, Table S1: A total of 1422 DEGs in IFP and IVP were identified,  Acknowledgments: The authors sincerely thank Rebecca Njeri Damaris at Hubei University and Tonny Maraga Nong'a at University of New Mexico for English writing optimization.

Conflicts of Interest:
The authors declare no conflict of interest.