Oligo-dT anchored cDNA-SRAP and cDNA-SCoT aided identification of transcripts differentially expressed during the early stages of recovery of resurrection plant Haberlea rhodopensis Friv. from freezing-induced desiccation

Abstract

Resurrection plants are a unique group of flora showing an extraordinary ability to survive extreme water shortages. Among them, the Bulgarian endemic plant Haberlea rhodopensis Friv. has been used as a model to study the mechanisms underlying desiccation tolerance induced by water and low temperature stress. Although extensive transcriptome data, including two RNAseq experiments, are available for drought induced desiccation and recovery, the gene expression during recovery after freezing-induced desiccation (RAF) remains largely unexplored. In this study we performed a differential screening using a new approach, based on oligo-dT anchored cDNA-SRAP (copy DNA sequence-related amplified polymorphism) and a modified oligo-dT anchored cDNA-SCoT (copy DNA start codon targeted) assays to provide new information on genes involved in the tolerance to low temperature stress and the recovery of H. rhodopensis. The putative functions of 19 identified transcript derived fragments were established by searching for homology with known sequences in public databases and the expression of 10 of them during the early stages of recovery was confirmed by quantitative reverse transcription polymerase chain reaction. Some of the identified genes have known functions in the response to abiotic stress in other plant species including resurrection plants, which confirms the effectiveness of the approach used in the present study. The results of this study will contribute to obtaining new information on the defence mechanisms in the recovery process of H. rhodopensis, especially at the early stages of recovery, which is the most vulnerable period for plants.

Keywords:

• Resurrection plants
• differential screening
• recovery after freezing-induced desiccation
• cDNA-SRAP
• cDNA-SCoT
• qRT-PCR

Introduction

The unique ability of resurrection plants to survive drought to an air-dry state makes them an ideal model system for studying the mechanisms of their tolerance to extreme drought conditions. Haberlea rhodopensis (Gesneriaceae) is a homoiochlorophyllous resurrection angiosperm species from the Balkan Peninsula, which is characterized with high ecological plasticity. Unlike many of the resurrection plants, European gesneriads including H. rhodopensis are exposed to low negative temperatures in winter, which also causes desiccation of plants, and they survive the harsh winter conditions in a dry state [1–3]. The protective strategies used by resurrection plants to cope with the detrimental effect of low temperatures are cell ultrastructural rearrangement, lowered glass transition temperatures, downregulation of photosynthesis, enhancement of the antioxidant system, changes in membrane lipid composition, accumulation of different protective sugars and proteins [1,3–5].

In order to better understand the survival strategy of H. rhodopensis it is important to study the protective mechanisms both during drought and in the process of plant recovery after rehydration. The rapid water uptake in the early hours of rehydration can cause cellular damage due to reactive oxygen species (ROS) generation [6]. The successful recovery of resurrection plants in response to drought- or freezing-induced desiccation depends on the involvement of protective mechanisms during dehydration of plants [4,6]. Our previous investigation of the physiological and biochemical alterations after rehydration from freezing-induced desiccation of H. rhodopensis showed that the water uptake during the first 15 h of recovery was slow, cell and chloroplast ultrastructure were restored, PSI activity recovered faster than PSII, the content of photosynthetic proteins related to thermal energy dissipation was elevated, alternative electron flows were involved in the process of recovery, antioxidant levels, proteolytic activity, stress-induced proteins and ELIP expression were kept high [7–10]. Up to date, the molecular mechanisms used by resurrection plants and H. rhodopensis in particular, to limit the ROS-induced damage during the early hours of rehydration after freezing-induced desiccation are scarcely investigated [9].

The differential display method, developed by Liang and Pardee [11] more than 20 years ago, is still a powerful approach for the identification of differentially expressed genes in a large number of cellular organisms and many modifications of this method have been adapted to specific cases. Sequence-related amplified polymorphism (SRAP) and start codon targeted (SCoT) methods were initially developed as simple and reliable gene-tagged molecular marker systems [12,13] and have been extensively used for cultivar identification and gene diversity studies in a number of crops as well as natural populations of herbs and medicinal plants [12–17]. The copy DNA (cDNA)-SRAP and the single primer cDNA-SCoT methods have also been employed to identify differentially expressed genes in a large number of plants [18–21].

The aim of this study was to identify new genes involved in the mechanisms of H. rhodopensis recovery from freezing-induced desiccation as a part of the survival strategy of this resurrection plant by harnessing novel differential screening strategies based on oligo-dT anchored cDNA-SRAP and a modified oligo-dT anchored cDNA-SCoT assays.

