2.1. Effects of drought stress on the growth of Jerusalem artichoke and the determination of physiological indices
This experiment simulated natural drought to cause stress to two genotypes of Jerusalem artichoke at the seedling stage, and the season in which Jerusalem artichoke was planted in Qinghai was selected for the experiment. The seeds were planted on April 28, and the Jerusalem artichoke plant was approximately 25 cm high in early June. Jerusalem artichoke with the same height and similar growth levels were selected for the experiment, and samples were collected from June 4 to June 10. Figure 1 indicates that, compared with the water treatment (WT), there was no significant change in both genotypes of Jerusalem artichoke plants under drought treatment (DT) at 2 d. QY1 wilted at 4 d, and most of the leaves of the plant withered at 10 d, while QY3 showed no significant wilting at 4 d, and most of the leaves remained normal at 10 d. Moreover, it can be seen from the color changes at the root that at 2 d of stress, QY1 became darker, while the color of QY3 did not become significantly darker until 10 d (Fig 1).
Fig 2a indicates that with the development of drought stress, the water content of DT group decreased to less than 10% in 10 days, while that of the WT group remained basically stable at 35%, and the water content of QY3 decreased faster than that of QY1. In contrast, Fig 2b shows that the content of dry matter in roots of two genotypes of Jerusalem artichoke increased significantly under drought stress. On the 2nd and 4th days after the plants were subjected to drought stress, there was no significant difference between the WT and DT groups. However, on the 10th day, the DT group reached its highest value. The difference between the WT and DT groups became highly significant, and the difference between the two genotypes of Jerusalem artichoke in DT group also became highly significant, with QY3 exhibiting higher levels of drought stress than QY1. The rate of accumulation of dry matter by the DT group was higher than that of the WT group. Next, we measured the contents of protein and chlorophyll. With an increase in treatment time, there was no significant difference in the contents of protein between DT and WT groups, but the DT group exhibited a trend of higher protein contents than those of the WT group. There were also significant differences between the two genotypes of Jerusalem artichoke in DT group on the 2nd and 10th days, with QY3 exhibiting higher trends than QY1. There was no significant difference between the WT and DT groups in their level of chlorophyll content on the 2nd and 4th days, but there was a very significant difference between these groups on the 10th day. QY1 and QY3 exhibited the same type of trend under drought stress. However, the range of change in trend of QY3 was larger than that of QY1 (Fig 2c and d). The results showed that continuous drought stress has a substantial effect on the growth of Jerusalem artichoke.
2.2. Characterization of RNA‑Seq
The leaves of two Jerusalem artichoke varieties were sequenced under normal water and drought conditions. A total of 89,550,751 raw reads were amplified. The sequenced clean reads and clean bases were 89,551,965 and 13,508,637,885, respectively. The average values of Q20 and Q30 were 97.66% and 93.265%, respectively, indicating that the sequencing quality was reliable. After data filtering and quality control, 89,551,965 clean reads were produced with an average of 3,198,284 clean reads per sequencing library and an average GC content of 45.585%. After assembly with Trinity, 309,323 unigenes and 455,817 transcripts were obtained. The average GC content of the unigenes and transcripts was 42.24% and 41.48%, respectively; the average N50 of the unigenes and transcripts was 727 bp and 788 bp, and the mean length of the unigenes and transcripts was 571 bp and 619 bp, respectively (Table 1). After splicing and assembling the clean reads, 455,817 transcripts were obtained with an average length of 619 bp and N50 of 788 bp. The number of transcripts between 200 and 400 bp was the most frequent at 200,512, accounting for 18.3% of the total (Fig 3).
The differentially expressed genes (DEGs) between the two Jerusalem artichoke varieties were screened using the DESeq method for subsequent gene function annotation and screening. The results of the difference data are shown in Figure 3. Under drought stress, genes in the QY1 and QY3 leaves diverged significantly. At the second day of drought stress, there were 5,613 genes differentially expressed in QY1/QY3, of which 4,091 genes were up-regulated, and 1,522 genes were down-regulated. At the fourth day of drought stress, there were 12,985 genes differentially expressed in QY1/QY3, of which 6,346 genes were up-regulated, and 6,639 genes were down-regulated. At the tenth day of drought stress, there were 24,923 genes differentially expressed in QY1/QY3, of which 12,081 genes were up-regulated, and 12,842 genes were down-regulated. With the increase in the trend of stress time, the number of up-regulated and down-regulated unigenes both increased, indicating that drought stress inhibited the expression of some unigenes and induced the expression of a number of other unigenes (Fig 4).
Fragments per kilobase of transcript per million (FPKM) > 0.3 was used as the standard pf expression to calculate the gene expression of the two Jerusalem artichoke varieties under normal water and drought stress, as shown in Fig 5. Under WT treatment, 39,569 unigenes were screened out in QY1, including 16,972 nonshared unigenes. A total of 3,635 unigenes were shared between two groups, 1,865 unigenes among three groups and 2,433 unigenes among four groups. In QY1, 9,601, 9,374, 8,632 and 11,962 unigenes were screened at 0 d, 2 d, 4 d and 10 d, respectively. The overall trend of expression of QY1 under the DT treatment was similar to that of WT group, but the number of unigenes increased significantly. A total of 43,650 unigenes were screened out, including 18,866 nonshared unigenes, 3,927 shared unigenes between two groups, 2,066 shared unigenes among three groups and 2,683 shared unigenes among four groups. In QY1, 9,601, 9,462, 11,132 and 13,485 unigenes were screened at 0 d, 2 d, 4 d and 10 d, respectively. More genes were expressed in the DT group than the WT. In QY3, there were more unigenes screened out under the WT and DT treatments than those in QY1. Among them, there were 44,666 unigenes under WT, including 19,091 in nonshared unigenes, 4,129 shared unigenes between two groups, 2,191 shared unigenes among three groups and 2,686 shared unigenes among four groups. In QY3, 13,031, 9,723, 11,084 and 10,828 unigenes were screened out at 0 d, 2 d, 4 d and 10 d, respectively.
According to the unigene detection results under drought stress, GO function classification was conducted for the DEGs. Since the difference was the most obvious at 10 d, the DEG of the two varieties at 10 d of the WT and DT treatments were selected for further functional analysis. After 10 d of drought stress in QY1, 9,394 DEGs were annotated with 43 functions, among which 19 belonged to biological process, 13 to cellular component and 11 to molecular function. After 10 d of drought stress in QY3, 10,796 DEGs were annotated with 41 functions, among which 18 belonged to biological process, 13 to cellular component and 11 to molecular function. The three pathways with the highest differential expression in biological process, cellular component and molecular function were analyzed. In biological process, the differences primarily existed in the metabolic process, cellular process and single organism process, with the DEGs in QY1 comprised of 561, 493 and 451 more than those in QY3. A similar situation also occurred in cellular component and mobile function. In the cellular component, most of the different pathways were membrane, cell part and cell, with the DEGs in QY1 numbering 294, 451 and 268 more than those in QY3, respectively. In molecular function, QY1 had 553, 299 and 64 more DEGs in the quantitative activity, binding and transporter activity channels, respectively. QY1 had more DEGs than QY3 in most pathways, but in general, there were more up-regulated genes than down-regulated genes in QY3. According to the data of differential expression, drought stress has a significant impact on the biological functions described above (Fig 6).
To further understand the metabolic pathway and other processes of drought stress genes in Jerusalem artichoke, MapMan was used to compare the metabolic pathways of differential genes of QY1 and QY3 under the DT and WT treatments based on the results of RNA-Seq. A visual analysis found that there were significant differences in the expression of drought-sensitive genes in the metabolic pathways. As shown in Figure 7, the metabolic pathway was divided into minor CHO metabolism, cell wall metabolism, lipids metabolism, secondary metabolism, amino acids metabolism, light reaction metabolism, carbohydrate metabolism, nucleus metabolism and tetrapyrrole metabolism. Through the comparison of the QY1 and QY3 DT/WT, we can observe that there was a large difference in the overall metabolic pathway between QY1 and QY3. The most different metabolic pathways were secondary metabolism, light reaction metabolism and cell wall. Compared with QY1 (Fig 8a.), the up-regulated genes of QY3 were primarily cell wall, secondary metabolism and light reaction metabolism. The number of up-regulated genes was 130, 98 and 83, respectively, and the down-regulated genes were also primarily concentrated in these three pathways. Secondary metabolism helps plants to resist environmental stress, improve their self-preservation and ability to compete and maintain their survival and development. Therefore, the study of this area is highly significant for the study of the secondary metabolic accumulation of industrial cash crops and their production. Drought, as one of the primary environmental stresses, affects the accumulation of secondary metabolites. We used the secondary metabolic pathway with the largest number of differential genes for further mining. In Fig 7c and d, we continue to use red for up-regulation and blue for down-regulation. The up-regulated genes in QY3 were significantly more abundant than those in QY1 and were primarily concentrated in flavonoids, phenylpropanoids, and the shikimate and terpenoids pathway. The results described above show that under the stress of continuous drought, two types of Jerusalem artichoke exhibited different drought resistance mechanisms, and the synthesis of cell wall, lipid, photosynthesis and secondary metabolism of QY3 were significantly higher than that of QY1. The primary component of the plant cell wall is cellulose, and the primary component of the intercellular layer of cell wall is pectin. Waxy synthesis, cellulose synthesis and the up-regulated expression of pectin methylesterase all indicated the enhancement of cell wall synthesis. QY1 has more up-regulated and down-regulated genes in carbohydrate metabolism than does QY3.
Endogenous hormones, as response factors to plant stress, play an important role in the adaptation to drought stress. ABA, as one of the primary endogenous hormones, plays an important role in plant stress-resistant growth, particularly under drought stress, and is even considered to be a potential indicator of crop drought resistance. In this study, 14 ABA signaling, 21 NAC, 7 DREB and 24 WRKY genes were selected for expression analysis. Some of them were both up-regulated in QY1 and QY3, including HAB2 (PP2C member of ABA signaling), ABA responsive element binding protein AREB1, NAC29 (NAC family), DREB2a (DREB family) and eight WRKY transcription factors: WRKY13, WRKY5, WRKY7, WRKY9, WRKY42, WRKY4, WRKY68, and WRKY8 (Fig 8.). NAC1 and some WRKY genes were down-regulated in both varieties of Jerusalem artichoke. On the whole, QY3 showed more variation in the NAC and WRKY genes than did QY1, particularly in the NAC family. In addition, QY3 showed an obvious trend of change in three SNRK genes, but the number of genes responding to drought stress was significantly higher in QY1.
To better understand the drought-tolerant gene expression regulatory network of Jerusalem artichoke, weighted gene co-expression network analysis (WGCNA) was used to analyze the drought-related genes. Based on the results of high mean connectivity and low scale independence, we finally chose β=14 as the threshold. The analytical results showed that 6,643 drought-related genes of QY1 could be clustered into 27 modules under drought stress, while 7,224 drought-related genes in QY3 could be clustered into 28 modules (Fig 9.). The WGCNA analysis showed that under drought stress, re-writing and de-writing existed in the co-expression network of a large number of genes of QY1 and QY3. We selected the dark slate blue (r = 0.45, p = 0.02), dark red (r = 0.38, p = 0.04), ivory (r = 0.42, p = 0.03) and brown (r = -0.39, p = 0.04) modules with the highest correlation with stress to perform a separate comparative analysis. The dark slate blue module had a significant correlation of expression in variety, stress and period traits (Fig 9d.). In the module of the trend of expression diagram that was drawn, QY1 and QY3 showed a consistent trend of expression in the dark slate blue module under the WT, but they showed an opposite trend of expression under the DT treatment. In addition, QY1 and QY3 had the same trend of expression in the ivory, dark red and brown modules in the DT group.
The functions of the genes in four stress resistance-related co-expression modules were analyzed. The GO enrichment results showed that 319 genes in the brown module were enriched in the items of molecular function (GO:0003674), catalytic activity (GO:0003824) and metabolic process (GO:0008152). A total of 130 genes in the dark red module were enriched in the items of cell component (go: 0005575), cell (go: 0005623) and cell part (go: 0044464). Four genes in the dark slate blue module were enriched in transferase activity (go: 0016740), and 24 genes in the ivory module were enriched in response to hormone stimulus (go: 0009725), response to endogenous stimulus (go: 0009719), response to organic substance (go: 0010033) and response to chemical stimulus (go: 0042221). In general, the genes enriched from the brown and dark red modules had the greatest number of molecular functions, while the primary GO items enriched from brown and ivory modules did not contain a molecular function (Fig 10.). To further explore the candidate genes involved in the drought stress of Jerusalem artichoke, the candidate genes with the highest correlation with drought stress were screened out from the GO classification of four modules. Through the function of gene annotation, it was found that the 16 genes selected by annotation belong to histone, ABA regulation and protein kinase among others. From the annotated species, there were eight candidate genes annotated to Cynara cardunculus, indicating that the mechanism of drought stress of Jerusalem artichoke is similar to that of C. cardunculus (Table 2).