Analysis of protein-protein interaction and weighted co-expression networks revealed key modules and genes in multiple organs of Agave sisalana

Carvalho, Lucas M.; Silva, Nicholas Vinícius; de Abreu, Luís Guilherme F.; Marone, Marina Püpke; Cardelli, Alexandra Russolo; Raya, Fabio Trigo; Araújo, Guido; Carazzolle, Marcelo Falsarella; Guimarães Pereira, Gonçalo Amarante

doi:10.3389/fceng.2023.1175235

ORIGINAL RESEARCH article

Front. Chem. Eng., 26 May 2023

Sec. Sustainable Process Engineering

Volume 5 - 2023 | https://doi.org/10.3389/fceng.2023.1175235

This article is part of the Research Topic The Role of Agave as Feedstock within a Sustainable Circular Bioeconomy View all 6 articles

Analysis of protein-protein interaction and weighted co-expression networks revealed key modules and genes in multiple organs of Agave sisalana

Lucas M. Carvalho^1,2

Nicholas Vinícius Silva¹

Luís Guilherme F. de Abreu¹

Marina Püpke Marone¹

Alexandra Russolo Cardelli¹

Fabio Trigo Raya¹

Guido Araújo^2,3

Marcelo Falsarella Carazzolle^1,2

Gonçalo Amarante Guimarães Pereira¹*

¹Laboratory of Genomics and BioEnergy (LGE), Department of Genetics, Evolution, Microbiology, and Immunology, Institute of Biology, State University of Campinas, Campinas, SP, Brazil
²Center for Computing in Engineering and Sciences, State University of Campinas, Campinas, SP, Brazil
³Institute of Computing, State University of Campinas, Campinas, SP, Brazil

Agave plants are well-known for their drought resilience and commercial applications. Among them, Agave sisalana (sisal) is the species most used to produce hard fibers, and it is of great importance for semiarid regions. Agaves also show potential as bioenergy feedstocks, as they can accumulate large amounts of biomass and fermentable sugar. This study aimed to reconstruct the A. sisalana interactome, and identify key genes and modules involved in multiple plant tissues (root, stem, and leaf) through RNA-Seq analysis. We integrated A. sisalana transcriptome sequences and gene expression generated from stem, leaf, and root tissues to build global and conditional co-expression networks across the entire transcriptome. By combining the co-expression network, module classification, and function enrichment tools, we identified 20 functional modules related to at least one A. sisalana tissue, covering functions such as photosynthesis, leaf formation, auxin-activated signaling pathway, floral organ abscission, response to farnesol, brassinosteroid mediated signaling pathway, and light-harvesting. The final interactome of A. sisalana contains 2,582 nodes and 15,083 edges. In the reconstructed interactome, we identified submodules related to plant processes to validate the reconstruction. In addition, we identified 6 hub genes that were searched for in the co-expression modules. The intersection of hub genes identified by both the protein-protein interaction networks (PPI networks) and co-expression analyses using gene significance and module membership revealed six potential candidate genes for key genes. In conclusion, we identified six potential key genes for specific studies in Agave transcriptome atlas studies, biological processes related to plant survival in unfavorable environments and provide strategies for breeding programs.

1 Introduction

Plants of the genus Agave are known for their drought resistance mechanisms and commercial applications. Agaves are native to semiarid and arid regions of North and Central America and possess several morphological, anatomical, and physiological drought resistance mechanisms (Eguiarte et al., 2021). One of the most important is CAM (crassulacean acid metabolism), the most water-use efficient type of photosynthesis (Borland et al., 2014; Yin et al., 2018). Also, several Agave species have established commercial uses, for instance, the production of tequila (Agave tequilana), mezcal (Agave angustifolia and others), and sisal fibers (mainly Agave sisalana). Sisal is the main source of natural hard fibers, and, due to its drought resistance, it has a great social-economical impact on semiarid regions (Broeren et al., 2017; Monja-Mio et al., 2019). In addition, A. sisalana has interesting characteristics to be used as a bioenergy feedstock in semi-arid regions (Davis et al., 2011), since only 4% of A. sisalana leaves are used to produce commercial fibers (Suinaga et al., 2007). The remaining waste could be used to produce biogas or second-generation ethanol (Mielenz et al., 2015). Brazil is the greatest sisal fiber producer in the world and has the capacity of increasing the area in which this crop is produced to keep up with the increasing demand for biofuels (International Energy Agency, 2021).

Transcriptomic atlases (Sarwar et al., 2019; Raya et al., 2021) and phylogeny and evolution studies of A. sisalana (Jiménez-Barron et al., 2020) are available and provide valuable information about its genetics. However, there are no analyses or available data on protein-protein interactions (PPI) in any Agave species. PPI networks represent physical interactions between proteins, which can lead to useful insights into cell physiology and help identify targets for genetic engineering (Ferreira et al., 2016; Van den Broeck et al., 2017; Khan et al., 2020; Yu et al., 2020; de Silva et al., 2022). PPI associations, for example, may represent a functional cell action and have already been used in several areas of the biological sciences (Braun et al., 2013; Yeger-Lotem and Sharan, 2015; Ding and Kihara, 2019; Yang et al., 2019). High-throughput technologies, such as mass spectrometers and genomic sequencers, generated a large amount of data allowing for an increase in the number of known protein interactions of various organisms, especially model organisms, and are available in several databases, such as STRING (Szklarczyk et al., 2011).

However, the interactions between proteins available to non-model organisms are very limited, so the complete mapping depends on bioinformatics strategies that aim to assign relationships through orthologs (Uhrig, 2006). These strategies for assembling PPI networks assume that proteins are conserved across different species, and therefore, the interactions between two proteins may also exist in those species (Wang and Jin, 2017). In addition, multi-omics analyses integrated with co-expression analysis in the understanding of complex phenotypes of non-model plants bring relevant information about the biological system, genes pattern, unexplored genes, and their relation with the environment (Abraham et al., 2016; Sarwar et al., 2019; Yu et al., 2020; Feng et al., 2021; Panahi and Hejazi, 2021; Li et al., 2022).