Materials and methods

Desiccation and rehydration of plants

Haberlea rhodopensis Friv. tufts were initially collected from the Rhodope Mountains and further cultivated under ex situ environmental conditions. Freezing-induced desiccation and rehydration of plants were performed as described by Georgieva et al. [7]. In brief, four tufts (4–6 plants in each) were kept under natural environmental conditions continuously thus allowing acclimation to low positive temperatures and exposure to freezing temperatures in the winter. The dehydration of the plants began at freezing temperatures below −10 °C, and they overwintered in an air-dry state. The rehydration of plants after freezing-induced desiccation [recovery after freezing (RAF)] was carried out in laboratory conditions at 21–23 °C, 16/8 light/dark cycle and light intensity 25–30 Ï»mol m⁻² s⁻¹. Initially, the soil substrate was well watered and a high constant air humidity was maintained throughout the recovery.

Leaf tissues for RNA extraction were collected from dry leaves (0 h), after 1, 3, 5, 7, 9, 24 h and 7 d of rehydration and immediately frozen in liquid nitrogen, and stored at −70 °C. At each time point, the leaf samples were pooled from 4 to 6 plants. Relative water content (RWC) was determined gravimetrically as described previously [7].

RNA extraction

Total RNA was extracted from 100 mg of leaf tissue, from all collected time points at desiccated state and after the start of watering of RAF plants. The extraction was done by using the Trizol based Direct-zol RNA miniprep kit, Cat. #R2050 (Zymo Research, Irvine, CA, USA) and in a column treated with DNAase I according to manufacturer’s instructions. RNA quantity and quality was determined using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA) and confirmed by ‘Bleach gel’ electrophoresis [22].

cDNA synthesis, oligo-dT anchored cDNA-SRAP and cDNA-SCoT differential display

The first strand cDNA synthesis for oligo-dT anchored cDNA-SRAP and cDNA-SCoT was done with the aid of PrimeScript™ RT reagent Kit with gDNA Eraser, Cat. #RR047A (Takara Bio, Shiga, Japan). The reverse transcription was performed according to the kit manufacturer’s protocol starting from 2 μg total RNA and primed with the oligo-dT18vn_M13_neg40 primer (Supplemental Table S1) to add anchor sequence at the 3′-end of the synthesized cDNAs for the annealing of the M13 (-40) primer used in oligo-dT anchored cDNA-SRAP and cDNA-SCoT assays. The reverse transcription reaction was performed at 37 °C for 15 min followed by deactivation of the RT enzyme for 5 s at 85 °C. The resulting first strand cDNA was stored at −20 °C until use for differential screening.

For polymerase chain reaction (PCR) amplification, the oligo-dT anchored cDNA-SRAP or cDNA-SCoT differential display methods use one SRAP or SCoT primer combined with the oligo-dT anchored primer (M13_neg40_F, Supplemental Table S1). Each 10 μL PCR reaction consisted of 5 μL MyTaq HS 2× PCR mix (Bioline, London, UK), 0.4 μL of each primer (10 μmol/L), and 1 μL of 10-fold diluted cDNA template. All PCR amplifications were performed on a Veriti 96 well Thermal Cycler (Thermo Fisher Scientific, Waltham, MA, USA) under the following cycling conditions: initial denaturation for 3 min at 95 °C; 35 cycles of 15 s at 95 °C, 60 s at 50 °C and 2 min at 72 °C; and final extension for 10 min at 72 °C. For the oligo-dT anchored cDNA-SRAP and cDNA-SCoT differential screening we used three bulked cDNA probes. The probe at freezing desiccated state (D) contained cDNA from dry leaf tissue at 10% RWC (0 h). For the early recovery (ER) bulk probe, equal amounts of cDNAs from 1 and 3 h of recovery were mixed, while the late recovery (LR) bulk combined the cDNAs from 24 h and 7 d of recovery.

Isolation, sequencing and functional identification of transcript derived fragments

The fragments amplified by oligo-dT anchored cDNA-SRAP and cDNA-SCoT methods were electrophoretically separated for 3 h at 50 V in 1.5% ТАЕ agarose gels using MUPID Mini Electrophoresis system (ATTO, Tokyo Japan). The amplification products were visualized on a Syngene™ UV Transilluminator 2020M (Syngene, Cambridge, UK) and documented using EC3 Darkroom W/Gel HR Camera (UVP LLC, Upland, CA, USA) after staining for 20 min in a 3× GelRed™ solution.

