Construction of drought stress regulation networks in potato based on SMRT and RNA sequencing data

Plant materials and drought stress treatments

A drought tolerant tetraploid potato “Désirée ”[45, 49], provided by Southwest University, was used in this study. In 2 weeks old vitro plantlets of about 10 cm were transferred and stablished into hydroponic solution (1/4 Hoagland nutrient solution) for 1 week at 22/18 °C (16 h light / 8 h dark), then the materials were divided into two parts, one for control and the other for drought stress treatment to simulate the drought stress as reported by different authors, which were treated with 20% polyethylene glycol (PEG)-6000 [25, 48, 76]. The aboveground tissues were taken at 0, 1, 3, 6, 12, 24 and 48 hours after drought stress treatment and immediately frozen at least 30 min in liquid nitrogen and kept at − 80 °C until use. The samples were collected for three biological replicates. The collected samples were sent to Biomarker for full-field transcriptome sequencing and second-generation transcriptome sequencing, including sample detection, library construction and sequencing.

Library construction for Illumina and PacBio sequencing

Total RNA of all samples was extracted using EZ-10 DNAaway RNA Mini-Preps Kit (Sangon Biotech, Shanghai, China) and the quality of total RNAs was evaluated using an Agilent 2100 Bio-analyzer (SA Pathology, Adelaide, SA, Australia). High-quality RNAs (RIN > 8.0) were used for the following cDNA synthesis and library construction. For Illumina RNA sequencing, 1.0 μg RNA per sample was used for cDNA synthesis with the NEBNext® Ultra™ II RNA Library Prep Kit for Illumina (NEB, USA) and sent to Biomarker Technologies Corporation (Beijing, China) for Illumina sequencing. For SMRT sequencing, equal amount of the total RNA of all samples were pooled together as the template for cDNA synthesis with the SMARTer PCR cDNA Synthesis Kit (Clontech, USA) in Biomarker Technologies Corporation (Beijing, China).

Analysis of the Illumina and SMRT sequencing data

