Identification of the shade-responsive core DEGs in M.sinostellata
A total of 181,902 genes were detected in this transcriptome sequencing, and the average comparison rate with the compared gene set was 80.49%. 246481 non-redundant transcript sequences were obtained. The total length of these sequences is 270,112,156, the longest sequence is 20445bp, the shortest sequence is 199bp, and the average sequence length is 3420.95bp. In this study, we analyzed the global gene expression profiles of M. sinostellata for shade response using five datasets, including CK-D0, CK-D5, ST-D5, CK-D15 and ST-D15, the details of which are shown in Table S2. Each dataset contains 3 replicates. Correlation analysis was performed among these 15 samples, which showed good reproductivity in the same group and a significance difference between control groups and shade treated groups (Figure S1). This analysis suggests that these samples can be used for subsequent screening and analysis.
From a total of 181,902 expressed genes, 4734 core DEGs were identified using cross-compared venn diagrams including CK-D0 vs ST-D5, CK-D0 vs ST-D15, CK-D5 vs ST-D5 and CK-D15 vs ST-D15 (Log2FC > 1, Qvalue < 0.05) (Fig. 1). Then, we characterized these 4734 core DEGs to get more insight into their contributing molecular pathways. First, 4734 DEGs were subjected to GO analysis, which were classified into three groups and 33 subgroups (Fig. 2A). The biological process group can be divided into 18 subgroups, among which cellular process and metabolic process were the top two sub-groups involved the most genes. Four subgroups were related to cellular component, among which cellular anatomical entity and intracellular were the main subgroups involved most genes. Eleven subgroups constitute the molecular function group, and the catalytic activity and binding involved the most genes. GO enrichment analysis showed that most significant GO terms are related to photosynthesis. The top five enriched GO terms were (GO:0009522); photosynthesis, light harvesting (GO:0009765); photosystem (GO:0009521); photosynthesis (GO:0015979) and photosystem II (GO:0009523) (Fig. 2B). Then, we performed KEGG analysis to explore contributing pathways of 4734 DEGs (Fig. 3A). For KEGG classification, genes were annotated into five groups and 19 subgroups. Among metabolism group, global and overview maps, carbohydrate metabolism and energy metabolism subgroups were significantly enriched. Photosynthesis-antenna proteins was the most significantly enriched KEGG pathway. The top five frequently enriched pathways were as follows: photosynthesis-antenna proteins, galactose metabolism, phenylpropanoid biosynthesis, starch and sucrose metabolism and alanine, aspartate and glutamate metabolism (Fig. 3B).
WGCNA of core shade responsive DEGs in M. sinostellata
To identify key shade responsive genes in M. sinostellata, WGCNA was performed to analysis 4734 core DEGs, which can identify modules of highly related genes. The selection of an optimal soft threshold power is an essential step to construct WGCNA. A network topology research of 1–20 was executed, and the scale independence and the average connectivity of the WGCNA relative equilibrium were determined. Threshold 18 was selected to construct a hierarchical clustering tree of DEGs (Fig. 4A). MEDiss Thres was set to 0.22 to merge similar modules and 4 modules was generated (Fig. 4B). Genes in grey module that cannot be assigned to any modules were not analyzed in the further study. Four different modules were generated by WGCNA including lightsteelblue1, paleturquoise, darkolivegreen and violet. The darkolivegreen module, lightsteelblues1 module, paleturquoise module and violet module including 2481, 2008, 75 and 71 DEGs, respectively (Table 1). To identify co-expression similarity of modules, characteristic genes were calculated and clustered according to their correlation (Fig. 4C). We found that these 4 modules are divided into two categories: the first included lightsteelblue1 and paleturquoise modules; the second included darkolivegreen and violet modules. Gene modules of the same category may have similar functions or contributing to the same regulatory mechanism. To investigate modules associated with shade treatment, we plotted module-trait relationship heat map (Fig. 4D). This result shows that lightsteelblue1 module positively correlated with two shade-treated groups (ST-D5, r = 0.7; ST-D15, r = 0.49) but negatively correlated with control groups (CK-D0, r=-0.40; CK-D5, r=-0.38; CK-D15, r=-0.41). Similarly, paleturquoise module had a strong correlation with ST-D15 (r = 0.68), while had a negative correlation with CK-D0 (r=-0.26), CK-D5(r=-0.25) and CK-D15 (r=-0.25). In contrast, darkolivegreen showed a positive correlation with control groups (CK-D0, r = 0.38; CK-D5, r = 0.35; CK-D15, r = 0.49) but illustrated a significant negative correlation with treated groups (cor<-0.5, p < 0.05). However, genes in violet module showed no significant correlation with treated groups or control groups. These results suggest that DEGs in lightsteelblue1 and paleturquoise modules mainly up-regulated and DEGs in darkolivegreen module down-regulated under shade treatment. These results showed that lightsteelblue1, paleturquoise and darkolivegreen modules were significantly correlated with shade treatment. To obtain further understand of the expression pattern of these three modules, heat maps of gene expression for these modules were generated along with eigengene expression values (Fig. 5). We observed that lightsteelblue1 module genes were significantly responsive to shade treatment and mainly up regulated. In paleturquoise module, genes were slightly induced in ST-D5 while markedly up-regulated in ST-D15. In contrast, darkolivegreen genes mainly downregulated under shade treatment.
Table 1
Details of four detected modules
Module
|
Gene number
|
lightsteelblues1
|
2008
|
paleturquoise
|
75
|
darkolivegreen
|
2481
|
violet
|
71
|
GO and KEGG pathways analysis of DEGs in key modules in M. sinostellata associated with shade
The functions of genes under shade treatment in specific modules and its contributing regulatory pathways were revealed by GO and KEGG analysis. GO analysis suggested that 'N-acetyltransferase activity' and 'acetyltransferase activity' as the most significantly enriched GO terms in lightsteelblue1 module, which means genes in this module mainly response to shade by regulating acetyltransferase activity (Fig. 6A). KEGG analysis in lightsteelblue1 module identified 'Valine, leucine and isoleucine degradation' and 'Amino sugar and nucleotide sugar metabolism' as the most enriched regulatory pathways (Fig. 6B). Concerning paleturquoise module, GO analysis displayed that 'hydrolase activity' and 'hydrolyzing O-glycosyl compound' were the top two enriched GO terms indicating that these genes regulating hydrolase-related metabolism (Fig. 6C). KEGG analysis demonstrated that 'Starch and sucrose metabolism' pathway was appreciably influenced by shade (Fig. 6D). These results matched GO and KEGG analysis conclusion of 4734 core DEGs, indicating that the results are reliable. In darkolivegreen module, GO analysis indicated that enriched GO terms of cellular component, biological process and molecular function were mainly related to photosynthesis, among which the top five GO terms were 'photosystem I', 'photosynthesis, light harvesting', 'photosystem', 'photosynthesis' and 'thylakoid' (Fig. 6E). KEGG analysis showed that 'Photosynthesis-antenna proteins' and 'Photosynthesis' were the most enriched pathways (Fig. 6F). These indicate that genes in darkolivegreen module mainly involved in photosynthesis.
Hub genes identification in key Co-Expressed modules in M. sinostellata under shade
To identify hub genes in these three modules, genes with a weight parameter over 0.4 were analyzed and visualized through Cytoscape 3.8.2 to construct interaction networks (Fig. 7). We identified 6, 7 and 8 hub genes in lightsteelblue1, paleturquoise and darkolivegreen modules respectively via integrated analysis results of MCODE, cytoHubba and Centiscape2.2 in Cytoscape3.8.2 (Table 2). Isoform_16555 (Anthesis Pomoting Factor 1, MsAPF1), isoform_15622 (DUF1644 domain-containing protein, MsSIZ1), isoform_210768 (Acyl-CoA N-acyltransferase protein, MsGNAT6), isoform_13861 (Detoxification 21, MsHMP21), isoform_16567 (Mitogen-activated protein kinase 10, MsCXIP4), and isoform_107196 (Unknown) were identified as hub genes in lightsteelblue1 module, indicating that these genes have significant functions in Amino acids and nucleic acids metabolism under shade. In paleturquoise module, isoform_10150 (Beta-glucosidase 18, MsBGL18), isoform_92874 (Basic 7S globulin, MsBg7S), isoform_192429 (Cytochrome P450 710A11, MsCYP710A11), isoform_238198 (Transcription factor TGA2.2, MsTGA2), isoform_55976 and isoform_152869 (Pathogenesis-related protein P2, MsPR4) might play central roles in carbon fixation related pathways under shade.
Table 2
the hub genes detected in three modules
Modules
|
Gene ID
|
Gene name
|
Arabidopsis Orthologs
|
Predicted functions
|
Lightsteelblue1
|
isoform_16555
|
MsAPF1
|
AT5G66240.1
|
Anthesis pomoting factor 1
|
isoform_15622
|
MsSIZ1
|
AT3G25910.1
|
DUF1644 domain-containing protein
|
isoform_210768
|
MsGNAT6
|
AT2G06025.5
|
Acyl-CoA N-acyltransferase protein
|
isoform_13861
|
MsHMP21
|
AT1G33110.1
|
Detoxifiction 21
|
isoform_114887
|
MsMAPK10
|
AT4G21970.1
|
Mitogen-activated protein kinase 10
|
isoform_16567
|
MsCXIP4
|
AT2G28910.2
|
Zinc finger protein
|
isoform_107196
|
unknown
|
AT1G23440.1
|
Unknown
|
paleturquoise
|
isoform_10150
|
MsBGL18
|
AT1G61820.1
|
Beta-glucosidase 18
|
isoform_92874
|
MsBg7S
|
AT1G03220.1
|
Basic 7S globulin
|
isoform_192429
|
MsCYP710A11
|
AT2G34500.1
|
Cytochrome P450 710A11
|
isoform_238198
|
MsTGA2
|
AT1G68640.1
|
Transcription factor TGA2.2
|
isoform_55976
|
unknown
|
AT5G57123.1
|
Unknown
|
isoform_152869
|
MsPR4
|
AT3G04720.1
|
Pathogenesis-related protein P2
|
darkolivegreen
|
isoform_10052
|
MsFLA15
|
AT3G52370.1
|
FAS1 domain-containing protein
|
isoform_12760
|
MsUGT73C7
|
AT3G53160.1
|
UDP-rhamnose: rhamnosyltransferase 1
|
isoform_121918
|
MsUGT91C1
|
AT5G49690
|
UDP-rhamnose: rhamnosyltransferase 1
|
isoform_108687
|
MsSBT3
|
AT1G66220.1
|
Subtilisin-like protein protease SBT3.3
|
isoform_11567
|
MsFLA17
|
AT5G06390.1
|
FAS1 domain-containing protein
|
isoform_28177
|
MsLECRK-V.1
|
AT1G70110.1
|
L-type lectin-domain-containing protein
|
isoform_119775
|
MsFMN
|
AT4G27270.1
|
NADPH-dependent FMN reductase
|
isoform_47022
|
MsGHL
|
AT4G31500.1
|
Geraniol 8-hydroxylase-like protein
|
In darkolivegreen module, isoform_10052 (FAS1 domain-containing protein, MsFLA15), isoform_12760 (UDP-rhamnose:rhamnosyltransferase1, MsUGT73C7), isoform_121918 (UDP-rhamnose:rhamnosyltransferase1, MsUGT71C1), isoform_108687 (Subtilisin-like protein protease SBT3.3, MsSBT3), isoform_11567 (FAS1 domain-containing protein, MsFLA17), isoform_28177 (L-type lectin-domain-containing protein, MsLECRK-V.1), isoform_119775 (NADPH-dependent FMN reductase, MsFMN) and isoform_47022 (Geraniol 8-hydroxylase-like protein, MsGHL) were identified as hub genes, suggesting that these genes exert vital functions in photosynthesis in M. sinostellata in response to shade (Table S3).
Verification of hub genes in M. sinostellata under shade by RT-qPCR
We executed a quantitative reverse-transcription PCR (RT-qPCR) analysis to verify the expression level of 21 hub genes under shade treatment after 0d, 5d and 15d in M. sinostellata. These results suggested that the expression level of all the hub genes were significantly changed during shade treatment, which indicating that these genes were all shade-responsive (Fig. 8). Interestingly, hub genes in lightsteelblue1 and paleturquoise modules were all significantly up regulated during shade treatment, while 8 hub genes in darkolivegreen modules were all down regulated under shade in M. sinostellata.