The transcript derived fragments (TDFs) identified as differentially expressed were cut off from the gel, purified by using a GeneJET^® Gel Extraction Kit (Thermo Fisher Scientific, Waltham, MA, USA) and after elution in sterile milliQ H₂O were amplified with the same primers and PCR conditions as those used in the oligo-dT anchored cDNA-SRAP or cDNA-SCoT analysis. The resulting products were purified using a GeneJET^® Gel Extraction Kit, according to the kit manual, eluted in 15 μL elution buffer and after concentration adjustment, mixed with sequencing primer and sent for sequencing. The capillary Sanger sequencing was performed by Macrogen-Europe (Amsterdam, NL). The resulting TDF sequences were compared to sequences published in the NCBI databases (NCBI; http://www.ncbi.nlm.nih.gov/BLAST; acc. date: 2 May 2023), using BLASTN and BLASTX algorithm [23,24]. TDF homologous to genes with already established function were classified (BLASTN and/or BLASTX) according to their putative function.

cDNA synthesis and qRT-PCR

The quantitative reverse transcription PCR (qRT-PCR) primers, with the exception of those selected from previously published papers, were designed using the Primer3plus webserver at https://www.primer3plus.com by selecting the server default settings for qPCR primers design (Supplemental Table S1). The designed primers were tested for specificity by the command line tool MFEPRIMER [25] against Haberlea RNAseq contigs [26,27].

To prepare the first-strand cDNA for use in qRT-PCR reactions, we used PrimeScript™ RT reagent Kit with gDNA Eraser, Cat. #RR047A (Takara Bio, Shiga, Japan) and 1 μg total RNA for each sample following the kit provider manual. The reverse transcription reaction was done at 37 °C for 15 min, using the oligo-dT/random primers mix provided with the kit, followed by deactivation of the RT enzyme for 5 s at 85 °C.

The quantitative PCR (qPCR) was carried out in 20-μL reactions with three technical replicates on an ABI 7300 Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA). The reaction mixes consisted of 1 x TB Green Premix Ex Taq (Tli RNase H Plus) (Takara Bio, Shiga, Japan), 2 μmol/L of each forward and reverse primers (Supplemental Table S1), 0.4 μL of 50× Rox dye and 2 μL of 20× diluted first-strand cDNA. The cycling conditions were according to the two-step protocol recommended by the manufacturer, and included initial denaturation for 30 s at 95 °C, followed by 40 cycles of 95 °C for 5 s to 60 °C 31 s. Data acquisitions were performed with the SDS v1.4 software (Thermo Fisher Scientific, Waltham, MA, USA). For relative quantification, the Actin 7 (HrActin, GB: GT270756) was selected as the reference gene.

Data analysis

Data are presented as mean values with standard deviation (±SD) from three technical replicates. Data analysis was performed using Genex Macro, a BioRad MS Excel Gene Expression macro v1.1 (BioRad, Hercules, CА, USA). Relative expression changes were calculated according to the Pfaffl method [28]. For each gene, the relative expression was scaled to the expression at the latest recovery stage (LR or 7 d of RAF) as a calibrator.

Results

Differential screening by oligo-dT anchored cDNA-SRAP and SCoT

The screening for TDFs that are differentially expressed at the early stages of H. rhodopensis recovery from freezing-induced desiccation was initially attempted by cDNA-SRAP with different combinations between ME and EM primers. Although the original cDNA-SRAP method allowed the identification of differentially expressed fragments, many of them were too short to be reliably identified by direct sequencing (data not shown). Therefore, to extend the differentially expressed fragments up to the 3′ poly A tail, we designed a novel oligo-dT anchor primer that adds an M13 (-40) primer sequence at the 3′-end of each cDNA. This allows using M13 (-40) primer in combinations with a single SRAP primer in PCR reactions to screen for differentially expressed TDFs. We named this method ‘oligo-dT anchored cDNA-SRAP’ by analogy with the previously reported oligo-dT anchored cDNA-SCoT method [29].