Here, we predict A. sisalana interactome using public transcriptome sequences generated from leaf, stem, and root tissues. Moreover, we developed a new strategy that integrates network analysis with weighted gene coexpression network analysis to identify correlations between genes across tissues based on coexpression relationships and focus on pointing out key genes and modules, which will aid current research and provide a framework for future Agave interactomics research.

2 Material and methods

2.1 RNA-Seq samples

In this study, we investigated the sisal using transcriptomic datasets previously generated by Raya et al. (2021). Our samples were extracted from sections of the leaves, stems, and roots from three different seven-year-old adult plants of A. sisalana collected at the EMBRAPA collection in Monteiro, in the state of Paraíba, Brazil (07°53′ S; 37°07′ W, elevation: 619 m). This region faced an irregular rainfall regime before sampling. For each cultivar, three biological replicates growing side-by-side were sampled. To maintain leaf maturity equivalency within the cultivars, we sampled the fifth leaf counted from the central spike of each plant. Although the plants were maintained in the field and exposed to long drought periods typical of the Caatinga biome, all the collected individuals were healthy, with the leaves showing homogenous green coloration, with no visible symptoms of diseases (necrosis, chlorosis, spots, etc.).

2.2 Transcriptomic analysis

The assembled transcripts and the raw RNA-Seq reads of three tissues (biological triplicates of leaf, stem, and root) of A. sisalana were downloaded from NCBI SRA (accession number PRJNA746623) (Raya et al., 2021). Poor quality reads and adapters were removed with the Trimmomatic v.0.39 (Bolger et al., 2014) software. Quantification of the assembled transcripts was performed by the software kallisto v.1.0.4 (Bray et al., 2016) and normalized by TMM using the edgeR package v3.38.2 (Robinson et al., 2010). The ORFs of the assembled transcripts were identified by the Transdecoder software with the parameter -m 300 (ORFs longer than 300 nucleotides). Transcripts that did not have ORFs and whose expressions were less than 1 TPM in all conditions were removed. After this removal step, a new round of quantification was performed. A. sisalana proteins were obtained from the translation of ORFs.

2.3 Grouping of orthologous families

Five organisms was selected to construct the orthogroups, together with A. sisalana proteins, based on their wide availability of PPI data in the STRING v.11 database: Arabidopsis thaliana, Musa acuminata, Oryza sativa, Sorghum bicolor, and Zea mays. Protein sequences from each organism were retrieved from the NCBI website. OrthoFinder v.2.3.3 software was used (Emms and Kelly, 2019) to construct orthogroups (gene families) from these organisms. From the gene families obtained, we selected those that are considered one-to-one or one-to-many. One-to-one gene families are those that have one copy of a gene in A. sisalana and one copy in every other organism. The one-to-many gene families mean that the gene of interest has one copy in A. sisalana and more than one ortholog in the other species (Zahn-Zabal et al., 2020). The reconstruction of the phylogenetic tree between the species was also generated by OrthoFinder from the default parameters.

2.4 Construction of A. sisalana interactome

For each gene family containing at least oneA. sisalana protein and proteins from other organisms, we identified the PPIs present in each orthologous organism in the STRING database, considering a cutoff score of 0.70 for each protein interaction. For example, if an A. sisalana gene A has as orthologs the gene B in S. bicolor and the gene C in M. acuminata, we assigned the interactions in the genes B and C to the gene A. Thus, the interactions are mapped according to their orthologs in A. sisalana. In the end, only interactions that have complete pairs with orthologs in A. sisalana were selected. To visualize the PPI network in A. sisalana, the predicted PPI data were loaded into Cytoscape v3.2.1 (Shannon et al., 2003).

2.5 Network analysis

After constructing the A. sisalana interactome and loading it into Cytoscape, the network metrics were analyzed with Network Analyzer v4.4.8. The top hub genes were selected from the intersection of the top 20 hub genes derived from three metrics: betweenness centrality, closeness centrality, and node connect degree, using Cytohubba v0.1 (Chin et al., 2014). The submodules were identified using MCODE v2.0.2 in Cytoscape with a cutoff of 5 in MCODE score.

2.6 Gene ontology (GO) analysis of PPI data

A. sisalana proteins were annotated using PANNZER2 software (Törönen et al., 2018). A table <gene> <go_id> was generated from the results, considering the Gene Ontology class “Biological Process (BP)" in the output of PANNZER2. The enrichment of biological processes was performed using the topGO package v2.48.0 (Alexa et al., 2006) using R 4.1.0.

2.7 KEGG orthology annotation and enrichment

The Kyoto Encyclopedia of Gene and Genomes (KEGG) orthology (KO) annotations were assigned using KofamKOALA 1.0.0 (Aramaki et al., 2020), which uses hmmsearch 3.1b2 against curated hidden Markov model (HMM) KO profiles. The enriched KEGG pathways and their modules have been identified through the enrichKEGG and enrichMKEGG functions of the ClusterProfiler package v4.4.4 (Wu et al., 2021) using R 4.1.0.

2.8 Tissue-specific analysis

Tau metric was used to measure expression specificity among all genes present in the A. sisalana interactome with the Tspex package (Camargo et al., 2020). Tau showed a robust behavior according to data normalization (Kryuchkova-Mostacci and Robinson-Rechavi, 2017), in our case, through the TMM normalization in gene expression data. We considered genes with Tau values greater than 0.90 as tissue specific.

2.9 Weighted correlation network analysis (WGCNA)

Gene expression matrix (normalized gene expression of biological replicates of leaf, stem, and root tissues) was used to construct an unsigned weighted gene co-expression network using the WGCNA R package v1.71 (Langfelder and Horvath, 2008). Hierarchical clustering of samples was conducted to remove outliers with a cut-off value of 80 to produce two stable clusters. Then, the soft threshold power β was determined to ensure a scale-free network. The power function was used to turn the Pearson correlation matrix into an adjacency matrix, which was then transformed into a topological overlap matrix (TOM). Using a dynamic cutting algorithm, hierarchical clustering was performed to cluster similar genes into the same module. Subsequently, we clustered the eigengenes according to the relationship and merged them into modules with an association >0.80. The association of each module with gene expressions in the tissues was evaluated based on Pearson correlation (|cor| > 0.5 and p-value <= 0.05). For each gene, module membership (MM) was characterized according to the association between module eigengene (ME) and its expression level. The association between gene expression and tissues was represented by gene significance (GS). Thus, for each module, GS and MM for each gene were computed and considered significant if |GS| > 0.2 and MM > 0.8.

