Changes in Ferns’ Relative Water Content and Maximum Quantum Efficiency During Desiccation-Rehydration Cycle
The two filmy fern species showed a similar rate of dehydration during the first three hours after cessation of irrigation (beginning of dehydration), both reaching ca. 60% of relative water content (RWC) (Fig. 2a). From 3 to 25 hours without irrigation, H. dentatum losses water faster than H. caudiculatum,reaching 18% and 30% RWC respectively. During this period of dehydration, the maximum quantum efficiency (Fv/Fm) drastically decayed from about 0.7 to 0.2 in H. dentatum but remained nearly 0.78 in H. caudiculatum (Fig. 2a insert). After a week without irrigation, both reached a RWC between 11 to 17% and a Fv/Fm near to 0.2 (Fig. 2a). When irrigation was reestablished, H. caudiculatum had a faster rehydration and recovery of Fv/Fm compared to H. dentatum. Nevertheless, both species reached similar maximum photochemical efficiency (Fv/Fm) by the end of the experiment. Recovery of photochemical efficiency at the whole frond level was confirmed by imaging fluorescence (Fig. 2b).
Transcriptional Profile and Transcripts Annotation
A total of 111,495,169 and 110,988,488 paired-end reads (101 bp) were obtained after sequencing libraries of H. caudiculatum and H. dentatum,on the Illumina HiSeq2000 platform (Additional file 1: Table S1). Following the removal of low-quality reads and duplicated reads, we performed a de novo transcriptome assemblies with Trinity software by using a set of ~85 million reads for H. caudiculatum and ~87 million reads for H. dentatum (Additional file 1: Table S1). The initial assemblies resulted in 161,689 contigs for H. caudiculatum and 332,003 contigs for H. dentatum, which were refined to remove low supported transcripts. Transcripts with an estimated abundance lower than 1 FPKM and highly similar or redundant transcripts with a sequence similarity higher than 95% were removed. The resulting transcriptomes are represented by 34726 contigs for H. caudiculatum and 69599 contigs for H. dentatum (Additional file 1: Table S2). Although the number of transcripts decreased significantly during the refinement, ca. 80% and 70% of high-quality reads were mapped to the H. caudiculatum and H. dentatum transcriptomes, respectively.
The final transcriptome assemblies were aligned to the SwissProt database for annotation, with an alignment rate of ca. 50% for the transcripts of each transcriptome. In spite of the low identification rate, most of the unknown transcripts (~80% in H. dentatum, and~65% H. caudiculatum; Additional file 1: Table S3) belonged to small size transcripts (< 1000 bp; see Fig. 3a). An insight into the taxonomic distribution of top blast hits of transcripts revealed that both H. caudiculatum and H. dentatum had among their top hits a high amount of sequences belonging to the model moss Physcomitrella patens and the lycophyte Sellaginella moellendorffii (Fig. 3b), which is consistent with their poikilohydry strategy and the regressive evolution hypothesis .
Differential Expression Analysis and Functional Annotation
After quality filtering and refinement, we examined the expression dynamics of annotated genes during the desiccation-rehydration cycle by pairwise comparisons by using a fold change ≥ 2 and a FDR < 0.05 as cut-off (Additional file 2: Dataset S1 and S2). In both species, few transcripts showed significant changes in differential expression (DE). For H. caudiculatum, the highest number of DE genes occurred during the dehydration process with a total of 265 DE genes, where most of them  showed an increase in their abundance (Fig. 4). In H. dentatum, the number of DE genes that increase and decrease their abundance were similar. However, among the different hydration states, the abundance of genes decreases significantly when transitioning from dehydration to rehydration. When comparing both species, H. caudiculatum presents ca. twice DE genes of H. dentatum and a higher proportion of both, increase and decrease abundance of genes under the dehydration process (Fig. 4).
From the differentially expressed transcripts we explored the function of the gene products by conducting a Gene Ontology analysis (GO) (see Material and Method for details). Both species showed similar enrichment pattern of sequences for each GO category (Fig. 5). For example, at the biological process category (BP), there was a high number of sequences in the metabolic process (> 4000 sequences) and in cellular process (> 2000). In the cellular component category (CC), organelle and cell part components showed the highest enrichment of sequences (e.g., ca. 1000 and 1500 sequences for organelle, for H. dentatum and H. caudiculatum, respectively). Finally, in the molecular function category (MF), the highest accumulation of sequences was found to be in the antioxidant activity, followed by binding process (Fig. 5).
Transcripts Clustering Patterns and Gene Co-Expression Network of H. caudiculatum and H. dentatum Across Dehydration-Rehydration Cycle
From the annotated DE genes, we studied the dynamics of gene expression in these two resurrection filmy ferns to identify, firstly, transcripts with similar accumulation patterns in response to a given hydration status, and secondly, the complexity of the interaction and co-expression of genes related with their desiccation and rehydration responses. The self-organizing maps (SOM) partitioned the DE genes into six clusters (hereafter nodes) arranged as a map (Fig. 6 a-c). The underlying topology of the SOM shows distinct accumulation patterns of genes membership across nodes (Fig. 6 a-c), and prominent densities of transcripts for a given hydration state within a node (Fig. 6 b-d), which in topological terms, reflects similar accumulation patterns. Thus, SOM Node 6 in H. caudiculatum and SOM Node 1 in H. dentatum have a prominent density of transcripts associated with the full hydrated state (FH). For the dehydrated and re-hydrated conditions, respectively, the enrichment of transcripts was observed in Nodes 3 and 2 in H. caudiculatum (Fig. 6b), and Nodes 2 and 6 in H. dentatum (Fig. 6d). In order to ascertain if SOM-based clustering yields biologically relevant information, we determined GO enrichment in those nodes. Specifically, we found that for H. caudiculatum, transcripts of the dehydrated state (Node 3) were enriched in functional categories related to stress signaling and response, photosynthesis and photosystem II stabilization and repair, unsaturation of fatty acid, and lignin biosynthetic process. Next, transcripts in the rehydration state (Node 2) were enriched in responses to oxidative stress, lignin biosynthetic process, photosynthesis, protein-chromophore linkage, cellular redox homeostasis, and translation. On the other hand, in H. dentatum, the dehydrated state (Node 2) was enriched in functional categories related mainly with antioxidant responses such as glutathione metabolic process and ROS detoxification systems, drought response, transcription, translation regulators, photosystems II stabilization, ATP synthesis and proton transport, photoprotection, and ABA non-regulated stress responses. For the rehydrated state (Node 6), transcripts were enriched in functional categories corresponding mainly to ethylene and abscisic acid related signaling, photosynthesis, proton transport, plasmodesmata-mediated intercellular transport, response to stress and to toxic substances.
From these SOM nodes (2,3 and 6 in H.caudiculatum, and 1, 2 and 6 in H. dentatum),, we carefully reviewed the annotated genes with high scale expression to construct weighted gene coexpression networks for each node (Fig. 7; see Materials and Methods for details). Based on the Fast Greedy modularity optimization algorithm for finding community structure, all gene-coexpression networks had two modules. The networks obtained from the full hydrated state of H. caudiculatum (Node 6, Fig. 7a), and from the dehydrated state of H. dentatum (Node 2, Fig. 7b) showed the highest gene connectivity, i.e., interactions between genes (6,216 and 866 connections, respectively). An overview of the resultant gene-coexpression network for each hydration state of H. caudiculatum (Fig. 7a) showed that at full hydration, twelve hub genes (> 200 connections) composed a core network connecting both modules, which were involved in: protective system against oxidative stress (GPX7, CAT3);; in light harvesting complexes and reaction centers of photosystems I and II (e.g., CB23, CP24, LHCA 4, LHCB 2, PSAA, PSBA);; lipid metabolism and transport (ACLA3, NLTP5).. Under desiccation, the resultant network showed six hub genes (> 100 connections), and they were part of one of the modules. Hub genes were involved in: cell wall reinforcement (WUN1);; glutathione metabolism (e.g., GSTUJ, GSTX4);; mitochondrial uncoupling protein (PUMP4);; glycosylation (U85A3);; nitrogen mobilization (NRTs). At the rehydration state, the resultant network contained only eighteen genes, all with the same number of connections among them. Two main process were represented by these genes, namely cell wall structure and architecture (e.g., PME53, PRP1, CSE),, and stress response and signaling (e.g., GRP, ASR1, DSP22, TET8)..
Regarding H. dentatum, there were no hub genes in none of the resultant gene-coexpression networks (Fig. 7b). At the full hydration state, the network contained twelve coexpressed genes, equally connected among them. They were grouped into different functional categories, such as plant defense and abiotic stress resistance (e.g., DIR5, ERF17, BURP16),, chloroplast development (e.g., PBP1),, and cell wall structure and architecture (e.g., PME53).. Under desiccation, the resultant network contained a total of forty-two coexpressed genes. Among them, we found genes involved in flavonoid biosynthesis (CHSY),, structure of photosystems I and II (PSAA, PSBS),, stress response (GL82),, ubiquitination (UBIQP),, aldo-keto reductases (ALKR4),, immune system and Salicylic acid homeostasis (DLO2),, non-symbiotic hemoglobins (HBL)which function in signal transduction pathways of several hormones, Jasmonic acid precursor (OPR7).. Lastly, at the rehydration state the resultant network contained only 6 coexpressed genes. They were mainly involved in membrane metabolism (e.g., GDPD),, stress response and signaling (e.g., DIR5, GRG1, DRE2A),, and plant immunity (e.g., PUB25, TGA2.1),, and sugar metabolism (GRG1)..