The oligo-dT anchored cDNA-SRAP method was tested in combinations of the M13 (-40) anchor primer with five ME and six EM SRAP primers (Supplemental Table S1). Among the tested combinations, those with ME primers produced clearer electrophoresis patterns and sharp bands allowing identification and excision of potential differentially expressed TDFs (Figure 1). Oligo-dT anchored cDNA-SCoT differential display was also performed with the newly designed anchor primer in combination with 22 SCoT primers, listed in the Supplemental Table S1, and examples of clear identification of differentially expressed TDFs that were selected for direct sequencing are shown in Figure 2. In both methods, bands with apparent increasing intensity in ER cDNA mix compared to dry condition (D) and with significantly reduced intensity or not present in LR were selected. In total, 29 potential differentially expressed bands, ranging in size from 200 to 1300 bp, were successfully eluted from the agarose gels, reamplified and selected for direct sequencing.

Figure 1. Representative example of oligo-dT anchored cDNA-SRAP analysis with three different SRAP primers (ME3, ME4 and ME5) during RAF: recovery after freezing-induced desiccation in H. rhodopensis. D: freezing desiccated state; ER: early recovery (1 + 3 h); LR: late recovery (24 h + 7 d); M: 100 bp Plus™; arrows indicate the differentially expressed bands that were selected for reamplification and sequencing.

Figure 2. Examples of oligo-dT anchored cDNA-SCoT analysis with five different SCoT primers (SCoT1, SCoT2, SCoT8, SCoT9 and SCoT13) during RAF in H. rhodopensis. D: freezing desiccated state; ER: early recovery (1 + 3 h); LR: late recovery (24 h + 7 d); NTC: PCR reaction without cDNA; M: 100 bp Plus™; arrows indicate the differentially expressed bands that were selected for reamplification and sequencing.

Sequencing, functional identification and annotation of TDFs

Among the 29 TDFs subjected to direct sequencing in both directions, 10 (35%) produced noisy sequencing results, suggesting that the bands contained mixtures of two or more TDF, and those were excluded from further analyses. The remaining 19 TDFs (8 with SRAP and 11 with SCoT primers) produced sequence reads of acceptable quality at least in one direction. The sequences were entitled HrERAF3 to HrERAF21 (H. rhodopensis early recovery after freezing) (Supplemental Table S5) and queried to BLAST searches. Using BLASTX analysis toward the non-redundant protein database at NCBI, 19 of the sequenced cDNA fragments showed clear homologs in the database and were annotated (Table 1). In most cases, sequencing in both directions produced the same sequence. However, in the case of HrERAF6 and HrERAF7 sequencing using the M13 (-40) and ME4 primers produced completely different sequences from the same electrophoresis band, suggesting that the band contained a mixture of two amplified fragments. More than half of the BLASTX annotations were confirmed by the hits identified in the nucleotide BLAST. Nevertheless, several TDFs (HrERAF8, HrERAF9, HrERAF12, HrERAF14, HrERAF16, HrERAF17 and HrERAF20) that showed homology to hypothetical proteins in BLASX (Table 1) were also highly similar to ribosomal RNA (rRNA) genes in the BLASTN searches (Supplemental Tables S2–S4). With both oligo-dT anchored SRAP and SCoT, there were cases where bands amplified with two different primers yielded the same or very similar sequence. These included HrERAF4 and HrERAF6, amplified with SRAP primers ME3 and ME4. Other examples are HrERAF16 and HrERAF17, which were amplified with the SCoT primers SCoT4 and SCoT8.

Table 1. Homology analysis of the selected differentially expressed TDFs to gene sequences in the NCBI non-redundant protein database using the BLASTX algorithm.