3 Results and discussion

A detailed outline of our study are summarized in Figure 1. Firstly, we conducted an analysis of orthologous genes between several plant species and the transcriptome assembly of A. sisalana. From the identification of ortholog families, we mapped the interactions of these pairs of orthologs in the STRING database. And then, we predicted the A. sisalana interactome. In addition, we identified the interactome submodules and performed its GO enrichment to search for modules with specific functions in A. sisalana from known databases. In addition, we performed the analysis of hub genes, which play an important role in the interactome.

FIGURE 1

FIGURE 1. Flowchart illustrating the study design. Agave sisalana interactome was predicted from five orthologs organisms based on STRING database interactions in order to identify hub genes. Also, a WGCNA analysis was performed to identify key modules based on module-tissue relationships from three organs of A. sisalana (leaf, stem, and root). Both analyses were integrated to find key genes. The dashed lines indicate software outputs or in-house scripting processes.

Likewise, we performed weighted gene co-expression network analysis from stem, leaf, and root tissues studied in the A. sisalana transcriptome. From this analysis, we can identify the key modules, which are co-expression modules correlated with each tissue, and make hypotheses about the functionalities observed in the plant.

Finally, we cross-referenced the lists of hub genes with their presence in co-expression modules, evaluating module membership (MM) and gene significance (GS). If these hub genes pass through the co-expression module evaluation cut-off, they were considered key genes.

All steps and results present in the flowchart are described in the next sections.

3.1 A. sisalana orthogroups

Initially, a total of 23,794 Agave sisalana proteins (extracted from the transcriptomic atlas) were parsed and grouped into orthologous families using five plant species (A. thaliana, M. acuminata, O. sativa, S. bicolor, and Z. mays), selected by the wide availability of interactome data. Table 1 summarizes the orthogroups containing at least one protein of plant species. From the orthogroups obtained, we selected those that are considered one-to-one or one-to-many. A total of 4,501 A. sisalana proteins present in those orthogroups were eligible for the search for interactions in the STRING database.

TABLE 1

TABLE 1. Overlap of A. sisalana orthologs between two genomes.

3.2 A. sisalana interactome

A. sisalana interactome was reconstructed through the mapping of the ortholog groups and protein interactions extracted from the STRING database. The number of interactions obtained from the STRING database was 772,055 in A. thaliana, 1,384,645 in M. acuminata, 588, 907 in O. sativa, 618,973 in S. bicolor, and 1,368,899 in Zea may. Here, it was assumed that interactions of proteins in one organism are expected to be conserved in other related organisms (Fraser, 2005; Uhrig, 2006).

After predicting PPIs in A. sisalana, the PPI network was visualized using Cytoscape v3.8.1. The final interactome contains 2,582 nodes and 15,083 edges (Supplementary File S1). The topological metrics of the PPI network were calculated using Network Analyzer (Table 2). Analysis of the interaction network showed the short path length distribution, the node degree distribution, the neighborhood connectivity distribution, and the clustering coefficient distribution (Figure 2). The declining neighborhood connection trend reveals a weak clustering coefficient among neighbor nodes with lesser connectedness. On the other hand, the shortest path length distribution showed that the edges have a frequency dominated by a short path length (3–5), as path length means edge. The node degree distribution shows that the A. sisalana interactome has many nodes with a low degree, suggesting this possible short path length. In the clustering coefficient distribution, we note the presence of nodes with a low degree (clustering coefficient equal to zero), but also a group of highly connected genes, with a clustering coefficient equal to 1.

TABLE 2

TABLE 2. Analysis of the interaction network topology of A. sisalana.

FIGURE 2

FIGURE 2. Agave sisalana interactome metrics. (A) short path length distribution (B) the node degree distribution, (C) the neighborhood connectivity distribution, and (D) the clustering coefficient distribution.

At the protein level, important genes are normally strongly interconnected hubs. The top 20 hubs were identified in the A. sisalana interactome (Table “hub genes” in Supplementary File S2) and most of them are functionally related to transcription, phosphorylation, response to stimulus, and nitrate assimilation. The D-fructose-1,6-bisphosphate 1-phosphohydrolase gene (AS_TRINITY_DN57407_c0_g1_i1) was identified as a hub and has an important regulatory role involved in photosynthetic CO₂ assimilation (Chehebar and Wolosiuk, 1980). Another hub is an ATPase, F1 complex, OSCP/delta subunit (AS_TRINITY_DN7474_c0_g1_i1), which has a role in inorganic ion transmembrane transport and participates in many metabolic processes. Protein-serine/threonine phosphatase (AS_TRINITY_DN57726_c0_g1_i1) was identified as a hub in the A. sisalana interactome. For Agave, Protein-serine/threonine phosphatases are related to diel expression patterns and exhibit inverted temporal shifts in abundance when compared to the C3 plant A. thaliana (Abraham et al., 2016).

3.3 GO and KEGG enrichment analysis

We performed GO enrichment analysis of all genes present in the interactome. We identified biological processes related to protein import, protein modification, and methylation, such as methylation (GO:0032259; p-value 5.0e-07), protein peptidyl-prolyl isomerization (GO:0000413; p-value 5.0e-10) and protein import into chloroplast stroma (GO:0045037; p-value 1.7e-04). Basal processes related to the photosystem and chloroplast were also identified, such as photosynthesis (GO:0015979; p-value 1.4e-21), light harvesting (GO:0009765; p-value 5.4e-06), chloroplast fission (GO:0010020; p-value 4.5e-04), chlorophyll biosynthetic process (GO:0015995; p-value 1.1e-07), photosystem I assembly (GO:0048564; p-value 9.8e-04), photosystem II repair (GO:0010206; p-value 9.8e-04) and photosynthetic electron transport in photosystem I (GO:0009773; p-value 5.3e-05).

