A comparative gene co-expression networks analysis of two congener filmy ferns identifies specific desiccation tolerance mechanisms associated to their microhabitat preference.

Background Filmy-ferns (Hymenophyllaceae) are poikilohydric, homoiochlorophyllous desiccation-tolerant (DT) epiphytes. They can colonize lower and upper canopy environments of humid forest. Filmy-ferns desiccate rapidly (hours), contrasting with DT angiosperms (days/weeks). We hypothesized that desiccation tolerance in filmy-ferns would be associated mainly with constitutive features rather than induced responses during dehydration. Hence, the inter-specific differences in vertical distribution would be associated with different dynamics of gene expression within the dehydration or rehydration phases. A comparative transcriptomic analysis with an artificial neural network was done on Hymenophyllum caudiculatum (restricted to lower canopy) and Hymenophyllum dentatum (reach upper canopy) during a desiccation/rehydration cycle.Results Raw reads were assembled into 69,599 transcripts for H. dentatum and 34,726 transcripts for H. caudiculatum . Few transcripts showed significant changes in differential expression (DE). H. caudiculatum had ca. twice DE genes than H. dentatum and higher proportion of increased-and-decreased abundance of genes occurs during dehydration. In contrast, the abundance of genes in H. dentatum decreased significantly when transitioning from dehydration to rehydration. According to the artificial neural network results, H. caudiculatum enhanced osmotic responses and phenylpropanoid related pathways, whilst H. dentatum enhanced its immune system responses and protection against high light stress.Conclusions Our findings provide a deeper understanding of the mechanisms underlying the desiccation tolerance responses of two filmy ferns and the relationship between the species-specific response and the microhabitats these ferns occupy in nature. Craterostigma plantagineum glycine-rich protein CpGRP1 interacts with a cell wall-associated

