WGCNA Identifies a Comprehensive Gene Co-Expression Network Associated with Development in autotetraploid potato tissues

1. Introduction

Potato (Solanum tubulosum L.) is the most consumed non-cereal food and vegetable crop and fourth largest food crop in the world after maize, rice, and wheat [1]. It is also an excellent source of energy, protein, vitamins, and minerals. Potato species with 24, 36, 48, 60, 72 polyploidy exist in nature [2]. Among them, most of the commercially cultivated potato species are autotetraploid. Autotetraploid potato has a high frequency of chromosome recombination and a low frequency of recessive gene expression, which can be fixed by asexual propagation. However, its genome is highly heterozygous, which hinders genetic analysis and improvement. Therefore, it is important to investigate and identify molecular markers and key genes for the genetic development of autotetraploid potato. High-throughput bioinformatics approaches can be used to map genes of interest more quickly and easily than traditional molecular biology techniques [3,4]. In 2011, the genome of the heterozygous diploid potato RH (S. tuberosum group) was sequenced and assembled by the International Potato Sequencing Team (PGSC) and it provided a reference database for discovering useful genetic information from massive biological data [1].

Weighted gene co-expression network analysis (WGCNA) is a method to specify modules based on the expression data of gene chip or RNA sequencing (RNA-Seq). It uses system biology to understand co-expression networks and explore connections between genes and target traits [4]. Moreover, WGCNA provides clues for screening the signaling pathways and regulatory factors involved in the target genes [5]. At present, a large number of literatures have used WGCNA to study potato disease resistance and defense mechanisms. Yan et al. used different hormones to treat late blight resistant potato genotype SD20 and used WGCNA to analyze transcriptome data to obtain nine defense genes, which provided reliable targets for functional verification [6]. Jing et al. used WGCNA analysis to find 13 TF-hub genes related to proline metabolism and salt stress resistance in potato [7]. Guo et al. analyzed gene expression in potato tissues under nitrogen-sufficient and nitrogen-deficiency conditions, and generated 23 modules from 116 differentially expressed hub genes involved in photosynthesis, nitrogen metabolism, and secondary metabolites by WGCNA [8]. Cao et al. identified number of potato hub genes by WGCNA, including genes that regulate plant immune responses such as flagellin-sensitive 2 and chitinases [9]. Additionally, studies have also used WGCNA to explore the molecular mechanisms of plant growth and development in other species. Liu et al., constructed a co-expression network of genes in rutabaga hypocotyl tuber, and WGCNA identified 59 co-expressed gene modules. Among them, two genes, Bra-FLOR1 and Bra-CYP735A2, related to tuber growth were grouped in the same module. This module also included genes involved in carbohydrate transport and metabolism, cell wall growth, auxin regulation, and secondary metabolism [10]. To date, only few studies have investigated the growth and development of potato by WGCNA.

In this study, we aimed to explore the co-expression network of genes involved in potato growth and development. Firstly, transcriptome sequencing was performed on roots, stems and leaves of the autotetraploid potato JC14 at three developmental stages, and differentially expressed genes (DEGs) were analyzed. Secondly, the Gene Co-expression Network was constructed by WGCNA, and the target gene modules related to growth and development traits were screened. Thirdly, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene ontology (GO) analysis were performed to reveal the major metabolic pathways involved in the target gene modules and the potential functions of these genes. The key genes in the module were then identified according to the connectivity of the genes in the corresponding network. Finally, real-time quantitative PCR (RT-qPCR) was used to verify the expression pattern of these key genes during growth and development in Solanum tubulosum L. This study was expected to construct a comprehensive and dynamic gene co-expression network of potato growth and development, reveal the spatio-temporal expression of genes and the regulation mechanism of traits at different developmental stages, lay a theoretical foundation for exploring the molecular mechanism of potato growth and development regulation, and provide new genetic resources for molecular breeding of autotetraploid potato.

2. Materials and Methods

2.1. Plant material culture and sampling

Potato genotype JC14 was cultivated in Nankou Pilot Base, Chinese Academy of Agricultural Sciences, Beijing (longitude 116.408490; Latitude 40.154860) in a greenhouse with natural light from 2022 to 2023. Cultivation and field management followed normal production protocols. Potato roots, stems, and leaf tissues and organs were sampled at the seedling, tuber formation, and tuber expansion stages in three biological replicates. After collection, the samples were washed with pure water, quickly froze in liquid nitrogen, and kept at -80°C for subsequent total RNA extraction.