Regarding KEGG enrichment, some basic pathways detected in plants were enriched, such as Plant hormone signal transduction (ko04075), MAPK signaling pathway - plant (ko04016), Photosynthesis (ko00195), Circadian rhythm - plant (ko04712), Ribosome (ko03010), Biosynthesis of cofactors (ko01240), Plant-pathogen interaction (ko04626), Carbon fixation in photosynthetic organisms (ko00710), Porphyrin metabolism (ko00860), Ubiquinone and another terpenoid-quinone biosynthesis (ko00130), Protein export (ko03060), Biosynthesis of amino acids (ko01230), Protein processing in endoplasmic reticulum (ko04141), Glycerolipid metabolism (ko00561), Glycolysis/Gluconeogenesis (ko00010), Photosynthesis - antenna proteins (ko00196), Carotenoid biosynthesis (ko00906) and Flavonoid biosynthesis (ko00941). In addition, an enriched metabolic pathway of Plant-pathogen interaction (ko04626; p.adjust = 0.00014) was identified, suggesting the presence of a phytopathogen. Recent studies using this same transcriptomic dataset have described several viral (Quintanilha-Peixoto et al., 2021) and fungal (Marone et al., 2022) interactions in A. sisalana, so this metabolic pathway might be a response to these interactions. Furthermore, the flavonoid biosynthesis pathway (ko00941; p.adjust = 0.03169) enrichment and the accumulation of these bioactive compounds in Agave species (Barriada-Bernal et al., 2014; Ahumada-Santos et al., 2013; Morreeuw et al., 2021) may indicate protection against these phytopathogens (Ferreyra. 2012) or may be considered a physiological adaptation to stress conditions (exposure to UV radiation, high temperature, high luminosity (photoprotection) and nutrient deficiency) (Agati et al., 2012, 2013; Taiz and Zeiger, 2003; Yonekura-Sakakibara et al., 2019; Zahedi et al., 2021; Zeng et al., 2021).

3.4 Organ-specificity analysis

Tissue-specificity analysis (Tau metric of 0.90, see methods for details) identified transcripts expressed exclusively in one of the three tissues studied (leaf, stem, and root). From a total of 2,582 genes in the A. sisalana interactome, 2, 18, and 101 were identified in the stem, leaf, and root, respectively (Supplementary File S2). We found a leaf-specific LIM zinc finger protein (AS_TRINITY_DN73400_c0_g1_i1; Tau = 1), which may be involved in plant growth and development, as well as regulating resistance mechanisms to diverse biotic and abiotic stresses as observed in other plants (Gupta et al., 2012).

As a stem-specific transcript, we found an auxin transporter regulating intracellular auxin homeostasis and metabolism (AS_TRINITY_DN46645_c0_g1_i1; Tau = 0.941) which is crucial for plant development (Rosquete et al., 2012). However, auxin has differential distribution (gradients) within plant tissues, thus this gene may be acting as a gatekeeper controlling auxin traffic in and out of plant cells. Finally, among the root-specific genes, we found regulation of jasmonic acid mediated signaling pathway (GO:2000022; p-value: 0.00019), transfer of electrons from cytochrome c to oxygen (GO:0006123; p-value: 0.02372) and photosynthesis (GO:0015979; p-value:0.00022).

3.5 Submodule identification

A total of 16 submodules were identified in the A. sisalana interactome. Each submodule has highly connected genes and can suggest specific processes and roles. The GO enrichment of each submodule (Supplementary File S2) identified two submodules related to methylation and DNA repair/recombination: submodule S9 has 5 genes related to histone methylation (GO:0016571; p-value 8.2e-04), whose expression values are similar in the three tissues, with maximum values ranging from 10 to 20 TPM. Submodule S5 has 21 genes related to DNA repair (GO:0006281; p-value 5e-04), DNA topological change (GO:0006265; p-value 3.7e-06), and DNA recombination (GO:0006310; p-value 5.3e-06), which are all related to DNA damage response. Most genes in submodule 5 have expression values above 10 TPM in all tissues, but genes related to DNA recombination and DNA repair have TPM values from 20 to 142, suggesting these processes were highly demanded by the plant at the moment of sampling. Indeed, the samples were collected at noon and UV light exposure and/or the presence of Reactive Oxygen Species (ROS) can induce the plant response mechanism to fix any damage (Nisa et al., 2019).

We found 3 submodules related to carbohydrate metabolism (S7, S11, and S15). Among them, S11 is the only submodule presenting only biological processes associated with carbohydrate metabolism (GO:0005975; p-value 6.3e-06). Interestingly, the S11 submodule contains the Raffinose synthase gene (AS_TRINITY_DN57582_c0_g2_i1) with TPM values of 28 in the leaf tissue. Raffinose has been previously described as an important carbohydrate in Agave metabolism (Raya et al., 2021) and it is an osmolyte that can act both in drought and oxidative stress (Nishizawa et al., 2008). As for the S7 submodule, we identified GO terms sucrose biosynthetic process (GO:0005986; p-value 7.9e-06) and starch biosynthetic process (GO:0019252; p-value 1.9e-14) as well as terms related to gene regulation, like transcription from plastid promoter (GO:0042793; p-value 2.7e-04), cytidine to uridine editing (GO:0016554; p-value 7.6e-04), and Group II intron splicing (GO:0000373; p-value 5.1e-07). Finally, submodule S15 presented the term regulation of starch biosynthetic process (GO:0010581; p-value 8.8e-05) and several terms related to abiotic stress. Several heat-shock proteins (HSP) were found in this submodule in categories like response to heat (GO:0009408; p-value 1.4e-07), chaperone mediated protein folding requiring cofactor (GO:0051085; p-value 6.5e-04), and protein folding (GO:0006457; p-value 7.9e-06). Among this HSP, Small HSP and DnaJ molecular chaperone were more present. Curiously, the term histone H3 deacetylation was also present in this submodule. The level of histone acetylation has been shown to regulate plant response to drought stress, including ABA-responsive element binding protein (AREB) (Li et al., 2021). AREBs can regulate Dehydration responsive element binding (DREB), which are known for inducing HSP and Heat shock factors (HSF) (Agarwal et al., 2017).