Name	Amp. primer	Seq. primer	Length (bp)	Hit_id	Hit definition	Hit length (aa)	E-value	Aligned length	Identity (%)
HrERAF03	ME3	M13 (-40)	1217	ref|XP_011089949.1	Heavy metal-associated isoprenylated plant protein 3 [Sesamum indicum]	295	4.371E − 08	113	52
HrERAF04	ME3	M13 (-40)	968	ref|XP_011100803.1	Calcium homeostasis endoplasmic reticulum protein [Sesamum indicum]	672	3.8739E − 71	145	84
HrERAF05	ME4	M13 (-40)	465	gb|KAF5831363.1	Ubiquitin-related domain-containing protein [Dunaliella salina]	157	5.1218E − 77	116	99
HrERAF06	ME4	M13 (-40)	1290	gb|KZV29496.1	Calcium homeostasis endoplasmic reticulum protein-like [Dorcoceras hygrometricum]	765	1.5789E − 94	187	89
HrERAF07	ME4	ME4	1058	gb|KZV58236.1	Early nodulin-75-like [Dorcoceras hygrometricum]	312	1.4947E − 13	53	53
HrERAF08	ME5	ME5	591	gb|KAJ0804047.1	Hypothetical protein HanOQP8_Chr00c219g0792781 [Helianthus annuus]	138	8.9543E − 28	66	79
HrERAF09	ME5	ME5	699	gb|KAF7800687.1	Protein TAR1 [Senna tora]	573	2.9141E − 09	151	38
HrERAF10	ME5	ME5	1054	gb|KAF4396242.1	Hypothetical protein F8388_019788 [Cannabis sativa]	189	5.9758E − 05	87	32
HrERAF11	SCoT1	M13 (-40)	433	ref|XP_034926389.1	3-Isopropylmalate dehydratase large subunit, chloroplastic-like [Populus alba]	513	3.7931E − 13	30	90
HrERAF12	SCoT1	M13 (-40)	370	gb|PHT29099.1	Hypothetical protein CQW23_31306 [Capsicum baccatum]	324	2.351E − 41	100	88
HrERAF13	SCoT1	SCoT1	1368	gb|KZV15906.1	Putative galactinol–sucrose galactosyltransferase 2 [Dorcoceras hygrometricum]	756	0	337	83
HrERAF14	SCoT2	SCoT2	1234	gb|PHT29099.1	Hypothetical protein CQW23_31306 [Capsicum baccatum]	324	4.5979E − 24	78	90
HrERAF15	SCoT2	SCoT2	1403	gb|KZV41512.1	Hypothetical protein F511_34418 [Dorcoceras hygrometricum]	532	3.645E − 175	266	90
HrERAF16	SCoT4	M13 (-40)	1026	gb|KAI3490027.1	Hypothetical protein L1887_45856 [Cichorium endivia]	816	9.754E − 41	40	93
HrERAF17	SCoT8	M13 (-40)	1197	gb|KAI3490027.1	Hypothetical protein L1887_45856 [Cichorium endivia]	816	3.5601E − 41	40	93
HrERAF18	SCoT9	M13 (-40)	2177	gb|UYK60745.1	Photosystem II protein D1, partial [Brugmansia suaveolens] >gb|UYK60747.1| photosystem II protein D1, partial [Scopolia parviflora]	111	6.7246E − 50	98	94
HrERAF19	SCoT13	M13 (-40)	101	gb|KZV45169.1	AREB-like protein [Dorcoceras hygrometricum]	436	7.8875E − 06	31	74
HrERAF20	SCoT13	M13 (-40)	704	dbj|GEW82138.1	Hypothetical protein [Tanacetum cinerariifolium]	983	7.8983E − 28	84	73
HrERAF21	SCoT13	SCoT13	1102	ref|YP_009242020.1	NADH dehydrogenase subunit 4 [Oryza minuta] >gb|AMQ23374.1| NADH dehydrogenase subunit 4 [Oryza minuta]	2992	2.0869E − 28	176	48

Validation of expression patterns using qRT-PCR analysis

Eleven of the functionally identified TDFs were selected based on their putative functions in the RAF, to analyse their expression through qRT-PCR. Seven of the selected TDFs (HrERAF3, HrERAF5, HrERAF7, HrERAF13, HrERAF18, HrERAF19 and HrERAF21), which encode proteins with known implications in abiotic stress, showed from 2- up to 20-fold higher expression in D and ER as compared to LR (Figure 3), which suggests their involvement in the ER events. The expression of HrERAF6, encoding a calcium homeostasis protein, which is a putative housekeeping gene, was almost unchanged (data not shown), while the expression of HrERAF5, encoding a ubiquitin protein, was 2.5-fold higher in D and 1.5-higher in ER (Figure 3), suggesting it might also have functions in the recovery of H. rhodopensis. No significant changes in the expression of HrERAF20 in the three stages of RAF were observed, indicating it might indeed be an rRNA gene as suggested by the BLASTN analysis (Supplemental Table S2). Notably, HrERAF15, which encodes a protein with unknown function from another resurrection plant Dorcoceras hygrometricum, showed the second highest increase in the expression in D and ER stages (Figure 4), implying its possible involvement in the resurrection process. Therefore, three of the TDFs (HrERAF7, HrERAF15 and HrERAF19) showing homology to proteins of D. hygrometricum, were subjected to more detailed qRT-PCR analysis (Figure 4). The expression of HrERAF7, which was 3.5-fold higher at the desiccated state (0 h), slowly declined but still remained 2.5-fold higher at 24 h compared to the full recovery at 7 d. The transcript abundance of HrERAF15 declined from 7-fold to 5-fold between 0 and 5 h followed by a sharp increase at 7 h (up to 8-fold higher, compared to transcript abundance at 7 d) and decreased to fully recovered levels at 24 h. The expression of HrERAF19 slowly decreased from 3-fold to 2-fold up to 9 h followed by a faster dropdown to about half of the fully recovered levels at 24 h.