Background Plant evolution has shaped several acclimation mechanisms and adaptations to withstand short and long-term periods of water deficit. However, most plants cannot survive desiccation, except for a small group of so-called resurrection plants [1,2].
Desiccation tolerance is a complex trait that involves the coordination of a cascade of molecular events divided into constitutive and inducible mechanisms to protect/repair cells and tissues against oxidative damage, and disruptions of metabolism and cell ultrastructure [2,3]. In general, the physiological and metabolic components of desiccation tolerance in resurrection plants resemble a combination of processes underlying drought stress responses and seed maturation [4,5]; however, the molecular "switches" and regulatory pathways behind desiccation tolerance are largely unknown [6]. Resurrection plants are classified into two groups according to their sensitivity, mechanisms of responses, and the velocity of water loss [7]. Those that can survive rapid water loss, such as mosses, possess constitutive (pre-existing) morpho-physiological desiccation tolerance mechanisms employing cellular repair mechanisms induced after rehydration, in contrast to those that survive only if water loss is gradual, which rely upon cellular protection mechanisms induced during dehydration.
Surprisingly, resurrection plants are also found in humid environments such as tropical and temperate rain forests [8]. One example is Lindernia brevidens, an angiosperm species from tropical rain forests [9]. The epiphytic ferns of the Hymenophyllaceae family (Pteridophyta) are another group of resurrection plants living in humid ecosystems [8]. The members of Hymenophyllaceae are called filmy ferns because they possess membranous fronds of a single layer of cells, normally lack cuticles, with undifferentiated epidermis, and lack of stomata [10,11]. Because the features of their leaves and the recurrent exposure to rapid desiccationrehydration events, these ferns evolved a poikilohydric and homoiochlorophyllous strategy, most typical of bryophytes, showing a high desiccation tolerance and the capacity to rapidly rehydrate [10,12,13].
The Hymenophyllaceae family is represented by 16 endemic species to the Chilean temperate rain forest. They are an important component of the epiphytic species diversity. Filmy ferns can colonize mostly all vertical strata of their host tree, although intraspecific habitat requirements have been reported [14,15,11]. For example, Hymenophyllum caudiculatum and Hymenophyllum pectinatum are restricted to the lower portion of the host trunk (from the ground to 1m of vertical height), where light availability is very low (10-100 μmol photons m -2 s -1 ), and humidity is high. Other species (e.g., H. dentatum, H. plicatum) extend their vertical distribution to heights above 10 m (were PDF eventually reach ≥ 1000 μmol photons m -2 s -1 ) [11,15,13]. Both light intensity (PFD) and vapor pressure deficit (VPD) increase with the height of the host plant, whereas relative humidity decreases.
Therefore, a species that reaches the top of trunks is prone to suffer frequent desiccation and photoinhibition. The wide spectrum of habitats, from very sheltered, high steady humidity environments to higher light, lower humidity along the vertical distribution over host trees offers a unique opportunity to address questions about how molecular and physiological mechanisms are shaped in congeneric populations differing in microhabitat preferences and ability to tolerate desiccation. Questions such as, are the molecular mechanisms responsible for protection/repair of these filmy ferns against tissue desiccation inducible during dehydration or constitutive (pre-existing), such as those exhibited by rapidly desiccating bryophytes? Is the dynamic of gene expression similar between filmy fern species from lower light, higher humidity environments with species from higher light, lower humidity environments? Do these filmy ferns invoke similar gene functions during a cycle of desiccation-rehydration? As filmy ferns are frequently and rapidly dehydrated, it has been proposed that desiccation tolerance in these ferns would be associated mainly with constitutive features rather than induced responses during dehydration.
However, we hypothesize that the inter-specific differences in vertical distribution would be associated with different dynamics of gene expression within the dehydration or rehydration phases. Published data partially support this idea [16,17]. Advances in Next Generation Sequencing tools have brought important advantages to uncover underlying mechanisms that control plant responses associated with natural constraints explored in situ or under experimental conditions. Here we used RNA-seq on the Illumina Hi-seq platform to study and compare the transcriptional responses of H. caudiculatum and H. dentatum, two Hymenophyllaceae species with contrasting vertical microhabitat preferences along host trees and different rates of water loss (Fig. 1). Specifically, we examined at the dynamics of gene expression and a Weighted Gene Co-expression Network (WGCNA) [18] coupled with neural artificial networks [19 -21] of fronds subjected to experimental desiccation-rehydration cycles to identify commonalities and differences on gene expression dynamics associated to its water status underlying their resurrection strategy.

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 lowquality 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.  (Fig. 3b), which is consistent with their poikilohydry strategy and the regressive evolution hypothesis [22].

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 [139] 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).
In parallel to the differential expression analysis, we also observed all those highly abundant transcripts based on the FPKM value but with a fold change vale < 2, i.e., constitutive highly abundant genes without significant change along the dehydration-rehydration cycle. A total of 156 transcripts for H. caudiculatum, and 199 for H. dentatum show that there are broadly three main processes involved in the desiccation tolerance of these two filmy ferns, namely translation, photosynthesis and antioxidant activity (See Additional file 2: Dataset S3 for the detailed list of the highly abundant constitutive transcripts in each species).

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 SOMbased 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 ( Regarding H. dentatum, there were no hub genes in none of the resultant genecoexpression 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

Discussion
Previous studies have indicated at different scales, from cellular to ecophysiological approaches, how Hymenophyllaceae species from temperate rain forest responds to desiccation [12,13,15,16,17,23]. However, the link among desiccation response, habitat preferences and gene expression have never been addressed before in Hymenophyllaceae until our present study. Both de novo assembled transcriptomes (the first for the Hymenophyllaceae family) exhibited good quality parameters, comparable to other non-model species transcriptomes sequenced with the same Illumina platform and using the same criteria for subsequent downstream assembly pipelines [24 -26]. Thus, the robustness of our transcriptomes allowed us to explore differential gene expression and to contrast the dynamical patterns of highand-low abundance of genes during the transition from dehydration to rehydration versus the constitutive gene expression, thereby contributing to a comprehensive view of the response mechanisms involved in a fast dehydration-rehydration process. Finding patterns in massive gene expression data sets consistent with a biological function is always a challenge. Despite the undeniable advantages of exploratory statistical models such hierarchical clustering [27] and k-means clustering [28], these methods are subjective due to human bias based on arbitrary statistical significance threshold and does not consider the topology between clusters [29]. The characteristics of the artificial neural networks, considering topology in clusters neighboring [19,20], provides an excellent tool to depict clustering patterns of gene expression across multiple factors [19,30,31]. With the differential gene expression analysis, we found similar patterns of molecular responses between both filmy ferns during the desiccation-rehydration cycle. Then, by combining the Self-Organizing Maps with weighted gene co-expression network analysis, we showed clear patterns of transcript accumulation and found a clear relationship between core genes that are co-expressed in a given hydration state.
Thus, we were able to elucidate both shared and species-specific molecular components such as mechanisms of signaling and responses to desiccation [6,32], disturbance of chloroplast homeostasis [33,34], and photoinhibition-related stresses [35 -37], associated with frond hydration state and the spectrum of microhabitat that H. caudiculatum and H. dentatum occupy along the vertical distribution over host trees.
At a fully hydrated state, the patterns of transcripts accumulation in both species reflect normal functioning of primary metabolic process such as photosynthesis and respiration. After dehydration, the clustering patterns and abundance of differentially expressed identified transcripts, shared mechanisms associated with the maintenance of redox homeostasis, stabilization and maintenance of the photosynthetic apparatus, and chloroplast operational signaling. In contrast, the species-specific responses were particularly associated to phenolic compounds biosynthetic pathways (e.g., phenylpropanoid metabolism) and photoinhibition related processes. Both general and species-specific responses agree with physiological, ultrastructural, and chemical changes, as with the homoiochlorophyllous strategy of these filmy ferns [11 ̶ 13, 15, 23].

Rehydration Cycle.
It must be noted that because the relatively few transcripts showing differential expression, the desiccation tolerance response of these two filmy ferns would rely mainly on constitutive mechanisms, as it has been suggested by proteomic analyses [16]. Based on our results, the both fern species have high abundance of transcripts are key components for desiccation tolerance [3,6,36]. , Given that GST enzymes inactivates a wide variety of toxic compounds (e.g., hydrophobic, electrophilic and cytotoxic substrates; see [38]), these two ferns appear to have evolved a mechanism to finely tune the GST responses against oxidative damage during the desiccated state. We found the same dynamic in genes encoding for ferritins which could be a more specific mechanism coupled with the aforementioned for GST. Because the fronds of these ferns ranges from one to a few cells layers, and given that plastid-to-nucleus signaling is essential for the coordination and adjustment of cellular responses to external and internal cues [33,46,47], environmental perturbations imposing restrictions for the chloroplast homeostasis would have a great impact in whole frond form and function. Fluctuations in frond water status will initiate an operational signaling to change the expression of thousands of genes to adjust chloroplast and the whole cell functioning [34]. According to this, both species showed an increase of transcript abundance that encode for AP2/ERF transcription factors during desiccation (e.g., DREB, ERF, RAP), which are known to be involved in responses to environmental stimuli, especially for redox, osmotic and drought stress [33,34]. Moreover, a coordination of chloroplastic and nuclear gene expression would be mediated by the redox status [48], activating a protein kinase phosphorylation cascades and AP2/ERF transcription factors. Our results are consistent with the above-mentioned mechanism, being a key mechanism contributing to the resurrection strategy used by the two studied filmy ferns.

Species-Specific Patterns of Transcriptional Responses During Dehydration and Rehydration.
One of the most ambitious goals of our study was to associate the different niche preference of these two filmy fern species with species-specific transcriptional responses to have insights into particular physiological conditions that may prevail under different frond hydration states. In general, based on the DE analysis and the pairwise comparison among hydration states (Table 1) caudiculatum during the dehydration-rehydration cycle, we observed an increase in the abundance of the DOG1 gene when full hydrated and rehydrated. This gene is involved in controlling the timing of seed dormancy in response to environmental signals such as temperature by altering gibberellin metabolism [49]. Given that during embryo development and maturation the acquisition of desiccation tolerance is a key process for the viability of the majority of seeds, it is logical that resurrection plants use part of the molecular components involved in the control of seed dormancy. However, the changes in the abundance of DOG1 gene in H.
caudiculatum suggest an induced mechanism linked to ABA and GA balance by a common intermediary such as DELLAs proteins [50]. Based on the genecoexpression network results, we observed for the dehydration state that core genes are related to minimize ROS formation and oxidative stress through NAD metabolism and mitochondrial uncoupling proteins [51,52], while for rehydration all co-expressed genes are involved mainly in cell wall and protein stability [53], membrane stability and signaling, and the reestablishment of glucose metabolism (Fig. 7a).
We observed an increase in the abundance of the ELI6 gene which encodes an early light inducible protein (Table 1) in H. dentatum when full hydrated and rehydrated.
The same patter was observed for the dehydration response element transcription factor DRE2A. ELIPs proteins are induced by light, and act as sinks for excitation energy under high light. The increase in abundance ELI6 and DRE2A only in H. dentatum (Fig. 7b) would be a response associated with its distribution toward more illuminated and low humidity micro-sites along the vertical direction of host trees.
Other specific responses are linked with the activation of defense responses driven by Ethylene and Salicylic Acid, and the scavenging of specific ROS such as hydrogen peroxide.

Conclusions
The present work provides a deeper understanding of the mechanisms underlying the desiccation tolerance responses of H. caudiculatum and H. dentatum. Here we have deciphered the transcriptional dynamics, gene co-expression networks and key genes involved in the signaling and response mechanisms during the process of dehydration and rehydration for both species. In addition, the identification of species-specific mechanisms contributes to explain their ecophysiological traits and microhabitat preferences they show in their natural environment. Thus, while the lower-canopy species H. caudiculatum seems to enhance osmotic responses and phenylpropanoid related pathways, the upper-canopy species H. dentatum enhanced its defense system responses and protection against high light stress.
Finally, much more work is needed to understand the roles and contribution of specific molecules (e.g., RNA-binding proteins) and posttranscriptional regulations that may have central roles in the desiccation tolerance strategy of these species when subjected to rapid loss-and-gain of water as occurs in the natural environment. Also, studies under an evolutionary developmental biology approach would help to decipher and understand the genetic basis of physiological, developmental and morphological variation associated with the microclimatic conditions in which these resurrection filmy ferns have evolved.

Individuals of Hymenophyllum caudiculatum (Mart.) var productum (K. Presl.) and
Hymenophyllum dentatum (Cav.) (Fig. 1) were collected from Katalapi Park (41°31'07.5" S, 72°45'2.2" W) Cordillera de Quillaipe, Los Lagos Region, Chile. The collecting site correspond to a coastal evergreen temperate rainforest. Details on microclimatic conditions of the site have been previously published [11]. Small pieces of bark or fallen trees covered with epiphytic filmy ferns were collected in late Spring (Nov-Mid Dec, Southern hemisphere) and transported to a shaded experimental nursery garden with automated irrigation at the Universidad de La Frontera, further description of nursery conditions and plant set up can be found in [17]. Briefly, plants were acclimated to nursery garden conditions under shade (25-30 μmol photons m −2 s −1) nebulized by an automated irrigation system (6 daily waterings, 2 min each) for 25 days prior to experimental procedures.
It is worth to note that the deposit of specimens to an herbarium was not necessary, since they already exist in the collection of the Herbarium of the Universidad de Concepción (http://www2.udec.cl/~herbconc/index.htm, Herbario CONC, Departamento de Botánica, Universidad de Concepción, Barrio Universitaro s/n, Casilla 160-C, Concepción, Chile.)

Experimental Design and Sample Collection
The experimental design was a desiccation-rehydration experiment under nursery conditions during late spring (Nov-mid Dec, southern hemisphere) in which fronds of H. dentatum and H. caudiculatum were studied at three different hydration states; full hydrated (FH), dehydrated (DH), and rehydrated (RH) (detailed in Fig. 2a). At each hydration state, we used attached fronds for monitoring the changes in maximal quantum efficiency (Fv/Fm), and two set of detached fronds; one to determine their relative water content (RWC), and the other for RNA isolation. To reach a fully hydrated state, ferns were subjected to irrigation pulses of three minutes at intervals of twenty minutes each, during four hours; a single frond from each plant (3 plants per species, 6 plants total) were taken at this hydration state (FH) to obtain its relative water content (RWC). Then, ferns went through a desiccation-rehydration process by adjusting the irrigation settings and monitoring the changes in RWC and Fv/Fm of the fronds during the next seven days. AT the first day without irrigation, the first DH sampling (fronds from three different individuals per species) was taken when fronds reached about 60% of RWC. A second DH sample was obtained when plants reached a value equal or below the critical relative water content reported for these species [11], achieved at day seven (Fig.   2a). Both DH samples were pooled in order to have coverage for early and late desiccation-regulated transcripts. Finally, the rehydration was carried out resuming irrigation until at least a 90 % of RWC was recovered. During rehydration (RH) two samplings were taken, one was obtained when plants recovered 60-70% of RWC, and the second sampling was taken at ca. 90% of RWC. Both samplings were pooled as RH in order to have coverage for early and late rehydration-regulated transcripts.
Samples were stored at -80° C until processed for RNA isolation for comparisons of the transcriptional responses associated to each hydration state, either within or between species.

Fluorescence Measurements
Maximal quantum efficiency (Fv/Fm) was calculated from chlorophyll a fluorescence signals obtained from attached fronds using a modulated fluorometer (FMS2, Hansatech, UK) as previously described in [17]. Briefly, attached fronds were carefully cover and dark acclimated for 30 min using FMS2 leaf clips. After dark acclimation the modulated light was turned on to obtain F0, then a saturating pulse (800 mS at ~3,000 µmol m -2 s -1 ) was applied to obtain the maximum fluorescence (Fm). Variable fluorescence (Fv) and Fv/Fm ratio were calculated according to [54] at ambient temperature (17 C° approx.). Data was analyzed after check normality assumptions under an ANOVA test with P value ≤ 0.05. When data did not meet the normality assumptions, we used the non-parametric test of Kruskal-Wallis.
Additionally, fluorescence images of Fv/Fm were obtained using detached fronds under each hydration condition FH, DH and RH using a Maxi-Imaging PAM (Walz, Effeltrich, Germany) to observe the rate of recovery of the quantum yield of fluorescence (YII) at the whole frond level of both species (Fig. 2b).
Reads from each sample were mapped to their corresponding final transcriptome using default RSEM parameters (RSEM = RNA-Seq by Expectation Maximization) written in the align_and_estimate_abundance.pl script. The resulting RSEMestimated gene abundances for each fern were merged in a matrix and analyzed with run_DE_analysis.pl script from Trinity, which involves the Bioconductor package edgeR in R statistical environment [58,59]. Transcripts with very low estimated counts (2 for combined groups), were not considered for edgeR pair-wise comparison of hydration states. To judge significance of gene expression, we used a False Discovery Rate value (FDR) lower than 0.05 and a minimum fold change (FC) of 2 as thresholds (Additional file 2: Dataset S1 and S2). To quantify transcript abundance, normalized RSEM-estimated counts were used for clustering assembled contigs based on expression patterns [29]. Finally, the resulting transcripts of each transcriptome were aligned into the SwissProt database using BLAST+ with an evalue filter of 1-e -10 as threshold [60]. Functional annotation, classification, and over-or under-represented groups of genes were performed using PANTHER ( [61], www.pantherdb.org).

Self-Organizing maps (SOM) Analysis
Normalized RSEM-estimated counts of both H. caudiculatum and H dentatum that met the expression values determined from the model described above were used for the SOM clustering method [29]. Specifically, only genes that vary significantly in expression across fronds hydration state of both filmy fern species were analyzed. To focus only on gene expression profile, and at the same time avoid biases with the differences in the magnitude of gene expression, expression values were mean centered, selected from the upper 50% quartile of coefficient of variation, and variance scaled using the scale function (R base package; [59]) separately for H. caudiculatum and H. dentatum. Scaled expression values were used to cluster genes in both species across fronds hydration states into a multidimensional 2 x 3 hexagonal SOM using the Kohonen package on R [19]. One hundred training iterations were used during clustering with a decrease in the alpha learning rate from ca. 0.0018 to 0.0010 (Additional file 3: Fig. S2). After the iteration process, the final assignment of genes to the winning units shaped the clusters of genes (termed nodes) associated to the hydration states. SOM outcome was visualized into pie charts for codebook vectors to obtain the counts number and mean distance of the genes assigned to each node ( [19], Additional file 3: Fig. S3).
The box plot option from the ggplot2 package on R was used to visualize the gene accumulation patterns associated to the hydration states of the fronds in each Node. Finally, the genes of each Node were analyzed for GO enrichments terms at a 0.05 false discovery rate cutoff.

Gene Co-expression Network Analysis
In order to further explore the gene cluster generated by the SOM analysis, a Gene Regulatory Network-based (GRN) approach was used to study the interactions between gene expression and light availability in the habitat. From the SOM clustering method, a subset of 67 and 183 annotated genes from H. dentatum and H. caudiculatum were selected, respectively. For the selection of genes, those transcripts from SOM nodes with the highest scaled accumulation patterns for a given hydration state were chosen. Then, only those genes with GO annotation were selected. (Additional file 2: Dataset S4 and S5). These genes were used to construct a weighted gene coexpression network by using the Weighted Gene Coexpression Analysis (WGCNA; [18]) and the igraph [21] R packages. In order to meet the scalefree topology criteria for WGCNA, a soft threshold (β) value of 9 was used for the calculation of the adjacency matrix. The resources and algorithms of the graph generator and community structure functions of the igraph R package were used to explore the network properties such as connectivity, centralization, modularity and community structure. Finally, custom graph functions of the igraph package were used for network visualization (available at https://github.com/eostria/Gavel).

Sequence Submission
The quality-filtered, barcode-sorted, and trimmed short read data set used for transcriptome assembly and gene expression analysis, was deposited in the NCBI

Ethics approval and consent to participate
Plant collection was made on the Katalapi Park Foundation, a non-profit institution for research and environmental education. Dr. Luis J. Corcuera, President of the Foundation provided the permission to access and collect the plant material from the Park. There were no need for special permissions from any other government or private agent for collecting the plant material used in our study. The fern species were identified by EO and LB based on the description of the species published on [62]. The experimental design was made following the institutional normative of the Ethics Committee of the University of La Frontera.

Availability of data and materials
The transcriptomes of the ferns species used in this study are available in the NCBI The custom script used for the analysis and construction of the gene co-expression networks is available at https://github.com/eostria/Gavel All data generated and analyzed in this study are available and included in this manuscript, as well as its supplementary information files.

Competing interest
The authors declare that they have no competing interest.      GO-category distribution of annotated genes for H. caudiculatum and H. dentatum among lev Figure 6 Mapping quality, clustering results for self-organizing maps (SOM) and GO enrichment of tran