Among the identified submodules, six of them were related to the photosynthesis process and were mostly composed of genes with higher expression in the leaf tissue (S1, S2, S3, S4, S10, and S13), which was expected due to the nature of the organism. Submodule S1 showed the following unique biological processes: nonphotochemical quenching (GO:0010196, p-value 1.4e-05), cellular response to high light intensity (GO:0071486, p-value 0.00024), photosystem II assembly (GO:0010207, p-value 1.4e-05), energy coupled proton transport, down electrochemical gradient (GO:0015985, p-value 2.0e-05), electron transport chain (GO:0022900, p-value 0.00020), starch biosynthetic process (GO:0019252, p-value 0.00014) and regulation of photosynthesis, dark reaction (GO:0010110, p-value 0.00024). The S2 submodule presented genes related to the unique processes protoporphyrinogen IX biosynthetic process (GO:0006782, p-value 0.00045), chloroplast-nucleus signaling pathway (GO:0010019, p-value 0.00058), DNA-templated transcription, termination (GO:0006353, p-value 0.00058) and cytidine to uridine editing (GO:0016554, p-value 0.00071). The biological processes unique to the S3 submodule were photosystem I assembly (GO:0048564, p-value 5.8e-06), protein import into chloroplast stroma (GO:0045037, p-value 0.00045), tetrahydrobiopterin biosynthetic process (GO:0006729, p-value 0.00096) and iron-sulfur cluster assembly (GO:0016226, p-value 0.00044). Submodule 10 has exclusives processes chlorophyll catabolic process (GO:0015996, p-value 2.5e-08), and a terpenoid biosynthetic process (GO:0016114, p-value 2e-04). The S13 submodule presented exclusive processes such as sucrose biosynthetic process (GO:0005986, p-value 2.6e-05), glycolytic process (GO:0006096, p-value 8.0e-05), isoprenoid biosynthetic process (GO:0008299, p-value 0.00036) and regulation of jasmonic acid mediated signaling pathway (GO:2000022, p-value 6.3e-10). By comparing the biological processes in common between these submodules, except for S10, we identified protein peptidyl-prolyl isomerization process, with the highest abundance of peptidylprolyl isomerase (PPIase) expression (AS_TRINITY_DN23174_c0_g1_i1) in leaf with 12 TPM. Interestingly, PPIase is responsible for the proper folding of proteins, being the only enzyme capable of catalyzing the cis-trans transition without energy requirement, unlike chaperones (Fischer et al., 1984; Fischer et al., 1989; Fanghänel and Fischer, 2004). Photosynthesis is a complex process involving several proteins, so the high expression of PPIase might be related to the correct folding of proteins without spending energy (Kirschbaum, 2004; Darko et al., 2014; Martin et al., 2018).

Lastly, submodule S14 has the processes of iron-sulfur (Fe-S) regulation (GO:0016226 and GO:0051176, p-values 4.6e-13 and 0.00038) and pyruvate metabolism (GO:0019450, p-value 0.00038). Plants have a high demand for iron in mitochondria and chloroplasts, mainly to ensure respiration and photosynthesis (Couturier et al., 2013). Iron-Sulphur is required in many metabolic pathways due to the presence of metalloproteins (Johnson et al., 2005). Moreover, Fe-S proteins play a wide range of functionality ranging from radical generation, electron transfer, participation in sulfur and nitrogen assimilation, DNA replication and repair, chlorophyll catabolism, and ribosome biogenesis (Johnson et al., 2005; Balk and Pilon, 2010). Regarding pyruvate metabolism, we found L-cysteine catabolic process to pyruvate (GO:0019450) was reported in this submodule, although the only gene present was L-cysteine desulfhydrase (LCD). LCD is responsible for cysteine degradation and is highly expressed in response to adverse environmental conditions and intracellular oxidative stress (Chen L Y et al., 2020; Wang et al., 2022), which agrees with agaves having a high capacity to be cultivated in places with more adverse environmental conditions.

3.6 Co-expression network analysis

We carried out a hierarchical clustering (WGCNA) to investigate the expression profile as a function of organs (Supplementary Figure S2A). During WGCNA analysis, we found β = 22 and R² = 0.80 as the optimal soft threshold parameters to guarantee a scale-free network (Supplementary File 4: Figure 2S–B). We set the clustering height cut-off to 0.20 to merge similar modules, which resulted in 73 modules (Supplementary Figure S2C). Specifically, modules darkgreen, purple, darkred, darkturquoise, magenta, bisque4, mistyrose, lightyellow, pink, plum1, magenta4, darkolivegreen, blue4, turquoise, darkorange, salmon1, mediumpurple1, brown4, skyblue1, and plum4 were identified as statistically significant (|GS| > 0.2 and MM > 0.8; see methods for details) in at least one organ (Figure 3).

FIGURE 3

FIGURE 3. Heatmap of associations among module eigengenes in organ samples (stem, leaf, and root) from weighted correlation network analysis (WGCNA). A black star marks the key modules, i.e., modules with |cor| > 0.5 and p-value < = 0.05 in at least one organ.

To explore potential genes and pathways associated with each organ, we conducted GO enrichment analysis on the modules with the highest correlation with the organs. (Supplementary File S3). A total of 8 modules showed enriched GO terms: modules darkgreen, magenta, darkred, and purple are positively correlated with root; pink and darkolivegreen are positively correlated with stem, and turquoise and blue4 are positively correlated with leaf. Leaf was the only organ presenting negatively correlated modules: darkorange and darkred (Figure 3). All results of the main GO processes and KEGG pathways associated with each organ of A. sisalana are summarized in Figure 4.

FIGURE 4

FIGURE 4. Co-expression modules and their associated biological processes in the leaf, stem, and root organs of the A. sisalana plant. For each co-expression module, we created summary tables highlighting the main enrichments observed, which included both correlation values and p-values obtained from the module-sample relationship analysis.