Figure 3. qRT-PCR analysis of the expression of 10 TDFs during early and late recovery stages of H. rhodopensis. D: freezing desiccated state; ER: early recovery (1 + 3 h); LR: late recovery (24 h + 7 d). Error bars on each column represent the SD of three replicates.

Figure 4. qRT-PCR expression analysis of three selected TDFs during the full course of H. rhodopensis RAF. 0 h: freezing desiccated state (8% RWC); 1, 3, 5, 7, 9, 24 h, 7 d: hours (days) after rewatering. Error bars on each column represent the SD of three replicates.

Discussion

Analysis of gene expression is an important step in many molecular and cellular biology studies as it allows the identification of differentially expressed genes that are potentially involved in stress response and tolerance. Although extensive transcriptome data, including two RNAseq experiments [26,27], are available for drought-induced desiccation and recovery of the resurrection plant H. rhodopensis, the gene expression during recovery from freezing-induced desiccation remains largely unexplored. The cDNA-AFLP (cDNA amplified fragment length polymorphism) method was previously successfully employed to identify and clone 33 TDFs that were differentially expressed in H. rhodopensis during drought-induced desiccation [30]. Despite its high specificity and efficiency in identifying differentially expressed TDFs, this method is quite demanding to apply as it requires second strand cDNA synthesis, restriction enzyme digestion, ligation and silver staining of the denaturing PAGE gels. Therefore, for this study, we were looking for fast, simple, reliable and cost-effective methods to explore and identify differentially expressed genes in the early stages of RAF in H. rhodopensis.

Luo et al. [29] proposed a new method called ‘oligo-dT anchored cDNA-SCoT’, which combines a single SCoT primer with an oligo-dT anchor primer in a PCR-based differential screening. The design of the method was aimed at the identification and cloning of full length cDNAs. Although, most of the identified TDF were not full length due to internal priming of SCoT primers used, it was demonstrated that this simple differential screening method was very efficient to identify differentially expressed genes in various stress treatments of mango [29].

In the present study, we have developed a new method named oligo-dT anchored SRAP. The method is based on the same idea as oligo-dT anchored SCoT which was modified by using a new design of the anchored oligo-dT primer in combinations with a single ‘ME’ SRAP primer. In the original SRAP method, the ‘ME’ primers were designed to target exon sequences, while the targets for ‘EM’ primers are introns [13]. The new oligo-dT primer design included the addition of M13 (-40) primer sequence at the 5′-end of the dT18 primer and the addition of the degenerate bases ‘VN’ at the 3′-end of the oligo-dT18 primer.

It has been demonstrated in well-designed experiments that most of the cDNAs synthesized during reverse transcription are products of internal priming by the oligo-dT primer, and the addition of a ‘VN’ sequence at its 3′-end leads to 3-fold reduction of the internal priming during cDNA synthesis [31]. Therefore, our novel PCR differential screening method should amplify the 3′-end of the respective genes including one or more exons, depending on the SRAP primer used. In our experiments, both oligo-dT anchored cDNA-SRAP and cDNA-SCoT with the new design of oligo-dT anchor primer were similarly efficient to identify TDFs differentially expressed in the early stages of RAF. The addition of M13 (-40) anchor sequence at the 5′-end of the oligo-dT allowed reliable direct sequencing of the identified TDFs.

