Network Topological Analysis for the Identification of Novel Hubs in Plant Nutrition

Network analysis is a systems biology-oriented approach based on graph theory that has been recently adopted in various fields of life sciences. Starting from mitochondrial proteomes purified from roots of Cucumis sativus plants grown under single or combined iron (Fe) and molybdenum (Mo) starvation, we reconstructed and analyzed at the topological level the protein–protein interaction (PPI) and co-expression networks. Besides formate dehydrogenase (FDH), already known to be involved in Fe and Mo nutrition, other potential mitochondrial hubs of Fe and Mo homeostasis could be identified, such as the voltage-dependent anion channel VDAC4, the beta-cyanoalanine synthase/cysteine synthase CYSC1, the aldehyde dehydrogenase ALDH2B7, and the fumaryl acetoacetate hydrolase. Network topological analysis, applied to plant proteomes profiled in different single or combined nutritional conditions, can therefore assist in identifying novel players involved in multiple homeostatic interactions.


INTRODUCTION
Living organisms are increasingly viewed as integrated and communicating molecular networks, thanks also to the diffusion of data-derived Systems Biology approaches (Barabási and Oltvai, 2004). Such approaches are well established in biomedical and pharmaceutical research (Zhou et al., 2014;Guney et al., 2016) but not widely used in plant science. However, a growing number of plant biologists are adopting them, also to achieve a refined comprehension of plant nutrition (Mai and Bauer, 2016;Di Silvestre et al., 2018;Tsai and Schmidt, 2020).
Plant systems biology studies are dominated by transcriptomic data and statistics that, by measuring the dependence between variables (transcripts), allow us to reconstruct and analyze co-expression network models (Rao and Dixon, 2019). Instead, fewer studies rely on proteinprotein interaction (PPI) networks, mainly due to the lack of accurate plant models (Di Silvestre et al., 2018). Nevertheless, the number of studies focusing on high-throughput profiling of plant proteomes (Vigani et al., 2017;Castano-Duque et al., 2018;Santos et al., 2019), on experimental identification of PPIs (Senkler et al., 2017;Osman et al., 2018;Huo et al., 2020) or their computational prediction (Ding and Kihara, 2019) are recently increasing. The computational prediction of PPI is usually inferred by transferring interactions from model plant orthologs, like Arabidopsis thaliana. This approach represents one of the most important resources for network inference in plant organisms, with some flaws, including the identification of false positives and inadequate coverage from associalogs (Lee et al., 2015).
The co-expression networks reconstructed from proteomic data may represent an alternative to the lack of accurate PPI models and a tool for handling, at the system level, largescale proteomic datasets related to the non-model plant. Similar to transcript networks, protein co-expression networks are defined as an undirected graph where nodes correspond to proteins and edges indicate significant correlation scores. By exploiting the laws underlying graph theory, the identification of hubs/bottlenecks and of differentially correlated proteins are sought, as well as the topological and functional modules related to specific biological phenotypes (Vella et al., 2017).
The studies that elucidate plant iron (Fe) nutrition are steadily increasing due to their impact on alleviating nutritional deficiencies in humans (Murgia et al., 2012(Murgia et al., , 2013. Given these premises, we performed the topological analysis of PPI and protein co-expression networks to identify new proteins involved in Fe and molybdenum (Mo) homeostasis in plants by taking into account subcellular compartmentalization. In particular, we reconstructed and analyzed networks from mitochondrial proteomes purified from roots of Cucumis sativus plants grown under single or combined Fe and Mo starvation (Vigani et al., 2017).

Protein-Protein Interaction (PPI) Networks
A C. sativus PPI network was reconstructed based on homology with A. thaliana by using the STRING v11 database (Szklarczyk et al., 2019). Retrieved homologous interactions were filtered by retaining exclusively those databases annotated and/or experimentally determined with a STRING Score of >0.15 and >0.35, respectively. Following these parameters, a fully connected network of 903 nodes and 10456 edges was built. Four different sub-networks were then extracted by considering the proteins identified in each condition analyzed: +Mo+Fe, 463 nodes and 4118 edges; −Mo+Fe, 466 nodes and 4064 edges; +Mo−Fe, 480 nodes and 4682 edges; and −Mo−Fe, 474 nodes and 4425 edges.

Functional Modules in C. sativus PPI Network
A set of proteins differentially expressed in C. sativus root mitochondria under different conditions (+Mo+Fe, +Mo−Fe, −Mo+Fe, and −Mo−Fe) were selected by Linear Discriminant Analysis (LDA, P < 0.05 and F ratio > 3), as previously reported (Vigani et al., 2017). Starting from this set of differentially expressed proteins, a C. sativus PPI network was reconstructed by homology with A. thaliana using STRING v11 database, as reported above. Proteins were grouped in functional modules by the support of BINGO (Maere et al., 2005) and STRING (Doncheva et al., 2019). Cytoscape's Apps; the network was visualized by the Cytoscape platform, and node color code indicates upregulated (red) and downregulated (light blue) proteins based on Spectral count (SpC) normalization (SpC normalized in the range 0-100, by setting to 100 the higher SpC value per protein).

Reconstruction of C. sativus Protein Co-expression Networks
Protein co-expression networks were reconstructed by processing C. sativus root mitochondria protein profiles characterized under different conditions: +Fe, n = 12; +Mo, n = 14, −Fe, n = 15; −Mo, n = 13. Only proteins identified at least in 51% of the considered samples were retained and processed for each condition. Protein data matrices were processed using Spearman's rank correlation coefficient; Spearman's rank correlation score >|0.7| and P < 0.01 were set as thresholds. All processing was performed using the statistical software JMP15.2, while networks were visualized and analyzed by the Cytoscape platform and its plugins (Su et al., 2014).

Topological Analysis of C. sativus PPI and Protein Co-expression Networks
Both PPI and co-expression networks were topologically analyzed by Centiscape Cytoscape's App (Scardoni et al., 2014). As for PPI networks, Betweenness and Centroid centralities were calculated, and nodes with values above the average calculated on the whole network were considered hubs (Di Silvestre et al., 2018). On the contrary, a set of differentially correlated proteins were selected in the co-expression network by filtering their Degree centrality.

PPI Networks
The analysis of root mitochondria of C. sativus plants grown under Mo and/or Fe starvation led to the identification of 1419 proteins (Vigani et al., 2017). Starting from this proteome, a C. sativus PPI network, connecting 903 nodes through 10464 edges, was built by exploiting the interactions between orthologs in A. thaliana; the average homology between C. sativus and A. thaliana was around 67% (Supplementary Figure 1). A subnetwork made of 118 nodes and 388 interactions was extracted by considering proteins (n = 134) whose level was significantly altered in at least one of the four nutritional conditions (control, single or combined Mo and Fe starvation) (Vigani et al., 2017). Proteins belonging to this sub-network were grouped in 21 functional modules, including redox homeostasis, amino acid metabolism, aldehyde dehydrogenases, protein folding, and electron transport chain (Supplementary Figure 2). The analysis of the behavior of their members under the various nutritional conditions indicated above can support the identification of novel candidate genes; in this paragraph, we will highlight some interesting examples.
Fe starvation, combined with Mo starvation or sufficiency, upregulates several proteins involved in amino acid metabolism, such as glutamate dehydrogenases GDH1 and GDH2, aminobutyric acid transaminase POP2 involved in the metabolism of γ-amino butyrate (GABA), which links C and N metabolism in mitochondria (Shelp et al., 2012), aspartate aminotransferase ASP1, the two subunits MCCA and MCCB of methyl crotonyl-CoA carboxylase involved in leucine catabolism, and beta-cyanoalanine synthase/cysteine synthase CYSC1, which detoxifies the cyanide produced as by-product of ethylene synthesis (Hatzfeld et al., 2000;Yamaguchi et al., 2000;Garcia et al., 2010;Arenas-Alfonseca et al., 2018;Supplementary Figure 2). Notably, ethylene is a well-established regulator of Fe deficiency responses (Romera et al., 1999;Lucena et al., 2015). An increased concentration of amino acids has been documented in the xylem sap of Fe deficient cucumber plants (Borlotti et al., 2012). These results suggest that, under Fe starvation, a larger part of the amino acid metabolism takes place in the radical apparatus, with respect to the aerial parts of the plants.
The redox homeostasis group includes formate dehydrogenase (FDH), the enzyme catalyzing the oxidation of formate into CO 2 (Alekseeva et al., 2011) and involved in Fe and Mo homeostasis (Vigani et al., 2017;Murgia et al., 2020). This group also includes various proteins upregulated under Fe starvation, such as glutathione peroxidase GPX6, alternative oxidase AOX2, glyoxalase II3/sulfur dioxygenase GLY3, At5g42150 coding for another glutathione peroxidase and aldehyde dehydrogenase ALDH2B7 (Supplementary Figure 2). The aldehyde dehydrogenase enzymes can act as "aldehyde scavengers" of reactive aldehydes produced during the oxidative degradation of membrane lipids (Brocker et al., 2013). In particular, ALDH2B7 converts acetaldehyde into acetic acid (Rasheed et al., 2018). The redox homeostasis group also includes SAG21, which is strongly upregulated under single Mo starvation (Supplementary Figure 2). Perturbation of SAG21 expression affects biomass, flowering, the onset of leaf senescence, the pathogens' ability to proliferate (Salleh et al., 2012) and alters the function or stability of mitochondrial proteins involved in ROS production and/or signaling (Salleh et al., 2012).
The fatty acid metabolism group includes dihydrolipoamide branched chain acyltransferase BCE2, also known as Dark Inducible 3 (DIN3), whose expression increases soon after exposure to darkness (Fujiki et al., 2000(Fujiki et al., , 2005. The perturbing effects of prolonged darkness on cellular Fe status and its impact on the expression of Fe homeostasis genes are well known: as an example, the light/dark circadian cycles as well as prolonged darkness influence the expression of the iron-storage ferritin, whose gene expression is activated by prolonged darkness (Tarantino et al., 2003) and is circadianregulated (Duc et al., 2009;Hong et al., 2013). BCE2 is upregulated under Fe starvation (Supplementary Figure 2); this finding is intriguing since darkness is associated with a rise in free Fe ions due to the dismantling of the photosynthetic structures. It would be interesting to investigate BCE2 regulation under Fe excess to establish if BCE2 is upregulated by Fe starvation only or, more generally, by Fe fluctuations.
The protein folding group includes the brassinosteroid (BR) biosynthesis gene DWF1, which mediates BR biosynthesis during the positive growth responses of the root system to low nitrogen (Jia et al., 2020). DWF1 is downregulated under Fe starvation and upregulated under single Mo starvation (Supplementary Figure 2).
The cell cycle/division group includes the mitochondrial GTPase MIRO1 (Supplementary Figure 2), which regulates mitochondrial trafficking and shape in eukaryotic cells. In particular, MIRO1 GTPase influences mitochondrial morphology in A. thaliana pollen tubes (Yamaoka and Leaver, 2008;Sørmo et al., 2011); interacting partners of MIRO GTPases are sought to better elucidate their functions in mitochondrial dynamics (Yamaoka and Hara-Nishimura, 2014). Notably, MIRO1 is upregulated under single Fe or Mo starvation but not under combined deficiencies (Supplementary Figure 2).
The translation group includes the pentatricopeptide protein PPR336, which associates with mitochondrial ribosomes (Waltz et al., 2019); PPR336 is downregulated under single Mo starvation, whereas it is upregulated under double Fe and Mo starvation (Supplementary Figure 2).
The electron transport chain group includes the mitochondrial carrier UCP1/PUMP1 transporting aspartate out /glutamate in (as well as other amino acids, dicarboxylates, phosphate, sulfate, and thiosulfate) across the mitochondrial membrane (Monne et al., 2018); such carrier is upregulated in the single or combined Fe or Mo starvation (Supplementary Figure 2).
Transmembrane transporter activity group includes OM47 related to the voltage-dependent anion channel VDAC family; members of this family are major components of the outer mitochondrial membrane and are involved in the channeling of the products of chloroplast breakdown into the mitochondrion and in the exchange of various compounds between the cytosol and the mitochondrial intermembrane space (Li et al., 2016). OM47 is upregulated under single Fe starvation (Supplementary Figure 2). The role of VDAC protein family also emerged by the topological analysis of PPI network. Indeed, 91 PPI hubs were selected in these conditions: 27 hubs in control (+Mo+Fe), 17 hubs in +Mo−Fe, 26 hubs in −Mo+Fe, and 21 hubs in −Mo−Fe (Supplementary Table 1), while 28 PPI hubs were specifically related to sufficiency or starvation of either Fe or Mo (5 hubs in +Fe, 10 hubs in −Fe, 2 hubs in +Mo, and 11 hubs in −Mo) (Supplementary Table 2). In this scenario, VDAC4 and VDAC2 turned out as hubs in PPI networks under Mo starvation (Supplementary Table 2).
The proteins downregulated under Fe or Mo starvation can represent interesting homeostatic regulators in the PPI networks, and their occurrence should not be neglected in future, more extensive analyses; as an example, rice OsIRO3 plays an important role in the Fe deficiency response by negatively   Protein selection was based on the differential node Degree: Proteins with a Degree above the Network Average Degree (in bold) were retained. For each selected C. sativus protein, the C. sativus gene name, the name of A. thaliana homologous gene and the corresponding homology (in percent value and score), Linear Discriminant Analysis (LDA) parameters (only for differentially expressed proteins), and node Degree in the various nutritional conditions are reported. Homologies were retrieved by STRING database. Co-expression hubs in ( regulating the OsNAS3 expression and, thus, nicotianamine NA levels (Wang et al., 2020). . Around 5 to 7% of the proteins that were retrieved as co-expressed were also physically interacting (Supplementary Table 3). Following the network topological analysis, proteins belonging to the bioenergetics and amino acids metabolisms were found topologically relevant under Fe sufficiency (18 protein hubs) ( Table 1A) or Mo sufficiency (12 protein hubs) (Table 1B). Moreover, 15 proteins emerged as network hubs under Fe starvation (Table 1C) and 16 proteins as hubs under Mo starvation (Table 1D); as examples, we then reconstructed the interactions for three hubs under Mo starvation, i.e. FDH (Table 1D and Supplementary Figure 3), ALDH2B7 (Table 1D and Supplementary Figure 4), CYSC1 (Table 1D and Supplementary Figure 5), and two hubs under Fe starvation, the Voltage-Dependent Ion Channel VDAC4 (Table 1C and Figure 1) and fumaryl acetoacetate hydrolase At4g15940 (Table 1C and Supplementary Figure 6). Formate dehydrogenase has been associated with stress responses in plants (Hourton-Cabassa et al., 1998;Shiraishi et al., 2000;Herman et al., 2002;Lou et al., 2016;Murgia et al., 2020), and it takes part in Fe and Mo homeostasis (Vigani et al., 2017;Murgia et al., 2020).
ALDH2B7 and CYSC1 are themselves hubs under Mo starvation (Table 1); indeed, under this condition, 21 different proteins are co-expressed with ALDH2B7, whereas, under the other conditions, a total of 10 proteins are co-expressed with ALDH2B7 ( Table 1). The majority of these are involved in the synthesis/utilization of acetyl CoA metabolic intermediate, such as (i) MCCB (3-methylcrotonyl-CoA carboxylase) involved in the branched chain amino acids metabolism and (ii) ATCS (citrate synthase) and SDH1-1 (succinate dehydrogenase) involved in TCA cycle (Supplementary Figure 4). The modulation of ALDH2B7 expression observed under stress might play a role in the aerobic detoxification of acetaldehyde in plants (Rasheed et al., 2018). The presence of several co-expressed proteins involved in amino acids and GABA metabolisms (GDH1, POP2), in formate metabolism (FDH), and in the nucleotide metabolism (FAC1, adenosine 5 monophosphate deaminase) suggests that such acetaldehyde detoxification might be targeted to the synthesis of precursor for the synthesis of nucleotides. Plants possess metabolic pathways for the de novo synthesis of purine nucleotides producing AMP as well as pyrimidine nucleotides producing UMP (Witte and Herde, 2020.) Nucleotides are synthesized from amino acids like Glutamine, Aspartic acid, Glycine, and formyl tetrahydrofolate; the latter is a metabolite related to FDH and formate metabolism (Witte and Herde, 2020).
CYSC1 is co-expressed, under Mo starvation, with 16 proteins whereas it is co-expressed with 11 proteins, under Mo sufficiency (Table 1); interestingly, two glutamate dehydrogenase isoforms, GDH1 and GDH2, are co-expressed with CYSC1, under Mo starvation and Mo sufficiency, respectively, but not under Fe sufficiency or Fe starvation (Supplementary Figure 5). GDHs are enzymes at the branch point between carbon and nitrogen metabolism; both isoforms are localized in mitochondria and GDH2 unless GDH1 is calcium-stimulated (Grzechowiak et al., 2020).
AC4 interacts with tRNAs, and it might be involved in their transport into mitochondria (Hemono et al., 2020); accordingly, atvdac4 mutant shows severely compromised FIGURE 1 | Differential VDAC4 correlation degree in co-expression networks. Cucumis sativus root mitochondria proteins identified under different conditions (+Mo+Fe, +Mo-Fe, -Mo+Fe, and -Mo-Fe) were analyzed according to their grouping into Fe sufficiency, Mo sufficiency, Fe starvation, and Mo starvation. Red and blue edges indicated positive and negative correlations, respectively. Edges width is proportional to Spearman's correlation score (P < 0.01).
growth (Hemono et al., 2020). Under Fe starvation, VDAC4 is a hub of two pentatricopeptide-repeat-containing proteins, i.e., PPR596 and At1g60770 (Figure 1). PPR proteins are nuclear-encoded RNA-binding proteins containing repeated motives of approximately 35 aa (Doniwa et al., 2010). PPR596 is involved in the partial C-to-U RNA editing events in the mitochondrial genome and, possibly, also in transcript stabilization and stimulation of translation (Doniwa et al., 2010). Both PPR596 and another PPR protein encoded by At1g60770 also belonging to the VDAC4 co-expression network under Fe starvation (Figure 1) are the closest homologs of the mitochondrial PPR protein PACO1, which affects the flowering time in A. thaliana (Emami and Kempken, 2019). VDAC4 is also the hub, under Fe starvation, of L-galactono 1-4 lactone dehydrogenase GLDH, the enzyme that catalyzes the last step in the biosynthesis of ascorbic acid in plants (Wheeler et al., 2015).
VDAC4 is a co-expressed protein of the fumaryl acetoacetate hydrolase (Supplementary Figure 6) and both of them are also hubs in +Fe−Mo and −Fe−Mo PPI networks (Supplementary Table 2). The fumaryl acetoacetate hydrolase catalyzes the formation of fumarate and acetoacetate in the L-tyrosine (Tyr) degradation pathway in plants (Dixon and Edwards, 2006). Altogether, 48 proteins are co-expressed with the fumaryl acetoacetate hydrolase under Fe starvation but none in the other conditions.

CONCLUSION
We hereby show how the topological network analysis applied to proteomes obtained from mitochondria purified from plants grown under Fe and/or Mo starvation suggests FDH as a hub of Mo nutrition in agreement with experimental observations (Vigani et al., 2017;Murgia et al., 2020). Such analysis could also suggest other potential hubs for Fe and Mo nutrition, such as VDAC4, CYSC1, ALDH2B7, and fumaryl acetoacetate hydrolase among others, which will be object of our future investigations. VDAC4 and fumaryl acetoacetate hydrolase were found to be hubs by both PPI and protein co-expression networks. The different origin of these network models (basically PPI reconstructed from literature, while co-expression networks from experimental data) strengthens the hypothesis that these two proteins may play a role in Fe and Mo homeostasis, and it confirms the complementarity and synergy of these two approaches in identifying candidate proteins with a role in the processes underlying the investigated phenotypes.
Although the approaches based on computational prediction present some intrinsic limitations, including false-positive interactions or the lack of true ones, they nevertheless promote and support new experimental avenues, including large-scale experimental identification of PPIs, which will improve the effectiveness and accuracy of the proposed approaches.

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
DD, PMo, GV, and IM conceived the work. DD reconstructed the PPI and co-expression networks with contributions from PMa and SH. IM, GV, and PMo analyzed all the networks and their biological relevance. IM wrote the manuscript with contributions from DD, GV, and PMo. All authors contributed to the article and approved the submitted version.