The darkgreen module genes (correlation of 0.99 with root expression profile) have enriched GO terms related to translation and cellular iron ion homeostasis (Supplementary Figure S1B). In addition, the magenta module showed enriched GO related to the pentose-phosphate shunt, GDP-mannose transmembrane transport, cytoplasmic translational elongation, and translation. The pentose-phosphate shunt may be acting on root ion transport systems. Studies suggest that the regulation of root nitrogen and sulfur acquisition by plant carbon status is governed by an unnamed oxidative pentose phosphate pathway-dependent sugar detection system, which coordinates the availability of these three elements for amino acid synthesis (Lejay et al., 2008).

The darkred module showed enriched GO terms related to the lignin biosynthetic process, protein localization to the membrane, intracellular auxin transport, S-adenosylmethionine biosynthetic process, and actin filament bundle assembly. The presence of the lignin biosynthetic process positively correlated with the root and negatively correlated with the leaf highlights the importance of lignin for different organs in the plant. This agrees with the low contents of lignin in leaf and stem in agaves (Raya et al., 2021), although there is no available quantification of lignin content in roots. Therefore, roots might possess more lignin to ensure mechanical support (Hoson and Wakabayashi, 2014), act as a defense barrier (Underwood, 2012), and a conduit for long-distance transportation of water and essential nutrients from the roots to the shoots (Boerjan et al., 2003; Zhong and Ye, 2015; Bonawitz and Chapple, 2010; Liu et al., 2018).

In the purple module, we observed the GO process related to response to farnesol, response to chitin, potassium ion transmembrane transport, transmembrane transport, and carbohydrate catabolic process. Here we highlight that the response process to farnesol is positively correlated with the root organ. Previous studies in Arabidopsis have shown that folk flowers, which lack farnesol kinase activity, accumulate farnesol and develop supernumerary carpels under water stress, providing evidence of a molecular link between the farnesol metabolism, abiotic stress signaling, and flower development (Fitzpatrick et al., 2011). Therefore, as the sampled plants were in a dry region, this group of genes was probably signaling the presence of this stress to the plant and activating another warning cascade. Moreover, we observed genes related to ABA signaling in submodule S15, which also shared a role with farnesol kinase in flower development (Fitzpatrick et al., 2011).

Regarding the modules highly correlated with the stem, pink module showed enriched GO terms related to protein phosphorylation, leaf formation, floral organ abscission, mitochondrial mRNA processing, auxin-activated signaling pathway, and regulation of transcription (Supplementary File S3). The process involving auxin signaling is responsible for ramet formation in Agave, involved in the survival and reproduction of a whole plant (Barreto et al., 2010). Another interesting process in the stem is the floral organ abscission. Abscission, in which plants shed unwanted organs, is a natural developmental program or in response to environmental stimuli (Patharkar and Walker, 2015). In A. sisalana, the genes related to abscission in this module, which are positively correlated with stem, could be floral repression genes, which can be part of further improvement studies to accelerate or repress flowering. Agaves are monocarpic, meaning they only flower once at the end of their lifespan, so better control of flowering time can help synchronize it in agave fields. Because agaves accumulate sugar in their stems to be used during flowering, this group of genes may be responsible for this accumulation. The relationship between each gene significance and module membership from the pink module presents a correlation of 0.9884 (Supplementary Figure S1C).

The GO enrichment of darkolivegreen module shows processes related to the regulation of defense response to fungus, brassinosteroid biosynthetic process, and positive regulation of transcription initiation from RNA polymerase II promoter. Indeed, fungal transcripts were described in the same dataset (Marone et al., 2022). Although most fungal transcripts were root-specific, some were present in the stem.

Regarding the modules highly correlated with the leaf, GO analysis of the turquoise module (positively correlated) demonstrated that genes in this module were primarily associated with photosynthesis, photosystem I and II, light harvesting, photosynthetic electron transport, and others (Supplementary File S3). Most of the processes are linked to photosynthesis processes, which makes the turquoise module a candidate module for the reconstruction of the interactome related to leaf genes. The findings of these photosynthesis-related genes revealed the strong linkage of these co-expressed genes to the cited GO terms in A. sisalana leaf. Previous findings in other plant transcriptomes also described that the primary function of light-harvesting complex (Lhc) proteins was related to light absorption by chlorophyll excitation and transfer of absorbed energy to photochemical reaction centers (Zhao et al., 2016; Ma et al., 2018). This relationship can be observed by the significance of the genes in relation to the module membership of the turquoise module, in which they present a correlation of 0.9989 (Supplementary Figure S1A).

The darkorange module is negatively correlated with the leaf (cor = −0.96 and p-value = 0.002). The GO enrichment shows processes related to membrane fusion and brassinosteroid-mediated signaling pathways. The brassinosteroid (BR) class of steroid hormones regulates plant development, growth, and physiology and is involved in mediating adaptation to abiotic stresses, such as drought, temperature changes, and salinity (Planas-Riverola et al., 2019). This result shows that these classes of genes play a key role in maintaining the delicate balance between growth and resilience to environmental threats. These mechanisms of hormonal responses and stress stimuli have already been observed in a previous study of A. sisalana (Sarwar et al., 2019).

Biological processes enriched in blue4 (positively correlated with leaf) are related to the cellular response to hypoxia, anthocyanin-containing compound biosynthetic, RNA biosynthetic process, and auxin-activated signaling pathway (Supplementary File S3). Anthocyanin and hypoxia-related genes may play a relevant role in the tolerance and adaptation of agave plants in stress conditions, such as high temperatures, high light, drought, and low O₂ availability. Genes linked to the anthocyanin-containing compound biosynthetic process may have a role in response to high light incidence (Steyn et al., 2002; Tattini et al., 2014; Petrella et al., 2016; Trojak and Skowron, 2017.; Zheng et al., 2020) since anthocyanins are essential against UV radiation present in light stress conditions (Takahashi et al., 1991; Gould, 2004; Guo et al., 2010; Landi et al., 2015) and mitigate DNA damage in UV-B-irradiated (Takahashi et al., 1991; Kootstra, 1994; Stapleton and Walbot, 1994).