In addition, in many cases sequencing data of good quality was produced also by using the respective SRAP or SCoT primer in the sequencing reaction. However, the usefulness of the SRAP and SCoT primers for direct sequencing appeared to depend on primer sequence and melting temperature. By applying both methods using 11 SRAP and 22 SCoT primers, 19 TDFs (8 with SRAP and 11 with SCoT primers) were successfully sequenced and identified. The observed high homology of HrERAF8, HrERAF9, HrERAF12, HrERAF14, HrERAF16, HrERAF17 and HrERAF20 with rRNA genes provides evidence that rRNAs are reverse transcribed by using the modified oligo-dT primer in a similar way to the generation of truncated cDNAs during reverse transcription [31] due to the presence of internal poly A stretches in the pre-rRNA sequences. Moreover, the internal reverse transcription from these ‘A-islands’ seems to be uneven and dependent on rRNA conformation in each sample leading to false selection of these fragments in the differential screening. To test this hypothesis, we examined the expression of ERAF20 by qRT-PCR, where the reverse transcription was performed using a random primers mixture, and its expression was confirmed to be nearly constitutive (Figure 3). Therefore, the efficiency of both methods for identifying truly differentially expressed TDFs could be improved by adding an rRNA depletion step before the reverse transcription.

The differential expression in the early stages of RAF was confirmed by qRT-PCR in nine out of the 19 sequence identified TDFs. Among them, HrERAF13, encoding a galactinol-sucrose galactosyltransferase 2 enzyme, showed the most elevated expression during the first hours of recovery. This enzyme (EC 2.4.1.82), also known as raffinose synthase, uses sucrose and galactinol to synthesize raffinose and inositol [32]. Raffinose is known to be synthesized and stored in seeds to protect the embryo from maturation associated desiccation and to accumulate in vegetative tissues in response to a range of abiotic stresses, including freezing. In addition, it has membrane stabilizing and antioxidant functions [33]. It has also been reported that raffinose exists in the chloroplast and may play a role in stabilizing photosystem II [34]. Our previous investigations demonstrated that inositol and raffinose accumulated during freezing-induced desiccation of H. rhodopensis [4]. Liu et al. [27] identified 22 RNAseq contigs encoding raffinose synthases that were more than 2-fold upregulated during desiccation of H. rhodopensis. Among them, CL15513.Contig6, which is most similar to HrERAF13 (Supplemental Table S3), was the most upregulated, reaching its maximum of 25-fold at 20% RWC and declining to 13-fold higher expression at desiccated state compared to the full recovery levels. It was found that galactinol-sucrose galactosyltransferase 2 transcripts were also upregulated under drought in chickpea plants [35]. The functions of raffinose synthase in enhancing drought tolerance in maize and Arabidopsis were also recently elucidated [32]. Taken together, these data strongly suggest the involvement of HrERAF13 in the RAF of H. rhodopensis. However, deciphering the mechanisms of its implication in this process requires further investigation.

The second largest increase in transcript abundance compared to fully recovered conditions was shown by HrERAF15 (Figures 3 and 4). It shows similarity to a hypothetical protein in Dorcoceras hygrometricum. Although this protein has no known function and no conserved domains have been identified, its expression pattern and the fact that it is present in another vegetative desiccation tolerant plant imply its likely involvement in the resurrection process.

Abundant transcript accumulation in early stages of RAF (4-fold) was also observed for HrERAF3 (Figure 3). It encodes a heavy metal-associated isoprenylated plant protein (HIPP). This group of proteins, consisting of metallochaperones that contain a metal binding domain and a C–terminal isoprenylation motif, have been implicated in heavy metal homeostasis and detoxification mechanisms, transcriptional responses to cold and drought and plant–pathogen interactions. A rice HIPP gene (OsHIPP41) that was highly expressed in response to cold and drought stresses was identified and its product was localized in the cytosol and the nucleus [36], suggesting a possible role of HrERAF3 in the recovery process.

One of the open reading frames of HrERAF7 showed homology with a D. hygrometricum protein annotated as early nodulin-75-like. Its expression at desiccated state (3.5-fold) slowly declined during RAF but at 24 h was still 2.5 higher than that at full recovery at 7 d (Figures 3 and 4). Early nodulin-like proteins are present also in non-nodulating plant species and have been implicated in the transport of nutrients, solutes, amino acids or hormones [37]. A recent study identified an early nodulin-like protein that was highly induced by cold stress only in a low temperature tolerant maize genotype and not in the susceptible one [38], suggesting its role in the alleviation of cold stress effects.

A big difference between transcript and protein accumulation and restoration of photosynthetic functions was observed during the recovery of H. rhodopensis after freezing-induced desiccation. The transcript levels of photosystem II (PSII) core protein D1, which is encoded by HrERAF18, remained high only during early hours of rehydration (Figure 3), while the content of D1 protein increased in the whole process of rehydration, and the electron transport through PSII was restored after 9 h of rehydration [7,9]. The observed difference in transcript and protein accumulation could be related to the protein and RNA stability and the higher demand for fresh protein synthesis in the early stages of recovery.