2.2. RNA extraction, Illumina sequencing and data analysis

Total RNA from 27 samples was extracted using TIANGEN's plant RNA extraction kit by magnetic bead method (item number: DP762-T1) (TIANGEN, Beijing, China) according to the manufacturer's protocol.Total RNA purity was detected using a Nano Photometer® spectrophotometer (IMPLEN, CA, USA) and total RNA samples were tested for concentration and integrity using Agilent 2100 RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA) in the Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). When RNA integrity RIN >= 7, the requirements of cDNA library construction and subsequent sequencing was fulfilled. After the total RNA samples were qualified, 1-3μg of which from each sample was used as starting material to construct the transcriptome sequencing library. According to the VAHTS Universal V6 RNA-seq Library Prep Kit for Illumina (No. NR604-01/02) (Vazyme, Nanjing, China), different index tags were selected to build the sequencing library. For qualified Total RNA samples, the eukaryotic mRNA was enriched by magnetic beads with Oligo (dT), and then the mRNA was broken into short fragments by adding fragmentation buffer. One-strand cDNA was synthesized with random hexamers, followed by adding buffer, dNTPs, RNase H and DNA polymerase I to synthesize double-stranded cDNA. The double-stranded cDNA was then purified using AMPure P beads purification or QiaQuick PCR kit. The purified double-stranded cDNA was subjected to end repair, A-tail addition and sequencing adapter ligation, followed by fragment size selection, and finally PCR enrichment to obtain the final cDNA library. After the library was constructed, the Qubit 3.0 was used for preliminary quantification, and the library was diluted to 1 ng/μL, and finally Bio-RAD KIT iQ SYBR GRN was used for QPCR to accurately quantify the effective concentration of the library (library effective concentration > 10 nM) to ensure library quality. The NovaSeq 6000 S4 Reagent kit V1.5 was used for clustering and sequencing on the NovaSeq 6000 S4 platform, and the double-end sequencing program (PE) was run to obtain 150 bp double-end sequencing reads.