Also, plant responses to oxygen limitation (anoxia and hypoxia) are modulated by common signaling pathways, which target metabolic adaptations (Chang et al., 2011). Under oxygen limitations, tolerant plants increased glycolytic flux or metabolic depression, reducing ATP consumption (Sasidharan et al., 2017). Hypoxia induces a decline in stomatal conductance with a rapid decrease in the rate of photosynthesis (Malik et al., 2001; Araki et al., 2012), however, it was shown that leaves with high CAM activity (crassulacean acid metabolism) have higher liquid photosynthesis under hypoxia (Pedersen et al., 2011). Moreover, CAM plants have their stomata closed during the day and open only at night, when they fix CO₂ (Males and Griffiths, 2017). Hypoxic pre-treated organs have been shown to have the ability to maintain a high glycolytic rate during prolonged periods of anoxia, as well as higher levels of ATP and energy load (Mugnai et al., 2011). There is a close correlation between anoxia tolerance and carbohydrate reserves (Zahra et al., 2021). Carbohydrates, specifically sugars, increase hypoxia tolerance in many species due to their close association as an energy-providing metabolite (Setter et al., 1997; Ram et al., 2002). In addition, drought stress quickly causes an osmotic imbalance (Gurrieri et al., 2020). As the severity of stress increases, plants face drought by accumulating high intracellular levels of osmoprotective compounds (i.e., sugars) to protect cellular components and restore osmotic balance (Singh et al., 2015; Sharma et al., 2019).

3.7 Selection of key genes

For the selection of key genes, we combined the results of interactome with the co-expression network analysis through the selection of hub genes with high significance and module membership in the module-organ relationship. Thus, a total of six hub genes were selected as key genes for participating in four WGCNA modules (turquoise, blue, and blue4) that are positively correlated to leaf, in addition to having a negative correlation with root and stem, revealing that they might also play roles in these two organs. The scoring of each hub gene in PPI and WGCNA is summarized in Table 3. For further validation of these potential key genes, we compared their expression values between organ samples in the dataset. Expression levels of these four key genes in the turquoise module were markedly elevated in leaf samples compared with stem and root samples (Figure 5).

TABLE 3

TABLE 3. Scores of six intersecting hub genes using different ranking methods in PPI and WGCNA.

FIGURE 5

FIGURE 5. Per organ normalized expression of six hub genes identified in PPI of A. sisalana. AS_TRINITY_DN57726_c0_g1_i1 (PSP); AS_TRINITY_DN57407_c0_g1_i1 (FBPase); AS_TRINITY_DN51514_c0_g1_i1 (RPL13); AS_TRINITY_DN7474_c0_g1_i1 (ATPase); AS_TRINITY_DN57016_c0_g1_i1 (PPIase); AS_TRINITY_DN57140_c0_g1_i1 (Lhc).

One of these transcripts is the Protein-serine/threonine phosphatase gene (AS_TRINITY_DN57726_c0_g1_i1). This transcript, which was also identified as a highly connected hub, is a protein phosphatase 2C (PP2C). The ortholog of this gene in Arabidopsis (AtPP2C19, At2g20050) (Schweighofer et al., 2004; Xue et al., 2008) is a possible interacting protein of [Ca²⁺]_cyt-assosiated protein kinase 1 (CAP1, At5G61350) (Yu et al., 2022) and it has been shown to contain the HSE motif, which is a cis-acting element involved in heat stress responsiveness (Xue et al., 2008). In A. thaliana, protein phosphatase 2Cs (PP2Cs) has been shown as a major negative regulator of ABA signaling (Umezawa et al., 2009), however, this gene does not group with genes related to this function (Schweighofer et al., 2004; Xue et al., 2008) and recent evidence indicates that it may be related to photosynthesis (chlorophyll molecule) and acquisition of magnesium ions, since it is co-expressed with magnesium transporters (Heidari et al., 2021). Furthermore, in Agave, the diel expression pattern of PP2C is inverted when compared to A. thaliana, and, for that, it has been proposed as one of the main regulators of stomata activity in CAM photosynthesis (Abraham et al., 2016). Moreover, D-fructose-1,6-bisphosphate1-phosphohydrolase (FBPase) (AS_TRINITY_DN57407_c0_g1_i1) is an important enzyme (EC 3.1.3.11) that acts mainly in the regulation of the Calvin cycle and glycogenesis. These reactions are involved in carbon fixation and sugar metabolism present in the chloroplast stroma and cytosol (Flechner et al., 1999). The ortholog in Arabidopsis encodes a chloroplastic FBPase (cFBP1, At3g54050) also known as HIGH CYCLIC ELECTRON FLOW 1 (HCEF1) involved in the regeneration of ribulose 1,5-bisphosphate (RuBP) and in the starch synthesis pathway (Livingston et al., 2010a; Rojas-González et al., 2015). The Arabidopsis mutant cfbp1/hcef1 shows a dwarf phenotype, chlorotic leaves, a low photosynthesis rate, high sucrose/starch rate (Rojas-González et al., 2015) and displays constitutively high cyclic electron flow (CEF) rates (Livingston et al., 2010b; Strand et al., 2015; Sharkey and Weise, 2016). Peptidylprolyl isomerase (PPIase) (AS_TRINITY_DN57016_c0_g1_i1) (EC 5.2.1.8), which are the only enzymes known that can catalyze cis-trans transition of peptide bonds, essential for the proper folding of proteins (Fischer et al., 1984) and belong to three major classes of proteins: cyclophilins, FK506-binding proteins or FKBPs, and parvulins (Singh et al., 2020). Interestingly, in our analyses, PPIase was reported to be present in all photosynthesis submodules. The ortholog in Arabidopsis encodes CYCLOPHILIN 38 (AtCYP38, At3g01480) and the mutant showed retarded growth and enhanced sensitivity to high light (Fu et al., 2007). In addition, Duan et al. (2021) reported on the importance of CYP38-dependent photosynthetic activity in supporting lateral root emergence. AtCYP38 in the thylakoid lumen functions in the assembly, maintenance of photosystem II (PSII; Fu et al., 2007; Sirpiö et al., 2008; Lepeduš et al., 2009; Vojta et al., 2019) and accumulation of PSII-LHCII supercomplexes (Zhu et al., 2022) but does not have PPIase activity in vitro (Vasudevan et al., 2012). Light-harvesting chlorophyll a/b-binding protein (Lhc) (AS_TRINITY_DN57140_c0_g1_i1) are a class of antennae proteins that play indispensable roles in capturing solar energy as well as photoprotection under stress conditions (Zhao et al., 2020). In Arabidopsis, the ortholog is known as PHOTOSYSTEM B PROTEIN 33 (AtPSB33, At1g71500) (Fristedt et al., 2015). The nucleus-encoded protein AtPSB33 has a crucial role in maintaining the stability of photosystem II complex and regulating photosynthesis, especially under fluctuating, high light and UV-A light stress conditions, and plants lacking PSB33 have a dysfunctional state transition (Fristedt et al., 2015; Fristedt et al., 2017; Kato et al., 2017; Nilsson et al., 2020). Lastly, ATPase, \F1 complex, OSCP/delta subunit (AS_TRINITY_DN7474_c0_g1_i1), which was also selected as a key gene, is present in the blue4 module of the WGCNA analysis. This specific module has a high correlation with the leaf (r = 0.84 and p = 0.04). cpATPase are essential for photosynthesis and their absence causes loss of photoautotrophy and increased photosensitivity. F-type ATP synthase uses the electrochemical proton gradient generated by photosynthesis to produce ATP from ADP and inorganic phosphate (Maiwald et al., 2003; Hahn et al., 2018).

