{'fraction_non_alphanumeric': 0.040000823485815455, 'fraction_numerical': 0.04181249227982048, 'mean_word_length': 5.292913589843244, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computational Systems Biology of Alfalfa – Bacterial Blight Host-Pathogen Interactions: Uncovering the Complex Molecular Networks for Developing Durable Disease Resistant Crop

Medicago sativa (also known as alfalfa), a forage legume, is widely cultivated due to its high yield and high-value hay crop production. Infectious diseases are a major threat to the crops, owing to huge economic losses to the agriculture industry, worldwide. The protein-protein interactions (PPIs) between the pathogens and their hosts play a critical role in understanding the molecular basis of pathogenesis. Pseudomonas syringae pv. syringae ALF3 suppresses the plant’s innate immune response by secreting type III effector proteins into the host cell, causing bacterial stem blight in alfalfa. The alfalfa-P. syringae system has little information available for PPIs. Thus, to understand the infection mechanism, we elucidated the genome-scale host-pathogen interactions (HPIs) between alfalfa and P. syringae using two computational approaches: interolog-based and domain-based method. A total of ∼14 M putative PPIs were predicted between 50,629 alfalfa proteins and 2,932 P. syringae proteins by combining these approaches. Additionally, ∼0.7 M consensus PPIs were also predicted. The functional analysis revealed that P. syringae proteins are highly involved in nucleotide binding activity (GO:0000166), intracellular organelle (GO:0043229), and translation (GO:0006412) while alfalfa proteins are involved in cellular response to chemical stimulus (GO:0070887), oxidoreductase activity (GO:0016614), and Golgi apparatus (GO:0005794). According to subcellular localization predictions, most of the pathogen proteins targeted host proteins within the cytoplasm and nucleus. In addition, we discovered a slew of new virulence effectors in the predicted HPIs. The current research describes an integrated approach for deciphering genome-scale host-pathogen PPIs between alfalfa and P. syringae, allowing the researchers to better understand the pathogen’s infection mechanism and develop pathogen-resistant lines.


INTRODUCTION
Alfalfa (Medicago sativa L.) or lucerne, a member of the Fabaceae family, is a common forage legume in the United States and other countries. Because of its broad adaptation, high biomass production, and perennial nature, it is known as the "Queen of forages." It also aids in soil and water conservation, biological nitrogen fixation, and insect pest interruption in crop rotations, making it an important component of sustainable agriculture (Nemchinov et al., 2017). The crop is widely cultivated in the regions of hot, dry summers and well-drained soils, and thus is successfully grown in similar environments of Asia Minor, Iran, Southern Europe, Mexico, and the Middle East (Putnam and Orloff, 2014). As per the reports of the USDA National Agricultural Statistics Service (NASS), the estimated value of alfalfa hay and haylage produced in the United States counts to $10.8 billion per year, with an area of 11.7 million acres under production and an average yield of 3.61 tons per acre in 2019 1 .
Many diseases affect alfalfa, but a recent outbreak of bacterial stem blight in the western and central United States has caused concern. There have also been cases of the disease recorded in Europe, Australia, and western Iran (Harighi, 2007). The disease, caused by Pseudomonas viridiflava, was first reported in California and Utah in April and May of 2016 and 2017, with the symptoms highly similar to those caused by Pseudomonas syringae pv. syringae, which causes a significant loss to crop quantity and quality (Lipps et al., 2019). P. syringae pv. syringae, a ubiquitous epiphyte, is a major disease-causing pathogen in a wide variety of cultivated plant species. It is a rod-shaped gramnegative bacterium with a well-sequenced genome. P. syringae species have a broad variety of virulence factors, such as a type III secretion system (T3SS), ice nucleation activity, toxic substances, cell wall degrading enzymes, and exopolysaccharides, which makes it a good model for better understanding the pathogen-host interactions (Morris et al., 2013). The strains of P. syringae produce a phytotoxin, syringomycin, which is thought to play a role in the virulence of the pathogen. Most strains are known to function as ice nuclei, causing frost damage to plants when temperatures drop below zero (Agrios, 2005). The disease infection occurs in two phases-localized foliar necrosis (blight) and systemic vascular wilt. The bacterium infects the host stem mainly at frost injury sites, where it penetrates and causes watersoaked lesions, as well as the emergence of spindly stems that blacken with age. Depending on the incidence rate, the disease can result in yield losses of up to 50% or more (Race et al., 2006;Lamichhane et al., 2015). The outer membrane (Lindow et al., 1989) of P. syringae is surrounded by a unique protein, icenucleation protein (INP), that mimics the crystalline structure of ice and thus serves as an initiator for ice formation, linking bacterial stem blight to frost. This protein is of high interest to the researchers not only because of its pathogenicity, but also because of its other potential applications such as frozen food preparation and snowmaking. Scientists also discovered that the disease-causing strain of P. syringae is a weak pathogen of crops such as sugar beet and snap bean, and a few genes that are unique 1 https://www.nass.usda.gov/ to the alfalfa pathogen have been identified (Li Q. et al., 2012;Proud sponsor of Midwest Forage Association, 2017).
Most of the disease-causing mechanisms include proteinprotein interactions (PPIs), which play an important role in the infection process as well as in initiating the defense responses against the disease. Therefore, studying the PPI network between the pathogen and plant proteins is a crucial step for understanding the underlying mechanism of the infection . Genomic sequencing reveals that nearly 3,000 proteins directly interact with the pathogen proteins or are involved in plant defense . On a genome-wide scale, the computational approaches reveal the relationships between predicted proteins. Recently, a wide range of computational methods for predicting hostpathogen interactions (HPIs) have been developed, based on diverse data types or properties such as protein sequence similarity (Matthews et al., 2001;Sun et al., 2017), protein domain interactions (Ng et al., 2003;Binny Priya et al., 2013), protein structural information (Keskin et al., 2008), and gene ontology (GO) annotations (Wu et al., 2006;Zhong and Rajapakse, 2020). Among the available computational methods, the interolog and domain-based methods (Shoemaker and Panchenko, 2007) are most used.
In the present study, we deciphered genome-wide PPIs between Medicago sativa and Pseudomonas syringae pv. syringae ALF3 using interolog (homology-based) and domain-based methods. We discovered that functional characterization of the HPIs can reveal the molecular mechanisms of pathogen infection. Furthermore, we identified novel protein hubs, enriched molecular function, biological processes, and pathways by conducting extensive functional annotation of the predicted interactome, which could be crucial in fully understanding disease infection mechanism. We also predicted the localization of P. syringae pv. syringae ALF3 protein in alfalfa cells. We assume that the predicted virulence factors described in this study will serve as a solid foundation for future experimental validations and provide a deeper understanding of pathogen infection in alfalfa.

MATERIALS AND METHODS
A detailed pipeline of the computational prediction and functional analysis of HPIs is depicted in Figure 1.

Data Collection
The whole proteomes of M. sativa and P. syringae pv. syringae strain ALF3 were downloaded from different sources. M. sativa proteome was downloaded from LegumeIP v3 (Dai et al., 2020), consisting of 87,892 protein sequences, while proteome of P. syringae was downloaded from JGI IMG/M (Taxon ID: 2609460275) and UniProt database (UP000028706) 2 with a count of 4,983 and 4,842 protein sequences, respectively. The redundant protein sequences were removed from the proteomes using CD-HIT (Fu et al., 2012) at 100% identity. In the case of P. syringae, the proteomes from both the sources were merged, followed by the removal of duplicate sequences.
Proteins whose subcellular position was predicted to be a cytoplasmic membrane or cytoplasm membrane were excluded from our analysis because they are not thought to be involved in host-pathogen interactions. Other proteins classified as periplasmic, extracellular, outer membrane, or unknown were considered to be positive candidates for interactions. Following that, the entire P. syringae proteome was analyzed in the EffectiveDB 3 (Jehl et al., 2011) to predict the proteins that are labeled as secreted for the proteins predicted as cytoplasmic. These proteins were also thought to be promising candidates for pathogen-host interactions. After removing the redundancy and cytoplasmic proteins in the proteomes, HPIs were identified using 87,156 alfalfa and 4,427 P. syringae protein sequences.
To implement the domain-based prediction, three domaindomain interaction (DDI) databases, namely 3did v2020_01 (three-dimensional interacting domains) (Mosca et al., 2014), 3 https://effectors.csb.univie.ac.at/ DOMINE v2.0 (Database of Protein Domain Interactions) (Raghavachari et al., 2008), and IDDI v2011.05.16 (Integrated Domain-Domain Interaction) (Kim et al., 2012) were downloaded and stored in individual files. For analysis, the Pfam v31.0 (Finn et al., 2014) database was downloaded to predict the domains for host and pathogen proteins. Further, using the Pfam database, the identified domains for the host and pathogen species were queried against the three databases (3did, DOMINE and IDDI) to search for the domain-domain interactions. Detailed information about the databases, along with the number of sequences and interactions is available in Supplementary Material 1, excel sheet 1.

Interolog-Based Prediction
Interolog is the method of determining the conserved interactions between two proteins based on their sequence similarity (Yu et al., 2004). When two interacting proteins X and Y in one species have interacting orthologs X' and Y' in another, the interaction pairs X-Y and X'-Y' are referred to as interologs. The proteomes of M. sativa and P. syringae were aligned against six protein-protein interaction databases, viz. IntAct, MINT, DIP, HPIDB, BioGRID, STRING, using BLASTp v2.7.1 with default parameters (e-value, sequence identity, sequence coverage) to make the prediction. An optimal combination of parameters to filter the BLAST alignments of M. sativa and P. syringae was determined using e-value (1e-04, 1e-05, 1e-10, 1e-20, 1e-25, 1e-30, 1e-50), sequence identity (30, 40, 50, and 60%), and sequence coverage (30, 40, 50, 60, and 80%). However, no gold standard cutoff criteria for the prediction of HPIs is reported. In the past, researchers predicted the cassava protein interactome using the interolog-based approach with at least 60% identity, 80% sequence overlap, and an e-value lower than 1e-10 (Thanasomboon et al., 2017). In another study on Human-Burkholderia pseudomallei , they used an e-value of 1e-50, 40% sequence identity, and 80% sequence coverage to decode the host-pathogen relationship. The investigation on Arabidopsis-Pseudomonas syringae (Sahu et al., 2014) for the prediction of plant-pathogen interactions used the e-value of 1e-04, sequence identity of 50%, and sequence coverage of 80%. In our study, we discovered that a sequence identity of at least 30%, an e-value of less than 1e-10, and sequence coverage of at least 80% was the best combination of alignment parameters for predicting HPIs.

Domain-Based Prediction
The domain-based approach takes into account the protein domain profiles of known intra-species PPIs, allowing it to predict host-pathogen protein-protein interactions (Dyer et al., 2007). The prediction in domain-based methods is based on protein structural knowledge found in various domain-domain interaction (DDI) databases. The host and pathogen domains were inferred from the Pfam database using HMMER v3.3.1 (Eddy, 2011), which were then used to predict PPIs using inhouse SQL queries on the local machine. For M. sativa, the results of hmmscan were filtered using an e-value and coverage of 1e-23 and 0.2, respectively, while for P. syringae, an e-value and coverage of 1e-18 and 0.35 were used. Hmmscan has been successfully used in a number of studies to infer host-pathogen interactions (Mondal et al., 2017;Cuesta-Astroz et al., 2019;Lian et al., 2020).

Integration of Interolog-and Domain-Based Approaches
After separately implementing interolog and domain methods as described above, the results from both the approaches were merged, and duplicate predictions were eliminated. Since there were a high number of predicted PPIs (more than 10 million), we considered the consensus predictions of both approaches, which reduces the chances of false-positive HPIs. A similar approach of consensus interactions has been considered previously while predicting HPIs from multiple interaction databases .

Subcellular Localization of Proteins Involved in Host-Pathogen Interactions
Plant pathogenic bacteria relies on type III secretion systems to infect host cells, thus affecting a variety of host cell mechanisms (Aung et al., 2017;Xue et al., 2019). So, to understand the infection mechanism of the pathogen, subcellular localization analysis is an essential step in the study. We evaluated various bioinformatics methods to determine the subcellular position of the proteins, and an appropriate benchmark was chosen for the study, based on the maximum positive matches with the subcellular localization annotation at UniProt 4 . We used PSORTb 3.0.2 (Yu et al., 2010), a commonly used tool for bacterial protein localization, to determine the subcellular localization of P. syringae proteins. To determine the subcellular localization of alfalfa proteins, a homology-based and machine learning model of the tool Plant-mSubP (Sahu et al., 2020) was implemented with default parameters.

Dataset Collection for Known Effectors in P. syringae for Validation
The effectors are deployed by the pathogenic bacteria to inhibit the plant immune response and disrupt the host cell mechanisms (Pelgrom et al., 2020). In this regard, we extracted effector proteins, T3SS-hop proteins, and virulence proteins from EffectiveDB (see text footnote 3), Pseudomonas syringae Genome Resources 5 , and Pseudomonas genome DB 6 , respectively. The following criteria was followed to obtain effector proteins from the above-mentioned sources.

Potential Effector Proteins
EffectiveDB was used to predict the possible effector proteins of P. syringae, using P. syringae pv. syringae B728 as a reference. EffectiveDB is a database of bacteria predicted secreted proteins that includes five methods for effector prediction (EffectiveT3, EffectiveELD, EffectiveCCBD, Predotar, and T4SEpre). EffectiveELD, EffectiveCCBD, and EffectiveT3 proteins were chosen because these methods predict "type III secreted" proteins. A total of 641 effector proteins were predicted.

Type III Secretion System-Hop Proteins
The known T3SS effector proteins belonging to P. syringae were obtained from Pseudomonas syringae Genome Resources. In literature (Dillon et al., 2019), it has been reported that there are 66 effector proteins, but we found only 26 effectors which were further processed.

Virulence Proteins
The virulence factor proteins were downloaded from Pseudomonas genome DB using Pseudomonas aeruginosa PAO1 as reference. The orthologs were searched and eight virulence factors were obtained for P. syringae.
In total, we obtained 655 unique effectors (known, potential, and virulence factors) combined from the above-mentioned sources. These effectors were also used in the selection of an appropriate combination of parameters for filtering BLAST alignments.

Dataset Collection for Ice-Nucleation Proteins for Validation
Ice-nucleation proteins are associated with bacterial stem blight as these are thought to be the initiator of the frost mechanism. The gene, Psyr1608, located on the outer membrane of the bacteria, belongs to Pseudomonas syringae pv. syringae B728a (Feil et al., 2005). This gene is known to contain the octamer repeat of the bacterial ice-nucleation proteins, and is thus involved in the disease. Taking this into consideration, we found the INPs in the projected HPIs. We obtained 30 orthologs of the INPs from Pseudomonas genome DB (see text footnote 6), which were then BLASTed against the P. syringae ALF3 protein dataset. We filtered the BLAST alignments at 40% identity to achieve high performance, which reduced the number of INPs to 17. These were then identified in the predicted HPIs.

Functional Enrichment Analysis
In order to classify the proteins in HPIs that have a common function/biological pathway, it is necessary to obtain the functional annotation of the proteins involved in the interactions. The R package clusterProfiler (Yu et al., 2012) was used to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. Benjamini and Hochberg test correction method (Benjamini and Hochberg, 1995) was used to calculate the adjusted p-values, and enriched GO terms were filtered based on adjusted p-value cutoff of ≤0.05. All the GO annotation ontologies (cellular component, molecular function, and biological processes) were used for gene ontology enrichment analysis. A GO database for alfalfa proteins was created using the makeOrgPackage function in the R package "AnnotationForge, " and GO enrichment analysis of alfalfa proteins was performed, while P. syringae proteins GO terms were enriched with org.Psyringae.eg.db Bioconductor package. KEGG pathway enrichment for alfalfa and P. syringae proteins was performed using enrichKEGG function in ClusterProfiler package at a p-value cutoff of 0.05.

Visualization of Interaction Networks and Identification of Protein Hubs
The study of protein-protein interaction network help in a deeper understanding of various cellular mechanisms involved in the process. To visualize the PPI networks, Gephi v0.9.2 (Bastian et al., 2009) was used in combination with various layout algorithms such as Hu (2006), OpenOrd (Martin et al., 2011), and others to improve the readability of the network (Pavlopoulos et al., 2017). We also analyzed the network on different parameters such as betweenness, closeness, and degree. Functional enrichment analysis of combined interactions is present in Supplementary Material 6 and Supplementary

RESULTS AND DISCUSSION
Alfalfa and P. syringae proteomes were randomly paired to identify the genome-wide protein-protein interactions. The interaction probability of each pair was estimated individually through computational approaches: domainbased and interolog-based models ( Table 1). The interactome, combined from both the computational approaches, contained a total of 14,186,848 putative PPIs, involving 50,629 alfalfa and 2,932 P. syringae proteins (details of these HPIs is presented in Supplementary Material 2, excel sheets 3-17). To reduce the number of false-positive interactions we selected 690,634 non-redundant HPIs predicted by both the approaches, involving 10,935 alfalfa and 1,386 P. syringae proteins, named as common interactions throughout the manuscript (detail of these HPIs is presented in Supplementary Material 3, excel sheet 3). The visualization of the common network is shown in Figure 2. We present a thorough review of protein hubs and the functional enrichment of proteins found in the "common" subnetwork in the following sections.

Protein Hubs
Protein network analysis is widely implied to investigate the significant nodes in an extensive network (Ashtiani et al., 2018). Simultaneous interactions between PPI hubs and proteins have been observed in biological networks. Identifying the role of such interactions could further bolster the understanding of the infection mechanism of the pathogen (Yang et al., 2019).
In this study, we discovered a large number of protein hubs. Protein hubs are an important component of a HPI network, providing deep insights into various biological pathways, molecular processes, cellular biochemistry, and physiology (Kuzmanov and Emili, 2013), which is useful in deciphering the virulence mechanism of a host-pathogen system. The detailed analysis information on protein hubs is provided in Supplementary Material 4, excel sheets 1, 2. Despite the fact that hubs were determined based on node degree, a greater collection of topological measures were calculated to provide a more detailed portrait of the information. The following is a brief description of the values of certain topological measures within the common network.

Degree
The metric used in this study to evaluate hubs and nodes was the degree. The P. syringae proteins had an average degree of 498, while the alfalfa proteins had an average degree of 63. Pathogen proteins, as predicted, have a higher degree (and a larger betweenness) than other proteins in the popular subnetwork. The ratio of the interacting proteins in common network is coherent with previous reports of computational PPI prediction, whereby a small number of pathogen proteins interact with the host interactome (Li Z. G. et al., 2012;Kurubanjerdjit et al., 2013). Studies in the past have established that during the process of infection in the host, a pathogen extensively mutates itself. While on the other side, the plant expands its gene families in response to the pathogen attack (Stahl and Bishop, 2000). Thus, the reason for the ratio of interacting proteins being maintained throughout. We reported a large number of protein hubs in this study; the top 20 hubs from P. syringae as well as the alfalfa are discussed below. Summary of the 20 hub nodes for P. syringae and alfalfa proteins, as well as the number of their interacting partners are presented in Figure 3.

Closeness Centrality
The closeness centrality of a node in a network shows how close it is to all other nodes in the network. The average of the shortest path lengths from the node to every other node in the network is used to measure it. The average closeness centrality value for alfalfa and P. syringae proteins in the common network was 0.333 and 0.332, respectively. The highest closeness centrality value of 1 was found between five alfalfa proteins (MSAD_214560, MSAD_257396, MSAD_258406, MSAD_298516, MSAD_310524) and 11 P. syringae proteins (Psyr_2612254914, Psyr_2612251287, Psyr_2612254775, Psyr_2612254957, Psyr_2612251834, Psyr_2612252057, Psyr_ 2612254066, Psyr_2612251255, Psyr_2612253067, Psyr_26 12254813, Psyr_2612255124).

Betweenness Centrality
The betweenness centrality measures how close a node is to other nodes. The number of shortest paths that pass through the target node is used to calculate this metric. The average betweenness centrality value for alfalfa and P. syringae proteins in the common network was 8.34e-3 and 1.49e-2, respectively. The highest betweenness centrality value of 1 was found between 1 alfalfa protein (MSAD_321294) and 1 P. syringae protein (Psyr_2612250718).

Medicago sativa Hubs
The protein hubs network revealed that the strongly interacting host protein (MSAD_216261) was discovered to be cyclic nucleotide-binding domain (CNBD) protein. CNBD proteins are mostly involved in pathogen effector (avirulence proteins) recognition, which further leads to downstream signaling activation, and hence pathogen resistance (DeYoung and Innes, 2006). A study on Arabidopsis experimentally identified 12 cyclic nucleotide binding proteins, of which eight proteins contain the cyclic nucleotide binding domain, and these proteins are involved in hydrogen peroxide signaling and immune response (Donaldson et al., 2016). The protein, MSAD_321294, was found to be sensitive to aluminum rhizotoxicity. Aluminum toxicity has been shown to reduce root growth, and alteration in abscisic acid levels in root apices indicate its role in the plant response to aluminum toxicity (Kopittke, 2016). MSAD_220741, a chaperone protein encoded by the DnaK gene, was present in the top 20 host hubs. This gene is involved in plant response to various stresses. In rice seedlings, DnaK was shown to be upregulated in response to heavy metal stress (Ul Haq et al., 2019). DnaK regulates hydrogen peroxide production and ABA-induced antioxidant response during heat and drought stress (Yu et al., 2015). The majority of the top 20 alfalfa hub proteins were found to be heat shock proteins, classified into heat shock cognate 70 kDa proteins (MSAD_230624, MSAD_317196, MSAD_230523, MSAD_230540, MSAD_245934, MSAD_ 256658, MSAD_219293), and heat shock 70 kDa protein (MSAD_289443, MSAD_310640, MSAD_290843, MSAD_218590, MSAD_274579, MSAD_308554, MSAD_ 315555, MSAD_328469, MSAD_218587). Heat shock cognate 70 (Hsc) proteins are involved in protein folding and intracellular targeting, while the heat shock 70 (Hsp70) proteins are involved in signal transduction and protein translocation during abiotic and biotic stresses (Kominek et al., 2013). Hsp70s play a role in both infection mechanisms and the host response to stress. Studies show that Hsp70 is the primary target for P. syringae virulence effector, HopI1. HopI1 hijacks plant Hsp70 and localizes itself to the chloroplast (Jelenska et al., 2010).

P. syringae Hubs
The superfamily II DNA and RNA helicases, which include Pseudomonas proteins Psyr_2612251937 and Psyr_2612254566, form the largest pathogen hub in the common interaction network. Superfamily II helicases are divided into DEAD and DExH box families based on the sequence of a conserved motif (Steimer and Klostermeier, 2012). These helicases play an essential role in bacterial replication. Studies in E. coli reveal that absence of DEAD-box helicase leads to growth defects under laboratory conditions (Jagessar and Jain, 2010). The pathogen protein "Psyr_2612250718" is a member of the Major Facilitator Superfamily (MFS) transporter family, which is ubiquitously found in prokaryotes. A MFS transporter (mfsG) was discovered in Botrytis cinerea that was involved in isothiocyanate (ITC) detoxification and showed upregulation on interaction with wild-type Arabidopsis in planta (Vela-Corcía et al., 2019). Five proteins (Psyr_2612255229, Psyr_2612253263, Psyr_2612252226, Psyr_2612252166, Psyr_2612252474) in the common network were identified as 3-oxoacyl-(acyl-carrierprotein) reductase (OAR). OAR is responsible for catalyzing 3-oxoacyl-ACP reduction in the fatty acid synthesis pathway. In Pseudomonas aeruginosa PAO1, 12 OAR-encoding genes have been identified, which play a key role in the production of specific quorum-sensing signals in the bacteria (Guo et al., 2019). The pathogen proteins (Psyr_2612253558, Psyr_2612250706, Psyr_2612253590, Psyr_2612252788, Psyr_2612252142) belonged to the NAD(P)-dependent dehydrogenase family.
A novel NADP+-dependent D-arabitol dehydrogenase enzyme, identified from rust fungus Uromyces fabae, showed increased concentration during pathogenesis and quenched reactive oxygen species in the host immune response (Link et al., 2005). Another chaperone protein Psyr_2612252608, encoded by HtpG, was found in the common network. High-temperature protein G (HtpG), a Hsp90 homolog, is reported to be involved in the cell protection against environmental stress (Grudniak et al., 2015).

Functional Enrichment Analysis
Functional enrichment analysis is a critical step in determining the biological importance of proteins involved in host and pathogen PPIs. We examined the functional compositions of the respective proteins by analyzing GO and KEGG pathways enrichment. The existence of enriched (over-represented) functional categories were found to be closely linked to host defense and pathogen infection validates predicted HPIs.

Gene Ontology
To determine the significant functions of the alfalfa and P. syringae proteins involved in PPIs, the proteins were investigated using GO enrichment analysis based on three functional categories of gene ontology viz. cellular part, molecular function, and biological process. The enrichment was performed using enrichment score [−log10(P-value)] of the GO terms.
A total of 573 GO terms were found to be enriched for P. syringae proteins in the interactions. The detailed enrichment analysis data is present in Supplementary Material 4, excel sheet 4. The most enriched GO terms in the cellular component category were GO:0043229 (intracellular organelle), GO:0043232 (intracellular non-membrane-bounded organelle), and GO:0005840 (ribosome). In the molecular function category, GO:0000166 (nucleotide binding), GO:1901265 (nucleoside phosphate binding), and GO:0016887 (ATPase activity) were found to be most enriched. While in the biological process category, GO:0006412 (translation), GO:0055114 (oxidationreduction process), and GO:0019538 (protein metabolic process) were found to be abundant. Top 15 most enriched GO terms from all three categories are depicted in Figure 4. Aside from these, a number of GO terms (GO:0071944, GO:0030313, GO:0031975, GO:0009279) were linked to the pathogen's external membrane and the ice-nucleation gene (Psyr1608). This gene has also been found in Pseudomonas syringae pv. syringae B728a, which aids in ice crystal formation (Feil et al., 2005). P. syringae protein enrichment revealed the importance of the host for pathogens to carry out different molecular and biological processes.
For alfalfa proteins, 2127 GO terms were enriched. The detailed enrichment analysis data is present in Supplementary Material 4, excel sheet 3. The most enriched GO terms in biological process were GO:0070887 (cellular response to chemical stimulus), GO:0006979 (response to oxidative stress), and GO:0032870 (cellular response to hormone stimulus). The highly enriched GO terms for cellular component were GO:0005794 (Golgi apparatus), GO:0005618 (cell wall), and GO:0030312 (external encapsulating structure). In molecular function category, the most enriched terms were GO:0016614 (Oxidoreductase activity), GO:0015291 (secondary active transmembrane transporter activity), and GO:0015297 (antiporter activity). For each category top 15 enriched GO terms are depicted in Figure 5.
Chloroplast stroma (GO:0009570) was also found to be enriched in alfalfa proteins. A total of 3,382 interactions were discovered whereby 94 alfalfa proteins were found to be interacting with 754 pathogen proteins (Figure 6). It has previously been established that homologs of many components of bacterial signal recognition particle (SRP) pathways are present in the chloroplast. These SRP systems assist in the transport of proteins to the thylakoid membrane, which serves as a bacterial infection site (Ziehe et al., 2017). Chloroplast also serves as a production house of various defense-related signals [reactive oxygen species (ROS)] and hormones [abscisic acid (ABA), jasmonic acid (JA)], which play a crucial role in the immune response of plants. These are a primary site for the pathogen effectors to attack and overcome plant immune signaling (Park et al., 2018). ABA synthesis occurs in the chloroplast and is an essential plant growth regulator that induces stomatal closure in response to pathogen attack, thereby restricting pathogen entry (Lu and Yao, 2018). ROS contributes to the induction of hypersensitive response (HR) and cell wall strengthening (Kretschmer et al., 2020). Reduced ROS production in the chloroplast caused the reduction in HR in tobacco on the attack with Xanthomonas campestris pv. vesicatoria (Zurbriggen et al., 2009).
Protein phosphorylation (GO:0006468) was found to be involved in about 39 alfalfa proteins. These proteins interact with 395 pathogen proteins, owing to 2,701 PPIs. Phosphorylation  is considered as one of the most common post-translational modifications, and about one-third of all eukaryotic proteins are believed to be phosphorylated (Olsen et al., 2006), which usually occurs on threonine (Thr) and serine (Ser) residues (de la Fuente van Bentem and Hirt, 2009). In past, phosphorylation has been demonstrated to play a role in the generation of immune responses in plants (Park et al., 2012). In Arabidopsis, the treatment of cell cultures with flg22 (flagellin peptide) or xylanase (fungal elicitor) lead to the identification of 1,170 novel phosphorylation sites from 472 phosphoproteins, which evoke an early immune response (Benschop et al., 2007), thus indicating the differential phosphorylation of proteins and the importance of phosphorylation in plant immunity. Additionally, the interaction networks generated for top biological process (GO:0042744; hydrogen peroxide catabolic process), and molecular function (GO:0015291; secondary active transmembrane transporter activity) are depicted in Supplementary Figures 6, 7, respectively.

Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment
Establishing a connection between the role of important genes (and/or proteins) and specific pathway provides more information about the importance of a gene in various mechanisms. KEGG pathway enrichment analysis was performed to better understand the biological aspects of the proteins in the predicted HPI network. We found 116 KEGG pathways enriched for alfalfa proteins versus 86 pathways enriched for P. syringae proteins in total. The detailed pathway enrichment analysis data is present in Supplementary Material 4, excel sheets 5, 6 for alfalfa and P. syringae, respectively. Top 20 enriched KEGG pathways for alfalfa and P. syringae proteins are depicted in Figures 7, 8, respectively.  immunity (PTI) being the primary response, while effectortriggered immunity (ETI) is known as the secondary. In Petunia, the plant-pathogen interaction pathway was found to be enriched in both pollinated and unpollinated corollas, and the defenserelated genes (mostly interacting with Ca 2+ ) in this pathway showed upregulation during flower senescence (Broderick et al., 2014). From the network (Figure 9), it was analyzed that the host protein, MSAD_331284 (with the highest degree), interacts with 396 proteins. Besides plant-pathogen interaction (mtr04075), this protein is also involved in a variety of plant defense-related pathways, including the MAPK signaling pathway (mtr04016) and plant hormone signal transduction. The protein was also found to interact with one of the ice-nucleation proteins (INPs), which serve as the nuclei for ice formation, and hence initiates the disease.
ABC transporters pathway (mtr02010) was found to be associated with 54 alfalfa proteins. ABC transporters play a significant role in plant-pathogen interactions (Kang et al., 2011). Pleiotropic drug resistance (PDR), a subfamily of ABC transporters, is associated with plant defense mechanisms. A study on Arabidopsis discovered a gene (AtPDR12) that putatively encodes PDR, and the expression of this gene against the pathogen requires sensitivity to jasmonates and salicylic accumulation, both of which play an essential role in plant immunity (Campbell et al., 2003). In another study, sclareolide (an antifungal diterpene) was used to treat the cell cultures of Nicotiana plumbaginifolia, which resulted in identification of NpABC1 gene that encodes ABC transporter and is also related to the secretion of secondary metabolites associated with plant defense. NpABC1-encoded protein was discovered to be localized in the plasma membrane (Jasiński et al., 2001).
In the MAPK signaling pathway (mtr04016), we found 102 alfalfa proteins involved in 5,566 interactions with 597 pathogen proteins (Figure 10). Mitogen-activated protein kinases (MAPKs) are implicated in signal transduction in response to plant stress and development (Bigeard and Hirt, 2018). It was reported that BWMK1, a rice MAPK, showed rapid activation of protein kinase activity in presence of different pathogen signals such as salicylic acid, jasmonic acid, ethylene, and hydrogen peroxide (Cheong et al., 2003). Plant defense signaling was further demonstrated by the activation of two MAPKs, SIPK and WIPK, in Nicotiana tabacum in response to pathogen-related signals and various abiotic stresses, thus indicating the role of MAPK in plant stress response (Seo et al., 1995). The activation of the orthologs of SIPK and WIPK in alfalfa and Arabidopsis has been observed under stress conditions (Zhang and Klessig, 1997;Nühse et al., 2000).
A total of 305 proteins in alfalfa are linked to protein processing in endoplasmic reticulum (mtr04141). The endoplasmic reticulum (ER) is considered as the production site for secretory proteins (Yamada et al., 2011), and regulates hormone biosynthesis (Eichmann and Schäfer, 2012). The interference of plant physiological conditions with protein folding or accumulation of misfolded proteins in the ER results in ER stress. ER quality control system regulates such conditions and generates the unfolded protein response (UPR), leading to the elimination of misfolded proteins from the secretory pathway or upregulation of components required for protein folding. UPR in plants is mediated by specific membrane-associated transcription factors such as the bZIP family (bZIP17 and bZIP28) (Srivastava et al., 2012). In Arabidopsis, a few transmembrane bZIP transcription factors (AtbZIP17, AtbZIP28, and AtbZIP60), residing in the ER have been reported (Urade, 2009;Liu and Howell, 2010), which are engaged in response to ER stress. Beta-Alanine metabolism (mtr00410) involves 49 alfalfa proteins, which has a key role in various plant metabolisms, and also acts as a defense compound against various physiological stresses in plants (Parthasarathy et al., 2019). Transcriptomic studies in Arabidopsis have identified a considerable number of genes involved in amino acid metabolism against avirulent Pseudomonas syringae pv. tomato (AvrRpt2) (Scheideler et al., 2002). Despite of reduced glutamine, alanine, and proline levels than wild type plants, Arabidopsis mutant lysine histidine transporter 1 (lht1) showed enhanced resistance to various fungal, oomycete, and bacterial pathogens (Rojas et al., 2014). The majority of alfalfa proteins are involved in plant hormone signal transduction (mtr04075), thus controlling the plant growth and development in response to biotic and abiotic stresses (Schepetilnikov and Ryabova, 2017). Several metabolic pathways were also found to be linked to the predicted proteins in HPIs, including pyruvate metabolism (mtr00620), glycolysis/gluconeogenesis (mtr00010), glycerophospholipid metabolism (mtr00564), and fatty acid metabolism (mtr00564) (mtr01212).
On the other hand, P. syringae proteins are involved in secondary metabolite biosynthesis (psb01110). Secondary metabolites in pathogens have been shown to mimic plant effector molecules (such as auxins and abscisic acid) and subvert the plant defense mechanism (Pusztahelyi et al., 2015). Furthermore, pathogen proteins are found to be enriched in ABC transporters pathway (psb02010). This class of protein families is suggested to be important for phytopathogenesis (Zeng and Charkowski, 2020).
The predicted HPIs using computational methods contain possible candidates for experimental validation of the HPIs, which further provide insights into the infection mechanisms of P. syringae, according to the functional enrichment.

Subcellular Localization
A pathogen directs effector proteins to the cytoplasm of host cells in order to inhibit host immunity. These effector proteins are carried to various subcellular locations after crossing the host plasma membrane and destabilize the host immune system, allowing pathogen development (Zalguizuri et al., 2018). Predicting the subcellular localization of P. syringae-targeted alfalfa proteins may thus provide insight into HPI mechanisms. The presence of targeted proteins in cellular components, a possible site for pathogen infection in an organism, supports this theory. We predicted the subcellular localization of alfalfa and P. syringae proteins to gain a better understanding of the position of host-pathogen interactions in the cell. The distribution of the subcellular localization of the predicted proteins in alfalfa and P. syringae is depicted in Figure 11, and the detailed list is available in Supplementary Material 5, excel sheets 1, 2. The subcellular localization of the proteins involved in the common interactions revealed that 23.22% of alfalfa proteins are localized in the cytoplasm, 18.23% are membrane-associated. While in Pseudomonas, 65.37% of proteins are localized in the cytoplasm and 19.62% in the plasma membrane. The results indicate that most of the interactions between host and pathogen proteins occur in the cytoplasmic region and membranes.
We were also curious to see where the P. syringae proteins were found after the pathogen had multiplied in the alfalfa cells. As a result, we determined the distribution of alfalfa protein localization categories for each interacting pathogen localization class. For example, 906 cytoplasmic proteins from P. syringae interacted with 10,004 alfalfa proteins in 491,034 interactions, out of which 24.23% (2,424 of 10,004 in 71,469 interactions) are localized in cytoplasm followed by other membrane-associated 18.03% (1,804 of 10,004 in 130,352), nucleus 11.55% (1,156 of 10,004 in 38,728 interactions) multi-target 9.61% (962 of 10,004, in 43,944 interactions), and 8.05% (806 of 10,004 in 74,576 interactions) as cell membrane (Supplementary File 5 excel sheets 1, 2). In the multi-target (66 out 1,386 in 42,395 interactions) category, the majority of alfalfa proteins were found to be targeted toward other membrane, cytoplasm, cell FIGURE 11 | Distribution of subcellular localization of host proteins (blue) and pathogen proteins (red) involved in the predicted interactions (consensus from interolog and domain-based methods). membrane and nucleus category. This shows that the cytoplasmic pathogen proteins can interact inside the cytoplasm or in other membrane or nucleus of the host cell. Similarly, cytoplasmic membrane proteins (272) from P. syringae interacted with alfalfa 7,385 proteins in 134,597 interactions, which are localized in 23.39% (1,728/7,385 in 53,528 interactions) other membraneassociated, followed by cytoplasm 18.06% (1,334/7,385 in 11,900 interactions), cell membrane 11.11% (821/7,385 in 29,529 interactions), nucleus 9.24% (683/7,385 in 12,138 interactions), and 8.66% (640/7,385 in 9,321 interactions) as multi-target, etc. The 38 P. syringae periplasmic proteins interacted with 2,807 alfalfa proteins in 6,766 interactions, out of which 25.50% (716/2,807 in 1,482 interactions) were localized in other membrane-associated followed by cytoplasm 14.03% (394/2,807 in 1,590 interaction), and nucleus 11.72% (329/2,807 in 672 interactions), etc. The seven P. syringae proteins interacted with 1,965 alfalfa protein in 2,550 interactions, out of which 23.05% (453/1,965 in 592 interactions) were localized in other membrane-associated followed by cytoplasm 15.01% (295/1,965 in 346 interaction), and chloroplast 11.85% (233/1,965 in 406 interactions), etc. From these results, we can infer that P. syringae proteins mostly interacted with other membrane-associated alfalfa proteins followed by cytoplasm, cell membrane, nucleus, multi-target chloroplast.

Identification of Novel Host Targets for Known and Potential Effectors of P. syringae
The novel host targets for known effectors were identified by searching the known 26 effectors in the consensus predictions, and it was found that four effectors interact with 1,787 host proteins, involved in 1,787 interactions. More than half of the identified host proteins were found to be the members of ABC transporter and kinase families, including LRR-receptor-like kinase, calmodulin-binding receptor-like cytoplasmic kinase, MAP kinase, calcium-dependent kinase, calcium-dependent protein kinase (CDPK), serine/threonine kinase, and others. DEAD-box RNA helicase was annotated in 56 host proteins. In rice, a gene (OsBIRH1) encoding a DEAD-box RNA helicase protein was cloned, which showed activated expression on treatment with defense-related signals in seedling leaves. Also, the OsBIRH1-overexpressing transgenic Arabidopsis plants were tolerant to oxidative stress (Li et al., 2008).
To identify the novel host targets for potential effectors, 646 potential effectors were searched in the consensus predictions, which identified 190 effectors with 8,266 host interacting partners, involved in 94,243 interactions. The host proteins were classified as disease resistance protein, DEAD-box RNA helicase, LRR and NB-ARC domain disease resistance, kinases, glycoside hydrolase family 1, and ubiquitin-conjugating enzyme. In higher plants, glycoside hydrolase 1 (GH1) is known to play an important role in plant defense against various biotic and abiotic stresses, as well as other processes such as lignification, cell wall hydrolysis (Xu et al., 2004). Ubiquitination has a great importance in response to plant defense against pathogens. A study in wheat identified a ubiquitin-conjugating enzyme (TaU4), which is localized in nucleus and cytoplasm, and acts as a negative regulator against Zymoseptoria tritici (Millyard et al., 2016).

Identification of Novel Virulence Effectors
Additionally, we searched for the novel P. syringae effectors in the predicted HPIs. For this, we analyzed the P. syringae protein dataset on EffectiveDB at default parameters. From EffectiveDB, we identified around 3,927 P. syringae proteins which were classified as EffectiveT3, EffectiveCCBD, and EffectiveELD. These proteins were then identified in our predicted HPIs and it was found that these effectors count for 686,184 interactions. We obtained 1,371 novel virulence effectors interacting with 10,934 host proteins. Each pathogen protein was found to interact with multiple host proteins.

Interolog and Domain Combined Interactions
The interactome, combined from both the computational approaches, contained a total of 14,186,848 putative PPIs, involving 50,629 alfalfa and 2,932 P. syringae proteins. The domain-based approach predicted 2,960,634 interactions, involving of 14,328 alfalfa and 2,515 P. syringae proteins, while the interolog-based approach predicted 11,916,848 interactions, involving of 49,995 alfalfa and 2,202 P. syringae proteins. Supplementary Material 2, excel sheets 1, 2, contains a catalog of both host and pathogen proteins, as well as their sequences. Due to the large number of interactions and computational limitations, we were unable to generate/visualize the proteinprotein interaction network. Functional enrichment analysis revealed that enriched GO terms are consistent with enriched terms from the "common" interactions as discussed above.

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
RKau formulated and designed the research and contributed to supervision, project administration, and funding acquisition. RKat analyzed the data, developed prediction models, and performed functional enrichment analysis, validations, literature mining, and contributed to writing-original draft preparation. RKau and ND contributed to writing-review and editing. RKat, ND, and RKau contributed to visualization. All authors have read and agreed to the published version of the manuscript.

FUNDING
The authors acknowledge the support to this study from Utah Agricultural Experiment Station (UAES), Utah State University, seed grant # A48025 to RKau. This research was also partially supported by UAES office and approved as journal manuscript number 9553. The funding body did not play any role in the design of this study; the collection, analysis, or interpretation of data; or in the writing of this manuscript.

ACKNOWLEDGMENTS
We are thankful to the members of KAABiL bioinformatics team for their valuable inputs throughout the research analysis. We sincerely thank the referees for all the suggestions and help in improving the research article.