Sequencing quality was assessed using FastQC and Trimmomatic software [11]. Raw reads were trimmed by eliminating adapters, poly-N, and low-quality reads to obtain clean data, which were then mapped to the potato reference genome using HISAT2 [12]. The reference genome was downloaded from Spud DB Potato Genomics Resource (http://spuddb.uga.edu/rh_potato_download.shtml). Transcription and expression level of genes were measured using FPKM (fragments per KB of transcript per million fragments mapped) in StringTie [13] and differential expression was analyzed using the DESeq R package [14]. In this study, a DEG was considered if the fold change was >2 and q-valuewas <0.05.

2.3. Construction of weighted co-expressed gene networks

By referring to the tutorial on the WGCNA official website, The WGCNA package of R software was used to construct the co-expressed gene network [4]. Firstly, samples with low correlation or that could not be clustered on the dendrogram were removed and cluster analysis was performed. Secondly, in order to meet the premise of scale-free network distribution, the pick soft threshold function in the WGCNA package was used to calculate the soft threshold. When the fitting curve was close to 0.9 for the first time, the threshold parameter β was selected. Then, according to the value of β, the correlation based association between phenotype and gene modules was performed to create the adjacency matrix. The adjacency matrix was further transformed into topological overlap matrix (TOM) to construct the gene connectivity network. Finally, based on the eigengenes of each module, the gene modules are identified and clustered by the dynamic tree cut method, and the modules that are closer to each other are merged into new modules. In this study, the module similarity threshold was set to 0.25, the expression threshold was set to 2, and the minimum number of genes in a module was set to 30 to partition modules.

2.4. Selection and enrichment analysis of key modules

To further investigate the gene modules related to potato growth and development, the correlation coefficients between the module characteristic genes and different samples were calculated. The module with a correlation coefficient greater than 0.5 was selected as the target gene module. DEGs in the target module were annotated using Blast2Go [15]. Blast2Go was used to retrieve gene ontology (GO) entries, including biological processes, molecular functions and cellular components. Pathway enrichment analysis of DEGs within the target module was performed using KEGG database [16].

2.5. Screening and functional analysis of hub genes

All DEGs of the four target modules were selected and using the WGCNA R package to calculate the intra-module gene significance, edges, and nodes [4] (Table S5). The information was exported to Cytoscape 3.9.1 software for gene network analysis and mapping [17]. DEGs with high connectivity in the target gene module were screened as candidate hub genes. Swiss-port databases [18] was used for homology annotation.

2.6. RT-qPCR was used to verify hub gene expression

In order to verify the reliability of RNA-seq data, 20 genes were used for RT-qPCR to determine the changes in their relative expression levels, and R package was used to calculate and plot the correlation between the relative expression levels and FPKM. And then, a total of 40 candidate hub genes were selected for RT-qPCR analysis to assess their expression levels. Primers for candidate hub genes were designed using Primer Premier 5.0 (Table S8). Potato stems from the seedling, tuber formation and tuber expansion stages were selected as raw materials. Gene expression analysis of candidate hub genes was performed on the LightCycler480 system (Roche, Rotkreuz, Switzerland) with a total reaction volume of 10 µL for RT-qPCR, this included 12.5 µL 5 µL 2X SG Fast qPCR Master Mix (Roche, Rotkreuz, Switzerland), 0.2 µL primer (10 µM), 1.0 µL template (10×diluted cDNA liquid), and 3.6 µL ddH2O. Using actin as the reference gene and the 2−∆∆ CT methods, the expression levels of selected hub genes were calculated [19]. Histograms were plotted by seaborn package in python. Three independent biological replicates were performed for each sample in this assay, significance and standard deviation was determined by independent t-test of statannot package (ns: p > 0.05; *: p < 0.05, **: p < 0.01; ***: p < 0.001; ****: p < 0.0001).

3. Results

3.1. Genome-wide transcriptome analysis

High-throughput RNA-Seq was performed on roots, stems, and leaves at three different developmental stages with three biological replicates per sample, resulting in 176.36 GB of clean reads. For each sample, the sequencing data was larger than 5.37 GB, the GC content ranged from 42.77% to 45.87%, and the percentage of Q30 bases was greater than 91.89% (Table S1). Compared with the potato (RH89-039-16) reference genome [20], the reads alignment of samples was between 66.47% and 79.37% (Table S1) and this indicated that the quality of the sequencing data was high and met the requirements for subsequent analysis. In potato root, 209, 953 and 1,018 up-regulated differential genes and 13, 380, 12, 765 and 924 down-regulated differential genes were identified between seedling stage and tuber formation stage, seedling stage and tuber expansion stage, tuber formation stage and tuber expansion stage, respectively (Figure 1A). In the stem, 113, 534 and 930 up-regulated genes and 144, 800 and 817 down-regulated genes were also discovered between any two of the three stages, respectively (Figure 1B). In leaves, 6, 701 and 982 up-regulated genes and 34, 783 and 765 down-regulated genes were screened between the two at the three stages, respectively (Figure 1C).

3.2. Construction of weighted gene co-expression network and identification of modules

After filtering genes with variance ≥ 0.5, 14, 475 genes were used to construct the weighted gene co-expression network. All 27 samples were clustered and the correlation coefficient of expression levels was calculated for each sample. Overall, the samples clustered well and no outliers appeared (Figure 2). The weight values were calculated using the pick soft threshold function included in the WGCNA software package, and the soft threshold β=28 was determined when the fitted curve was first close to 0.9 (Figure 3A). Modules with similar expression are then merged by the dynamic tree cut method. Genes with high correlation are assigned to the same module, and different branches of the cluster tree in the figure represent different gene co-expression modules, and different colors represent different modules, which are divided into 12 modules in total (Figure 3B). A total of 1, 546 genes were obtained from the 12 co-expressed gene modules. Among them, there were eight modules with 150 genes, including black, cyan, darkgreen, darkorange, grey60, midnightblue, red and steelblue and the plum1 module contained the lowest number of genes of 41 (Table S2).

The correlation coefficients between characteristic gene modules and traits were tested to explore the relationship between the identified modules and potato growth and development. Based on the correlation between the trend of gene expression modules and tissue site enrichment, the Pearson correlation coefficient (r2 >0.5) and the magnitude of significance p-value (P<0.05) were used as the screening criteria. Genes in midnightblue module was significantly enriched in roots at the tuber formation stage and plum1 was significantly enriched in roots at the tuber expansion stage. In comparison, darkturquoise, grey60, and darkgreen were significantly enriched in stems at the seedling stage, and darkgreen was significantly enriched in stems at the tuber formation stage, whereas steelblue and paleturquoise were significantly enriched in stems at tuber expansion stage. In addition, lightsteelblu1 was significantly enriched in leaves at the seedling stage, while cyan was significantly enriched in leaves at the tuber expansion stage. These results suggest that the significantly enriched modules differ across organs during the same developmental period and that the regulation of gene expression is tissue specific. Furthermore, the same tissue also has different enrichment modules and differential gene expression patterns at different developmental stages (Figure 4). Subsequently, differentially expressed genes in the darkturquoise, grey60, darkgreen, steelblue, and paleturquoise modules related to potato stem growth and development were selected for GO annotation and KEGG enrichment analysis. However, since the darkturquoise module at the seedling stage did not enrich any pathways, its subsequent analysis was therefore not conducted.

3.3. GO and KEGG enrichment analysis of key modules in different tissues

Biological process (BP), cellular component (CC) and molecular function (MF) are the three GO domains. For the four target gene modules grey60, darkgreen, steelblue and paleturquoise, the three most significantly enriched GO terms in the BP entry were metabolic process, cellular process and monomer process. CC entry contains cells, cell parts and organelles, whereas MF entry has binding, catalytic activity, and transporter activity (Figure 5A, Figure 6A, Figure 7A, and Figure 8A and Table S3). The grey60 module is enriched in olefin biosynthetic process (GO: 1900674), defense response (GO: 0006952), and chitin metabolism (GO: 0006030). The darkgreen module was enriched in carbohydrate metabolism (GO: 0005975), cell growth (GO: 0016049), cell wall polysaccharide biosynthesis (GO: 0070592), acyltransferase activity (GO: 0016746), methylammonium transmembrane transporter activity (GO: 0015200), etc. steelblue module was enriched in carbohydrate metabolism (GO: 0005975), growth regulation (GO: 0040008), thiamine metabolism (GO: 0006772), hydrolase activity (GO: 0016787), glucosidase activity (GO: 0015926), etc. Paleturquoise module was mainly enriched in amino acid metabolism (GO: 0006520), primary metabolic process (GO: 0044238), steroid metabolism (GO: 0008202), carbon-nitrogen ligase activity (GO: 0016884), acyltransferase activity (GO: 0016746), etc. In general, DEGs in the four target gene modules could be enriched in metabolic pathways related to growth and development, such as carbohydrate metabolism, cell growth, and carbon-nitrogen ligase activity.

KEGG enrichment analysis of DEGs in the four target gene modules (grey60, darkgreen, steelblue, and paleturquoise) identified 3, 5, 13, and 2 enriched pathways (FDR < 0.05), respectively (Figure 5B, Figure 6B, Figure 7B, and Figure 8B). The pathways enriched by grey60 module were mainly zeatin biosynthesis (map00908), MAPK signaling pathway (map04016) and plant-pathogen interaction (map04626). The darkgreen module was enriched in fatty acid elongation (map00062), phagosome (map04145), amino sugar and nucleotide sugar metabolism (map00520), biosynthesis of nucleotide sugars (map01250), and ascorbate and aldarate metabolism (map00053). The steelblue module was mainly enriched in thiamine metabolism (map00730), tyrosine metabolism (map00350), tropane, piperidine and pyridine alkaloid biosynthesis (map00960), citrate cycle (TCA cycle) (map00020), phenylalanine metabolism (map00360), glycine, serine and threonine metabolism (map00260), amino sugar and nucleotide sugar metabolism (map00520), isoquinoline alkaloid biosynthesis (map00950), starch and sucrose metabolism (map00500), β-alanine metabolism (map00410), etc. The main enriched pathways in the paleturquoise module were carbon metabolism (map01200) and biosynthesis of amino acids (map01230) (Table S4). In general, DEGs in seedling stems were enriched in defense-related metabolic pathways, such as MAPK and hormone signaling, plant-pathogen interaction, and so on. In tuber formation stage, DEGs were enriched in amino sugar and nucleotide sugar metabolism pathways, whereas in tuber expansion stage, DEGs were enriched in growth-related metabolic pathways, such as starch and sucrose metabolism, carbohydrate metabolism, and amino acid metabolism.

3.4. Screening and functional annotation of hub genes in the module

In each of the four target gene modules, 30 genes with high connectivity were selected as candidate hub genes for co-expression network mapping, in which the color from red to orange indicates the connectivity from high to low (Figure 5C, Figure 6C, Figure 7C and Figure 8C). Ranked according to the connectivity level, 10 candidate hub genes out of 30 candidates from each of the four modules (40 in total) were selected for subsequent analysis and discussion (Table S6). Swiss-port databases was used for homology annotation, most of which were related to carbohydrate metabolism and plant growth and development. In addition, the expression levels of hub genes obtained by RNA-Seq were made into a heat map (Figure 5D, Figure 6D, Figure 7D and Figure 8D and Table S7). The four modules had genes involved in carbohydrate metabolism and cell division, such as reduced wall acetylation 3 (RWA3, RHC02H1G1413.2) and polygalacturonase (PG, RHC08H1G2453.2) and zeatin O-glucosyltransferase (RHC05H1G02002.2). Hub genes related to defense response included glucanendo-1, 3-β-glucosidase (GLU, RHC10H1G2629.2) and Cysteine-rich repeat secretory protein 55 (CRRSP55, RHC09H1G1673.2), and hub genes related to transcription regulation included transcription factors MYB61 (RHC01H1G3797.2) and PCF3. In addition, DUF642, a recurrent conserved gene in the stems of many plants.

3.5. Validation of the hub genes by RT-qPCR

In order to verify the reliability of RNA-seq data, 20 genes were selected for RT-qPCR to determine changes in their relative expression levels. There is a good correlation between the results of RT-qPCR detection and RNA_seq (R = 0.70), which proves the reliability of the transcriptome analysis results (Figure 9). In addition, RT-qPCR was used to amplify the aforementioned 40 hub genes. Compared with the seedling stage, 14 hub genes were up-regulated at the tuber formation stage and 27 at the tuber expansion stage (Figure 10 and Table S7), which was consistent with the results of RNA-Seq (Figure 5D, 6D, 7D and 8D and Table S7). This indicates that these genes play a promoting role during the growth and development of potato stems. In comparison, 12 hub genes were down-regulated at the tuber formation stage and 10 at the tuber expansion stage (Figure 10 and Table S7), suggesting that these genes may be negative regulators of the growth and development of potato stems.

4. Discussion

Genetic analysis and improvement are obstructed by the high degree of genomic heterozygosity in cultivated cultivars of autotetraploid potatoes. Identifying key candidate genes related to growth and development can provide target gene resources for molecular breeding. Previous literatures have studied the mechanisms of potato cold resistance [21], disease resistance [22], and nitrogen utilization [8] through transcriptome approach. The mechanisms involved in the growth and development of autotetraploid potato, however, are extremely complex, gene expression changes in multiple metabolic pathways occur at different developmental stages. In this study, WGCNA was used to construct a co-expression network of weighted genes related to potato growth and development traits based on transcriptome datasets from roots, stems and leaves of potatoes at different developmental stages. A total of 12 gene modules were identified, of which 4 were significantly associated with stem growth and development. Thirty candidate hub genes with high connectivity were selected from each of the four modules to make a network diagram (Figure 5C, Figure 6C, Figure 7C and Figure 8C), and the top 10 hub genes in each module with the highest connectivity were selected for subsequent analysis. These hub genes may play important roles in potato growth and development through carbohydrate metabolism, cell growth, defense response and other related metabolic pathways.

Carbohydrate is an important energy substance and plays a crucial role in plant growth and development. Xylan is the most abundant hemicellulose, and the accuracy and location of xylan acetylation are critical for xylan function and plant growth and development. RWA protein family plays an important role in regulating cell wall acetylation after xylan synthesis and are essential for the efficient binding of xylan to cellulose [23]. Over-expression of Dendrobium officinale RWA3 significantly increased the acetylation level of polysaccharides in seeds, leaves and stems [24]. PG is involved in carbohydrate metabolism and is very important for cell wall metabolism and maintenance of plant tissue morphology. Literature suggests that PG is associated with quality processes such as fruit ripening, flavor, firmness [25]. Moreover, PG has a negative regulatory role in promoting fruit softening. PG is highly expressed in immature fruits and plays an important role in maintaining the firmness of fruits [26]. In addition, the up-regulation of genes related to cellulose and pectin metabolism pathways, including PG, will increase the tensile stress of Catalpa bungei, which is conducive to the formation of high-quality wood [27]. In this study, RWA3 and PG were highly expressed at all three stages, indicating their important involvement in the whole growth and development stages of the potato stem.

Auxin and cytokinins (CKs) are essential hormones that regulate plant growth and development. O-glucosyltransferase is involved in the metabolic process of plant CKs and catalyzes the reversible inactivation of O-glucosylation. It is an important mechanism for the homeostasis of plant CKs in vivo. Both O-glucosyltransferases and UDP glucosyltransferases are subfamilies in the glycosyltransferase superfamily. Among them, zeatin O-glucosyltransferase belongs to the CK O-glucosyltransferase superfamily [28]. Zeatin O-glucosyltransferase is a negative regulator of growth and development, and its up regulation may be the reason for the long life and slow growth of Dracaena [29]. In this study, the expression of O-glucosyl transferase was low in the stem at the tuber formation stage and high at the tuber expansion stage, which may contribute to the gradual slowing of potato growth as development progressed. However, UDP-glycosyltransferases (RHC01H1G4299.2) showed a low expression level at all three developmental stages.

Pathogenesis-related proteins are specific proteins induced under pathological conditions and play an important role in plant disease resistance. GLU is a pathogenesis related protein in the pathogenesis-related family. Its expression is highly induced several hours after soybean infection with Phakopsora pachyrhizi. Soybean PEPPA-7431740 may trigger immunity by interfering with the activity of -1, 3-β-glucosidase within glucan and inhibiting pathogen-associated molecular patterns [30]. The pesticide GLY-15 has a good antiviral activity against tobacco Mosaic virus in plants. It induces a stress response in plants, in which the expression of 12 key genes including glucan-1, 3-β-glucosidase is significantly up-regulated [31]. Other studies have shown that glucan endo-1, 3-β-glucosidase is involved in plant response to soil alkaline stress [32]. 1, 3-β-glucan is a polysaccharide widely distributed in the cell walls of several phylogenetically distant organisms, such as bacteria, fungi, plants and microalgae, and is involved in cell wall degradation [33]. CRRSP55 is a disease-related protein. MgMO237 interacted with rice defense related protein 1, 3-betaglucan synthase component (OsGSC) and OsCRRSP55 to inhibit plant basal immunity during the late parasitic stage of nematode. Thus, it promotes parasitism in gramineae [34]. In this study, GLU and CRRSP55 genes were expressed in stems at all three developmental stages (Figure 5A), suggesting that these genes may play an ongoing role in potato defense against pests and diseases.

In addition, MYB61 and PCF3 were selected as transcription factors in the module. Sugarcane MYBs play an important role in stem development and various stress responses, among which MYB61 positively responds to drought stress [35]. Seedlings of barren-tolerant wild soybean are more adaptive to phosphorus stress than common wild soybean, which absorb more phosphorus by increasing root length, and MYB61 is a key transcription factor for resistance to low phosphorus stress in barren-tolerant wild soybean [36]. Rice MYB61 may regulate the efficiency of crop nitrogen usage, which in turn affects rice growth and yield [37]. TCP (teosinte branched1/ccincinnata/proliferating cell factor) is a set of specific transcription factors plays an important role in plant growth and development, and PCF3 is a TCP family transcription factor [38]. TCP transcription factors regulate leaf curvature, flower symmetry and the synthesis of secondary metabolites. Most TCPs in Ginkgo biloba also respond to exogenous hormones such as ABA, SA and MeJA [39]. In this study, the expression of MYB61 was higher in stems at the seedling and tuber formation stages and decreased at the tuber expansion stage, suggesting that MYB61 may play a role in response to drought and low phosphorus stress and nitrogen fixation in the early growth stage of potato. As growth and development progressed, the expression of PCF3 in the stem gradually increased and was highest at the tuber expansion stage, suggesting that PCF3 may promote potato stem growth.

5. Conclusions

In this study, weighted gene co-expression networks of roots, stems, and leaves at different developmental stages were constructed, and 12 target modules were identified. Four target gene modules were associated with stem growth and development. GO and KEGG enrichment analysis showed that DEGs were mainly enriched in defense-related metabolic pathways, such as MAPK and hormone signal transduction, plant-pathogen interaction, and growth-related metabolic pathways including starch and sucrose metabolism, carbohydrate metabolism, and amino acid metabolism. Carbohydrate metabolism-related genes RWA3 and PG were highly expressed at all three stages, indicating their important involvement in the whole growth and development stage of stem. The low expression of zeatin O-glucosyl transferase in stems at tuber formation stage and high expression in tuber expansion stage may contribute to the slow growth rate of potato during development. GLU and CRRSP55 genes were expressed in stems at all three developmental stages, suggesting that these genes may play an ongoing role in potato defense against pests and diseases. In addition, MYB61 and PCF3 transcription factors were identified in the target module, and their regulation of downstream target genes affects the life process of potato. The results of this study will help to further understand the molecular mechanism of potato growth and development and provide new genetic resources for future breeding of autotetraploid potatoes.