4 Conclusion

In this study, we predicted 2,582 interactome components in A. sisalana using model organism PPI databases, with the numbers of A. sisalana interactome components determined from the numbers of model organism PPIs and A. sisalana orthologues. Also, we used RNA-Seq data from A. sisalana organs to identify co-expressed genes and, consequently, co-expression modules that were positively correlated with leaf, stem, and root organs in an analysis of the module-organ relationship. In the co-expression analysis with WGCNA, we found 72 co-expression modules, of which turquoise, blue, and blue4 were closely associated with leaf organ.

From the co-expression analysis of the module-organ relationship, we were able to identify key modules and their associations with three different organs of A. sisalana (leaf, stem, and root). These associations identified genes that play the role of abscission which can be part of further improvement studies to accelerate or repress flowering. Also, we found anthocyanin and hypoxia-related genes which may play a relevant role in the tolerance and adaptation of Agave plants, and genes with a role in adaptation to stress through brassinosteroid-mediated signaling pathways.

We identified key genes from the interactome network analysis through hub genes identification metrics and performed a correlation with their association in the module-organ relationship. Further analysis suggested that six candidate genes were positively significantly related to leaf organ and negatively related to stem and root, revealing that the key genes in the leaf have a role in maintaining the stability and organization of photosynthetic apparatus during stress conditions.

In conclusion, the predicted PPI network of A. sisalana expands the possibility of comparative analyses with other Agave species, thus providing additional insight into network evolution among species. Furthermore, the identification of co-expression modules and their module-organ relationship, together with the identification of key genes, suggest target genes for specific studies in Agave transcriptome studies, biological processes related to plant survival in unfavorable environments and provide strategies for breeding programs.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

LC: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing-Original draft preparation, Writing-Review and editing, Visualization. NS: Investigation, Data curation, Writing-Original draft preparation, Writing-Review and editing, Visualization. LA: Investigation, Data curation, Writing-Original draft preparation, Writing- Review and editing. MM: Investigation, Data curation, Writing-Original draft preparation, Writing-Review and editing. AC: Software, Writing-Review and editing. FR: Investigation, Data curation, Writing-Original draft preparation, Writing- Review and editing. GA: Writing- Review and editing. MC: Investigation, Methodology, Data curation, Writing-Original draft preparation, Writing- Review and editing. GG: Writing-Original draft preparation, Resources, Writing-Review and editing. All authors contributed to the article and approved the submitted version.

Funding

This work was financed by the Center for Computational Engineering and Sciences–FAPESP/Cepid (2013/08293-7) and the São Paulo Research Foundation (FAPESP) through grants 2019/12914-3 and 2020/02524-0.

Acknowledgments

We would like to thank the São Paulo Research Foundation (FAPESP) and the Center for Computational Engineering and Sciences-FAPESP/Cepid (2013/08293-7) for all infrastructure and financial support. Also, we would like to thank Shell Brasil and ANP (Brazilian National Agency of Petroleum, Natural Gas and Biofuels) for the strategic support provided through the regulatory incentive for Research, Development & Innovation.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fceng.2023.1175235/full#supplementary-material

References

Abraham, P. E., Yin, H., Borland, A. M., Weighill, D., Lim, S. D., De Paoli, H. C., et al. (2016). Transcript, protein and metabolite temporal dynamics in the CAM plant Agave. Plants 2, 16178. doi:10.1038/nplants.2016.178

PubMed Abstract | CrossRef Full Text | Google Scholar

Agarwal, P. K., Gupta, K., Lopato, S., and Agarwal, P. (2017). Dehydration responsive element binding transcription factors and their applications for the engineering of stress tolerance. J. Exp. Bot. 68, 2135–2148. doi:10.1093/jxb/erx118

PubMed Abstract | CrossRef Full Text | Google Scholar

Agati, G., Azzarello, E., Pollastri, S., and Tattini, M. (2012). Flavonoids as antioxidants in plants: location and functional significance. Plant sci. 196, 67–76. doi:10.1016/j.plantsci.2012.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahumada-Santos, Y. P., Montes-Avila, J., de Jesús Uribe-Beltrán, M., Díaz-Camacho, S. P., López-Angulo, G., Vega-Aviña, R., et al. (2013). Chemical characterization, antioxidant and antibacterial activities of six Agave species from Sinaloa, Mexico. Industrial Crops and Products 49, 143–149. doi:10.1016/j.indcrop.2013.04.050