The transcript abundance of HrERAF19, which encodes an AREB-like protein, was 3-fold higher in the D, slowly declined afterwards but still remained 2-fold higher up to 9 h of recovery compared to full recovery at 7 d (Figures 3 and 4). AREB-like proteins are bZIP transcription factors that regulate ABA biosynthesis and could be related to the expression of LEA genes [39]. Dehydrins, a LEA2 group of proteins, were detected during freezing-induced desiccation as well as after rehydration of H. rhodopensis [4,9]. Moreover, upregulation of a TDF encoding a bZIP transcription factor during the late phase of dehydration of H. rhodopensis was previously reported [30]. The distinct expression pattern and its possible involvement in regulation of the recovery process after freezing-induced desiccation makes HrERAF19 a good candidate for further functional characterization.

The elevated expression of HrERAF5 in D and ER (Figure 3) was surprising as it encodes a ubiquitin protein. This gene is known as a housekeeping gene and is often used as a qRT-PCR reference. Ubiquitin is linked to proteins targeted to degradation. This pathway has been demonstrated to be required for both the bulk degradation of cellular proteins and the targeted proteolysis of specific regulatory proteins, including many proteins of the ABA abiotic stress response signal cascade [40]. Therefore, the observed expression pattern HrERAF5 could be caused by increased demand for protein degradation in early stages of RAF and such hypothesis merits further investigation.

The expression 3-isopropylmalate dehydratase large subunit, chloroplastic-like (HrERAF11) was slightly increased during ER (Figure 3). This enzyme is involved in the amino acid leucine biosynthesis and is located in chloroplast stroma. Wheat chloroplastic 3-isopropylmalate dehydrogenase transcripts were significantly upregulated by cold [41]. A natural mutation in the Arabidopsis gene severely impaired the growth of the plants at higher than the optimal temperatures and light intensities [42]. Therefore, this enzyme could have important functions in protecting the chloroplast and photosystems from damage during the recovery process.

HrERAF21, which encodes a mitochondrial NADH dehydrogenase subunit 4, also showed higher transcript abundance in early stages of RAF. A rice mutant of the nad4 gene was seeds and seedlings lethal [43], which implies its important role in the germination of seeds and suggests it might act in a similar way during the recovery of H. rhodopensis.

HrERAF4 and HrERAF6 amplified with SRAP primers ME3 and ME4 yielded the same or very similar sequence. Both were identified as calcium homeostasis endoplasmic reticulum proteins. Although the qRT-PCR analysis revealed that the expression of HrERAF6 was almost unchanged during the recovery process (data not shown), the encoded protein could still be involved in the RAF process. Calcium signalling is related to the relocalization of proteins and regulatory components that could strongly impact the cell functioning [44] and is an essential element for cell wall stability [45]. Chloroplasts are the primary target of ROS generated in response to stress such as dehydration and subsequent rehydration of plants. The endoplasmic reticulum is a storage compartment for calcium ions, and processes in chloroplasts are calcium-dependent [44,46].

Taken together, the presented data show that the differentially expressed TDFs identified in the present study, have high potential to function in H. rhodopensis RAF, thus confirming the efficiency of the applied differential screening approaches. However, further investigations are needed to uncover the mechanisms of their action in this process.

Conclusions

In the present study, we developed a novel, simple and reliable method for differential screening, named oligo-dT anchored cDNA-SRAP, which was used along with previously reported oligo-dT anchored cDNA-SCoT to identify differentially expressed genes during early stages of recovery after freezing-induced desiccation in H. rhodopensis. Both methods showed similar efficiency for identifying differentially expressed TDFs, and to our knowledge, are applied for the first time in H. rhodopensis. The identification of genes associated with regulation of abiotic stress underscores the effectiveness of the oligo-dT anchored cDNA-SRAP and cDNA-SCoT assays for the identification of differentially expressed genes in H. rhodopensis. The identified differentially expressed TDFs will contribute to our understanding of the defence mechanisms in the recovery process of H. rhodopensis, especially at the early stages of recovery, which is the most vulnerable period for plants.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Bulgarian National Science Fund project ‘Mechanisms of recovery from drought and frost-induced desiccation: strategies for survival of the resurrection plant Haberlea rhodopensis’ [Grant No. КП-06-H21/8].