The raw Illumina data was filtered to generate clean reads using FastQC tool (http://www.bioinformatics.babraham.ac.uk/projects /fastqc/) (v0.11.5) and remove the adaptor contamination and low-quality reads using Trimmomatic (v0.35) [7]. The high-quality clean reads were then mapped to the reference genome database (http://solanaceae.plantbiology.msu.edu/rh_potato_download.shtml) using HISAT2 version 2.1.0 and the default parameters [33]. For SMRT sequencing data, the circular consensus sequences (CCS) were generated from the PacBio Iso-Seq raw reads with smrtanalysis_2.3.0 (http://www.pacb.com/) and then processed into error corrected reads of insert (ROIs) with a minimum full pass of > 3 and minimum accuracy of > 0.9. These ROIs were classified as full-length, non-chemiric (FLNC) transcripts using smrtanalysis_2.3.0 based on the presence and location of the 5′ and 3′ cDNA primers and poly (A) tail. All full-length reads were aligned to the S. tuberosum group Tuberosum RH89–039-16 reference genome (http://solanaceae.plantbiology.msu.edu/rh_potato_download.shtml) using GMAP software [78]. Finally, all mapped reads were collapsed using a Python script in the ToFU package (https://github.com/PacificBiosciences/cDNA_primer/) with min-coverage = 85% and min-identity = 90%. All raw data was deposited in the National Center for Biotechnology Information (NCBI) sequence read archive with accession no. PRJNA728834.

Gene functional annotation, enrichment analysis, and weighted gene co-expression network (WGCNA) analysis

Gene functions was predicted using the following databases: NR (NCBI non-redundant protein sequences), Pfam (protein family), Swiss-Prot, KEGG (Kyoto Encyclopedia of Genes and Genomes), and GO (gene ontology). The GOseq package in R (version 3.4.3) and the KOBAS software were used to perform the GO and KEGG pathway enrichment, respectively. The gene expression levels were calculated using FPKM (fragments per kilobase of exon per million fragments mapped). Genes with FPKM values > 1 in each sample were used for WGCNA analysis using the R package WGCNA (v1.42).

Differential expression analysis and candidate gene determination

Differentially expressed genes (DEGs) were analyzed by DESeq and detected with the parameter FDR < 0.01 & |Log2Ratio| > 1. The plant transcription factors were predicted by iTAK software, and the differentially expressed transcription factors were screened by mapping to the DEGs database and the genes related to plant hormones in DEGs were screened based on Gene Ontology (GO) functional annotation. Twenty-nine DEGs were randomly selected and verified using the real-time (RT)-PCR to detect the accuracy of transcriptome data (Supplementary Table S1). Toolbox for Biologists v1.087 [9] was used to construct the heatmap.

Identification of lncRNAs and AS events

In order to analyze the coding potential in the transcripts, four of the most widely used computational approaches including coding potential calculator (CPC) (Li, [85]), coding-non-coding Index (CNCI) [5], coding potential assessment tool (CPAT) [13], and pfam protein structure domain analysis were combined. Transcripts were regarded as lncRNA candidates using the parameter: lengths > 200 nt & > 2 exons.

To identify the AS (Alternative splicing) events, the ASTALAVISTA algorithm [17] based on the assembled GTF file from the PacBio Iso-Seq and Illumina data was conducted and five types of AS events including IR (Intron retention), ES (Exon skipping), A5 (Alternative 5′ splice site), A3 (Alternative 3′ splice site) and MEE (Mutually exclusive exon) were classified with the default settings.

Isoform identification and reverse transcription (RT)-PCR

To distinguish AS events, the reverse transcription (RT)-PCR was performed. Briefly, the RNA used for Iso-Seq was isolated using EZ-10 DNAaway RNA Mini-Preps Kit (Sangon, https://www.sangon.com) and synthesized to complementary DNA using the Hifair® II 1st Strand cDNA Synthesis SuperMix for qPCR gDNA digester plus (Yeasen). Isoform-specific primers were designed using Vector NTI software (Supplementary Table S1). PCR amplification was performed using 2 × Taq PCR Master Mix (K1034, https://www.apexbt.com/) and the PCR products were showed in agarose gel stained with GoodView™ Nucleic Acid Stain (Apr-13-2021, HGV-II, https://www.sbsbio.com).

BlazeTaq™ SYBR® Green qPCR Mix 2.0 (GeneCopoeia, USA) was used for RT-qPCR reactions on a CFX96 Real-time System (BIO-RAD, USA). The amplification program was briefly described in previous study [73]. To ensure the precise and reproducible results, each sample was replicated three times and Ef1α in potato were used for each sample as an endogenous control (Supplementary Table S1). The formula F = 2 ^–ΔΔCt was used to calculate the relative expression levels of selected genes.

Alternative polyadenylation analysis

Alternative splicing and alternative polyadenylation (APA) of pre-mRNAs greatly contribute to transcriptome diversity, coding capacity of a genome and gene regulatory mechanisms in eukaryotes. TADIS pipeline (Transcriptome Analysis Pipeline for Isoform Sequencing) was used to further analyze the full-length non-chimeric sequence (FLNC) to identify APA [1] and MEME analysis was also performed to identify poly (A) sequence signals in our data.

Illumina and PacBio sequencing data analysis

To construct a comprehensive regulation network in response to drought stress, control and treated samples were processed through SMRT sequencing on the PacBio platform and next-generation sequencing (NGS) on the Illumina Hi-Seq platform. High-quality RNA was extracted from above tissues of “Désirée” after PEG-6000 treatments (0 h, 1 h, 3 h, 6 h, 12 h, 24 h and 48 h) for Illumina Hi-Seq and pool together for construction of Iso-Seq library. After removing low-quality reads and adapter sequences, 38,347,790-46,349,352 reads were obtained and more than 80% sequences were mapped to the reference genome (Table 1). The Pearson’s correlation coefficient showed that all correlation values between the three replicates ranged from 0.826 to 0.925, indicating a perfect positive correlation and that the sequencing results could be used for further analysis (Fig. 1). For Iso-Seq, the PacBio platform generated 33.24 Gb clean data. Among these data, 464,332 circular consensus (CCS) reads were obtained using Iso-seq pipeline with min full pass > 3 and min predicted accuracy > 0.9. Then, 375,387 full length non chimeric (FLNC) sequences were determined by searching for the polyA tail signal and the 5′ and 3′ cDNA primers and 121,434 high-quality consensus sequences were obtained by clustering the FLNC sequences. SMRT sequencing can offer better transcript integrity (up to 50 kb in length) than NGS and has greatly contributed in maize [74], populous [8] oilseed rape [84] and sorghum [1, 44], but no literature has been reported in potato. Significantly longer isoforms were detected by SMRT than by NGS assembly, demonstrating better isoform integrity (Fig. 2a). In addition, the overall expression levels of transcripts detected by SMRT were higher than those detected by NGS, indicating better quantification accuracy (Fig. 2b).

Analysis of gene expression and transcriptomic changes in response to drought treatment

DESeq was used for differential expression analysis between control versus treatments. False Discovery Rate (FDR) was used as statistical parameter for differential expression gene screening. Fold Change≥2 and FDR < 0.01 were used as screening criteria. In total, 31,637 DEGs were detected between control versus treatments, and the number of DEGs in D24h sample was the largest. (Fig. 3a and Supplementary Table S2). Moreover, almost all of DEGs in the libraries showed moderate expression levels with only a small percentage of the genes expressed in high levels (Fig. 3b and Supplementary Table S3). More than 88% of the DEGs had a − 3<fold change ≤3. Transcription factors (TFs) play important roles in plant abiotic stress responses (Fig. 3c) [29]. To further analyze the function of the DEGs, we used iTAK software to predict transcription factors in DEGs. In total, 664 DEGs belonging to 28 TF families were identified. The TF family members were unevenly distributed. Seven TF families, bHLH (69), MYB (60), ERF (53), bZIP (46), NAC (43) WRKY (35) and HD-ZIP (33) accounted for more than 51% (339/664) of all TFs (Fig. 3d). In total, 919 DEGs belonging to eight classifications of plant hormones, and ABA, IAA and CK-related genes accounted for the majority, accounting for about 94% of all hormone-related genes (Fig. 3e). To combine the hormone profiling data with the expression profiles of genes involved in hormone metabolism and signaling, we analyzed 231 genes involved in plant hormone metabolism and visualized them as heat maps after functional classification (Fig. 4). Heat map results showed that DEGs were most involved in plant hormone signaling.

Gene profiling in response to drought stress in potato

To explore drought stress response mechanisms in potato, the expression patterns were calculated (Fig. 5a). The hierarchical cluster analysis was done and then the K means algorithm used for 98,624 genes. Fourteen clusters were obtained and 63 pathways were enriched, such as carbon metabolism, carbon fixation in photosynthetic organisms, porphyrin and chlorophyll metabolism, protein processing in endoplasmic reticulum and oxidative phosphorylation (Fig. 5b and c). Furthermore, WGCNA was also performed to screen genes response to drought stress by detecting molecular modules. Thirteen modules were generated and more than one module was significantly correlated with each sample. Blue, darkred, midnightblue, orange, cyan, lightgreen, and darkgreen were significantly correlated with DCK, D1h, D3h, D6h, D12h, D24h and D48h, respectively (Fig. 5d). Then, KEGG enrichment analysis was conducted to explore the enriched pathways in each module correlated to each sample. Galactose metabolism (ko00052) were enriched in more than two treated samples. Plant-pathogen interaction (ko04626) and glutathione metabolism (ko00480) were enriched in early stages after drought stress (Fig. 5e and Fig. S1). In addition, through the comparative analysis of the control group and the other six samples processed in different periods, a total of 12,798 differentially expressed genes after deduplication were found (Fig. S2 and Supplementary Table S4). The functional enrichment of these DEGs during the six drought treatment period was conducted (Fig. 6a) and 29 DEGs were selected to confirm the transcriptome sequencing data using the real-time (RT)-PCR (Fig. 6b).

Putative roles of AS for drought stress regulation in potato

AS plays critical roles in abiotic stress response by increasing genetic diversity in plants [36, 42]. However, little information about AS involved in drought stress response has been obtained in potato. In this study, 21,756 AS events were categorized into five patterns, including 12,453 retained intron (RI), 2493 skipping exon (SE), 236 mutually exclusive exons (MX), 4481 alternative 3′ splicing sites (A3), and 2093 alternative 5′ splicing sites (A5) using Isoform sequencing technology, which yields long reads without assembly (Fig. 7a). There is similar number of AS events for control and drought treatments, among these AS events, RI events predominated, accounting for 57.2%, followed by A3, SE, A5 and MX (Fig. 7b), which was consistent with findings of previous reports in other plant species [83]. Four random selected genes with AS events were successfully validated by reverse transcription (RT)-PCR (Fig. 7c).

To identify the co-expression network of genes involved in AS events, the WGCNA package was used and three significant modules were identified (brown, blue and turquoise). The DCK sample correlated to the brown module (r = 0.79, P = 2e^− 5), whereas the turquoise (r = 0.71, P = 3e^− 4) and blue (r = 0.77, P = 5e^− 5) modules were significantly correlated with samples after 6 and 24 h treatments, respectively (Fig. 7d). KEGG enrichment of the significant modules was also performed. Six, eleven and four pathways were significantly enriched in brown, turquoise and blue modules. The galactose metabolism (ko00052) was significantly enriched in turquoise and blue modules, which were treated with PEG-6000 after 6 and 24 h, respectively (Fig. 7e).

Putative roles of lncRNA involved in drought stress regulation in potato

LncRNA is a major component in transcriptome sequencing, and plays important roles in the regulation process after transcription [34, 39, 50, 58]. However, few information about lncRNA in potato has been obtained. In this study, 3345 lncRNAs with various lengths were identified using composite computational algorithm based on the PacBio sequencing data (Fig. 8a, b). More than 54% of lncRNAs were located in the intergenic region while less than 7% were located in the antisense region according to reference annotations (Fig. 8c). To identify the co-expression network of target genes of lncRNAs, the WGCNA package was used and yellow, brown, blue and turquoise modules were significantly enriched (Fig. 8d, Supplementary Table S5). The lncRNA target genes of KEGG analysis indicated that terpenoid backbone biosynthesis (ko00900) and phenylalanine metabolism (ko00360) pathway were most enriched in DCK and D6h; Glycerophospholipid metabolism (ko00564) and ubiquitin mediated proteolysis (ko04120) pathway were most enriched in D12h and D24h, respectively (Fig. 9).

Alternative polyadenylation analysis

In eukaryotes, polyadenylation is an important transcriptional modification in mRNAs [15, 61]. Alternative polyadenylation can easily increase the transcripts variety and regulate gene expression by producing different transcripts [18, 60, 62]. The iso-seq data allowed us to detect polyadenylation sites in potato using the program TAPIS [1]. In this study, 26,153 poly (A) sites from 13,010 genes detected in the Iso-Seq data were identified (Supplementary Table S6). On average, 2.01 poly (A) sites per gene were found and 856 genes had at least five poly (A) sites (Fig. 10a). In addition, a clear nucleotide bias of the flanking region of all poly (A) sites was found. To confirm this phenomenon, we investigated the motifs of 50 nucleotides upstream of predominant sites using MEME analysis. Furthermore, three conserved polyadenylation signals, i.e. TTTTGT, GCCCCC and GGGGCG were overrepresented in upstream of poly (A) cleavage sites (Fig. 10b).

DEGs involved in drought stress response

In this study, SMRT and NGS technologies were combined to conduct a global survey of drought-treatment potato transcriptome. With the intensification of drought stress, the number of DEGs increased. The number reached the highest at 24 h after drought stress treatment (Fig. S1 and Supplementary Table S2). which was consistent with the findings previously reported in potato [10]. Plant hormones, such abscisic acid (ABA), auxin (IAA), and cytokinin (CK), play vital roles in abiotic stresses [72]. Transcription factors and plant hormone-related genes in differentially expressed genes were analyzed, and found that drought stress-related transcription factors such as bZIP and plant hormone-related genes such as ABA-related genes may played a dominant role. The functional annotation of differentially expressed genes confirmed their important role in drought stress, and 29 DEGs we verified using real-time (RT)-PCR (Supplementary Table S1 and Fig. 6b). Moreover, the WGCNA screening and KEGG enrichment analysis of drought stress response genes suggested that plant-pathogen interactions and glutathione metabolism pathways were enriched in the early stage of drought stress, and galactose metabolism pathways were enriched in at least two treatment samples, which indicated these pathways play an important role in the potato drought response process. Previous researches reported that galactose is involved in response to drought, high salt and low temperature tolerance [70], and transgenic studies demonstrated that glutathione helped Tortula ruralis adapt to rehydration following rapid desiccation [14]. It was clear that the number of genes down regulated was greater that the upregulated which means that the plant stops regular active metabolism as main mechanism of response to the drought stress.

AS events involved in drought stress response

One of the advantages of PacBio sequencing is its ability to detect AS events. (Wang, [75, 87]). For different biological responses, AS is an effective mechanism to increase the complexity and flexibility of the transcriptome [41]. Previous studies have confirmed that AS events play critical roles in potato growth and abiotic stress response [43, 52]. AS events identified during the potato drought stress response process were analyzed in this study. The proportion of different AS events was consistent with the findings previously reported in other plant species [80, 83]. Our experiment verified the AS events detected above, and the results showed that the size of the gel band fragments was consistent with the splice isoforms identified from the Iso-Seq data (Supplementary Table S1 and Fig. 7c) and the functional annotations of these genes also indicated that they were involved in drought stress response. In summary, we have verified that the four randomly selected genes undergo alternative splicing under drought stress, and their functional annotations are also related to plant stress tolerance. Our experiments results and sequencing data indicated that potatoes may increase their protein diversity under drought stress through alternative splicing, and thereby improve their environmental adaptability.

LncRNA plays an important role in drought stress response

As a type of non-coding RNA, lncRNA act as regulators in a series of biological processes [26]. Previous research suggested that lncRNA has potential functions in potato defense response and development [34, 56]. In this study, a total of 3445 lncRNAs of different lengths were identified through four methods: CPC analysis, CNCI analysis, pfam protein domain analysis, and CPAT analysis (Fig. 8a and b). The WGCNA and KEGG analysis results of lncRNA target genes show the term ‘phenylalanine metabolism’ was enriched in the brown module, indicating that lncRNAs respond to drought stress by targeting genes related to amino acid metabolic pathways (Fig. 9), it also provided additional effective candidate genes for future functional characterization.

APA is essential for drought resistance of potato

Polyadenylation refers to the covalent linkage of polyadenylic acid to messenger RNA (mRNA) molecules. In the process of protein biosynthesis, this is part of the way to produce mature mRNA ready for translation. In eukaryotes, polyadenylation is a mechanism that interrupts mRNA molecules at their 3’ends. The polyadenylic acid tail (or poly A tail) protects mRNA from exonuclease attack, and is very important for the termination of transcription, export of mRNA from the nucleus, and translation. A comprehensive genome-wide APA map draft consisting of 26,153 poly (A) sites from 13,030 genes was constructed in this study (Supplementary Table S6). The nucleotide bias was evaluated and obvious nucleotide bias in the flanking regions of all poly (A) sites were found (Fig. 9b). These findings are consistent with those of previous studies in Brassica napus [84]. Two conserved polyadenylation motifs in B. napus were identified, while three conserved polyadenylation motifs in potato were identified, indicating that the polyadenylation motifs are species-specific [40]. Previous studies have shown that alternative polyadenylation of RNA is critical for gene function by increasing the complexity of the transcriptome and regulating gene expression [79].

Ethics approval and consent to participate

The potato material used in this study has been approved by Southwest University, we comply with the IUCN Policy Statement on Research Involving Species at Risk of Extinction and the Convention on the Trade in Endangered Species of Wild Fauna and Flora, All methods for this study are in compliance with the provisions of Chinese law and international guidelines and regulations.

Consent for publication

Not applicable.

Competing interests

The authors declared no conflict of interest in the authorship and publication of this document.