Protein Markers for the Identification of Cork Oak Plants Infected with Phytophthora cinnamomi by Applying an (α, β)-k-Feature Set Approach

Abstract

Cork oak decline in Mediterranean forests is a complex phenomenon, observed with remarkable frequency in the southern part of the Iberian Peninsula, causing the weakening and death of these woody plants. The defoliation of the canopy, the presence of dry peripheral branches, and exudations on the trunk are visible symptoms used for the prognosis of decline, complemented by the presence of Phytophthora cinnamomi identified in the rhizosphere of the trees and adjacent soils. Recently, a large proteomic dataset obtained from the leaves of cork oak plants inoculated and non-inoculated with P. cinnamomi has become available. We explored it to search for an optimal set of proteins, markers of the biological pattern of interaction with the oomycete. Thus, using published data from the cork oak leaf proteome, we mathematically modelled the problem as an α, β-k-Feature Set Problem to select molecular markers. A set of proteins (features) that represent dominant effects on the host metabolism resulting from pathogen action on roots was found. These results contribute to an early diagnosis of biochemical changes occurring in cork oak associated with P. cinnamomi infection. We hypothesize that these markers may be decisive in identifying trees that go into decline due to interactions with the pathogen, assisting the management of cork oak forest ecosystems.

1. Introduction

Cork oak decline in Mediterranean forests is a complex phenomenon, observed with remarkable frequency in the southern part of the Iberian Peninsula, causing the weakening and death of these woody plants [1,2,3,4]. Mediterranean oak forest dynamics mirror the impact of anthropogenic factors and environmental variables over time, with the expectation of a more pronounced decline in the long term [5,6,7]. Canopy defoliation, dry peripheral branches, black spots, and exudations on the trunk are visually observable symptoms used for decline prognostics, sometimes complemented by the identification of Phytophthora cinnamomi in the rhizosphere of trees and the adjacent soils [8,9]. This oomycete invades the roots of the cork oak leading to cell death of the infected tissues, triggering a localized and distal defence response [10,11,12]. In the aerial part of the plant, the evidence of the plant’s interaction with the parasite depends on the extent and degree of the necrotized roots [10]. With global warming, cork oaks are exposed to extreme climate variations that do not favour resilience to biotic stresses such as the interaction with P. cinnamomi [13]. The resilience patterns observed in natural populations are diversified [14], in accordance with the genetic variability typical of this long-lived plant, which constitutes an asset in a long-term process of adaptation [15,16]. Based on the sixth National Forest Inventory data, only 5% of the cork oak stands in Portugal do not show any damage such as defoliation and/or discoloration of the canopy [17]. According to these data, trees are moving away from an optimal vegetative state and the early physiological/biochemical assessment may be decisive in order to understand the tree’s level of resilience to overcome adversity and to recover. Recently, a large proteomics dataset obtained from the leaves of cork oak plants inoculated and not inoculated with P. cinnamomi was made available [12] and can be explored in order to search for an optimal set of proteins, markers of the biological decline pattern. This work has yet to be done.

Statistical methods are commonly used to reduce dimensions and classify data, e.g., principal component analysis (PCA) and exploratory factor analysis (EFA). They are applied to extract features and knowledge from data but require high-dimension data or correlation and distribution conditions.

Combinatorial optimization approaches generally resort to model the research question as a covering problem, either based on features or in patterns, and have been widely applied to solve real world problems in physical and biological sciences such as molecular biology and biochemistry, as well as in engineering and computer sciences [18]. Among the cover problems, we highlight the (α, β)-k-Feature Set Problem (α, β-k-FSP), which can be reduced to the Red-Blue Bipartite Dominating Set problem and is a generalization of the k-Feature Set Problem (FSP). It is also used to extract features and knowledge from data as well as to reduce dimensions and classify data, having been mostly applied in medical sciences helping with the explanation and diagnosis of disease. For example, it was applied to some cancer types [19,20,21,22], Alzheimer’s disease [23,24,25], and epilepsy [26]. As far as we know, this is the first time that it is applied to plant biology, namely to the cork oak–P. cinnamomi interaction. The k-FSP attempts to reduce a given data set by selecting a relevant subset of k-features (minimizing k) that contributes to a better explanation and understanding of the whole data, eliminating unrelated, redundant, or conflicting features. The (α, β)-k-Feature Set Problem (α, β-k-FSP), proposed by Cotta et al. [27], aims to determine the minimum subset of features that enlightens the dichotomy within the samples while reducing the data size. It intends to minimize the number of features (k) and simultaneously to maximize the differences between samples of different groups (α) and the similarities between samples of the same group (β). The selected feature set is as reliable as the greater α and β, although a larger feature set will be obtained. The α, β-k-FSP is NP-hard, since it is a generalisation of a problem that is NP-complete, α = 1 and β = 0 k-FSP. The latter can be reduced to the k-Vertex Cover Problem [28], known to be NP-complete [29]. Therefore, there is no known algorithm to solve exactly the α, β-k-FSP in polynomial time, meaning that it is not probable that there is an algorithm able to find the optimum solution for any instance of the problem in a computational time proportional to a polynomial function of the size of the input, being necessary to use approximate methods. Nevertheless, for small and even for some medium size instances, α, β-k-FSP may be solved exactly, in a reasonable time. Using the published cork oak protein dataset, we mathematically model the problem as an α, β-k-FSP and analyzed the obtained results. We select leaf proteomic markers, finding a set of proteins (features), as small as possible, that represent the dominant effects on the metabolism of the host resulting from the action of the oomycete on the roots, in order to contribute to an early diagnosis of the biochemical changes associated with cork oak infection by P. cinnamomi. Variations in the abundance of these proteins have an impact on the processes of protein production, essential for cellular reprogramming in situations of biotic stress, beta-oxidation of fatty acids in peroxisomes, photosynthesis, and photorespiration. We hypothesize that these markers can be decisive to detect trees that go into decline due to the interaction with the pathogen, helping to manage forest ecosystems in which the cork oak predominates.

2. Materials and Methods

2.1. Public Cork Oak Protein Dataset

In this study, we used 80 proteins out of 424 proteins quantified by SWATH-MS (sequential windowed data-independent acquisition of the total high-resolution mass spectra) in cork oak leaves inoculated in the roots with P. cinnamomi [12]. The biological assay was carried out with 12 cork oak plants, with six being inoculated with P. cinnamomi. Proteins included in the selected group of 80 were identified using the Arabidopsis proteome database as reference showed a p-value below 0.05 or a fold change greater than 2 or less than 0.5 in their ratio levels between inoculated and non-inoculated (control) samples, then being considered as altered. The aforementioned mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE [30] partner repository with the dataset identifier PXD021455 and is publicly available.

2.2. Modeling the Cork Oak Protein Dataset as an (α, β)-k-Feature Set Problem

As mentioned before, in the introduction section, we aim to determine a dominant subset of proteins (features) that robustly explains the dichotomy within the cork oak samples. Next, we will present the general definition of the α, β-k-FSP [20,27] followed by its mathematical formalization for our case and a small example.

Definition 1.

Consider a given data set with m samples and n features, represented by a discrete matrix D and an array T. Each element d_if of matrix D specifies the value of feature f in sample i. Each element t_i of array T specifies the type or class of each sample. α, β, and k are positive integer parameters. Then the (α, β)-k-Feature Set Problem seeks to minimize the size k of a subset S of features (columns of matrix D), such that:

(a) for all pairs of samples (i, j), with i ≠ j, if they do not belong to the same type, t_i ≠ t_j, there is a subset S′ ⊆ S such that |S′| ≥ α and they have different values, d_if ≠ d_jf, for all feature f of S′;

(b) for all pairs of samples (i, j), with i ≠ j, if they belong to the same type, t_i = t_j, there is a subset S″ ⊆ S such that |S″| ≥ β and they have the same value, d_if = d_jf, for all feature f of S″.

So, the α, β-k-FSP seeks to determine a subset of features, with minimum cardinality, having at least α features for differentiating between any two samples of different classes and at least β features for coupling any two similar samples of the same class.

In our study, considering our experimental data set, we want to determine the minimum subset of proteins that consistently explains the differences between any pair of control and inoculated samples and the similarities between any pair of control samples or of inoculated samples. This can be modelled as an α, β-k-FSP, as described below.

We used the following notation for the experimental data set:

$C=\left\{c_1, c_2,\dots ,c_m \right\}$ is the set of cork oak samples;

$T=\left\{t_1,t_2,\dots ,t_m\right\}$, $t_i\in \left\{\mathrm{control},\mathrm{inoculated}\right\},\forall i=1,\dots ,m$, holds the type of the samples;

$P=\left\{p_1, p_2,\dots ,p_n \right\}$ is the complete set of biological or statistical relevant proteins obtained in the biological experiment;

$Data=\left(d{t}_{c_i p_l}\right)\in {ℝ}^{m\times n}$ is the matrix with the experimental data set holding the quantification level of each protein $p_l$ in each cork oak sample $c_i$;

$D=\left(d_{c_i p_l}\right)\in {\left\{0, 1\right\}}^{m\times n}$ is the binary matrix got from the data matrix classifying each protein quantification level as less abundant or more abundant represented, respectively, by 0 and 1.

Table 1. Data for a small example.

D		Proteins (P)
		${\mathit{p}}_1$	${\mathit{p}}_2$	${\mathit{p}}_3$	${\mathit{p}}_4$	${\mathit{p}}_5$	T
Samples (C)	$c_1$	0	0	1	1	0	Control
	$c_2$	0	1	1	0	1	Control
	$c_3$	1	0	1	0	1	Inoculated
	$c_4$	1	1	0	1	1	Inoculated

shows the above sets for a small example.

We also denoted the minimum proteins set obtained in the solution of the α, β-k-FSP as $S=\left\{p_{s1}, p_{s2},\dots ,p_{sk}\right\}\subseteq P$.

We defined matrices $A=\left(a_{\left(c_i, c_j\right) p_l}\right)$ and $B=\left(b_{\left(c_i, c_j\right) p_l}\right)$, respectively, for the pairs of samples that belong to different types and for the pairs of samples that belong to the same type as follows:

$$a_{\left(c_i,c_j\right)p_l}=\left\{\begin{array}{c}1 , \mathrm{if} d_{c_i{p}_l}\ne d_{c_j{p}_l} \mathrm{and} t_{c_i}\ne t_{c_j}\\ 0, \mathrm{otherwise}\end{array}\right.;$$

$$b_{\left(c_i,c_j\right)p_l}=\left\{\begin{array}{c}1 , \mathrm{if} d_{c_i{p}_l}=d_{c_j{p}_l} \mathrm{and} t_{c_i}=t_{c_j}\\ 0, \mathrm{otherwise}\end{array}\right.,$$

matrices A and B for the data of Table 1 are presented in

Table 2. Matrices A and B for the data of Table 1.

		Proteins (P)
		${\mathit{p}}_1$	${\mathit{p}}_2$	${\mathit{p}}_3$	${\mathit{p}}_4$	${\mathit{p}}_5$
Sample pairs	A
	$\left(c_1,c_3\right)$	1	0	0	1	1
	$\left(c_1,c_4\right)$	1	1	1	0	1
	$\left(c_2,c_3\right)$	1	1	0	0	0
	$\left(c_2,c_4\right)$	1	0	1	1	0
	B
	$\left(c_1,c_2\right)$	1	0	1	0	0
	$\left(c_3,c_4\right)$	1	0	0	0	1

.

Then, we formalized the α, β-k-FSP as a binary linear program as shown below, where the binary variable $x_{p_l}=1$ if the protein $p_l$ is selected, otherwise $x_{p_l}=0$.

$$\min k={{\displaystyle \sum }}_{l=1}^n{x}_{p_l}$$

(1)

subject to:

$${{\displaystyle \sum }}_{l=1}^n{a}_{\left(c_i,c_j\right)p_l}{x}_{p_l}\ge \alpha ,\forall \left(c_i,c_j\right)\in C:t_{c_i}\ne t_{c_j}$$

(2)

$${{\displaystyle \sum }}_{l=1}^n{b}_{\left(c_i,c_j\right)p_l}{x}_{p_l}\ge \beta ,\forall \left(c_i,c_j\right)\in C:t_{c_i}=t_{c_j}$$

(3)

$$x_{p_l}\in \left\{0,1\right\}$$

(4)

In the above mathematical formalization, the objective function (1) minimizes the number of selected proteins (which will belong to set S). Constraints (2) ensure that for all pairs of samples that belong to different types, there are in set S at least α proteins having different values and constraints (3) ensure that for all pairs of samples that belong to the same type, there are in set S at least β proteins having identical values. Constraints (4) force variables to be binary.

It is very useful to represent the α, β-k-FSP as a graph G(V, E), with $V=A\cup B\cup P$ nodes, and a set of edges, E, linking the nodes according to the 1 values of matrices A and B. So, we obtain the graph of Figure 1 for the example of Table 1 and Table 2.

Clearly, if we consider α = β = 1, then k = 1 is the optimum solution and $S=\left\{p_1\right\}$ because protein 1 dominates every node, since it is linked to all of them, see Figure 2a. This is due to the fact that it presents the same value to the same type pairs and opposite values to the different type pairs (see Table 1), being completely coherent. But as mentioned before the selected feature set, S, is as more reliable as the greater α and β are, so we can see in Figure 1 that the largest value that α and β can take is 2, because it is the number of edges linked to nodes $\left(c_2,c_3\right)$,$\left(c_1,c_2\right)$, and $\left(c_3,c_4\right)$, therefore, those nodes cannot be covered more than two times. For α = β = 2, it is clear that the proteins linked to those three nodes must be in the optimum solution, then $S^{\prime }=\left\{p_1,p_2\right\}$ covers node $\left(c_2,c_3\right)$ and $S^{″}=\left\{p_1, p_3,p_5 \right\}$ covers the other two nodes and no more proteins are needed since every node is covered twice (see Figure 2b), being k = 4 and $S=S^{\prime }\cup S^{″}=\left\{p_1,p_2, p_3,p_5 \right\}$ the optimum solution of α, β-k-FSP for this example. For the other values, α = 1 and β = 2 or α = 2 and β = 1, the solutions will be, respectively, k = 3, $S=\left\{p_1, p_3,p_5\right\}$ and k = 3, $S=\left\{p_1,p_2, p_4\right\}$. Note that the higher α and β, the higher the number of proteins selected (value of k). So, usually, for high-dimensional data sets where it is necessary to use approximate methods, that do not guarantee the optimum solution, the problem is solved for some different values of α and β and then the best compromise solution is chosen.

Based on the above formalization, we have implemented the procedures needed to prepare the data and solve the respective α, β-k-FSP in MATLAB (MATrix LABoratory, which is a multi-paradigm programming language). We briefly describe them below.

Input data: It is known, and has been confirmed in our experimental dataset, that statistical significance can be very different from biological meaningfulness. Thus, in this study, we considered the set of 80 proteins quantified for 12 samples, divided into two groups, control and inoculated, each with 6 samples, above-mentioned (in Section 2.1). So, using the above notation, the input was:

$C=\left\{c_1,c_2,\dots ,c_{12}\right\}$, $T=\left\{t_1,t_2,\dots ,t_{12}\right\}$, where $t_1,t_2,\dots ,t_6=\mathrm{control}\mathrm{and}{t}_7,t_8,\dots ,t_{12}=\mathrm{inoculated}$ and $P=\left\{p_1,p_2,\dots ,p_{80}\right\}$.

Construction of matrices: Our biological data set, matrix data, has continuous values for the protein level and it was necessary to consistently separate less abundant quantification level samples from more abundant quantification level samples. Thus, we established for each protein a range, centred on its mean, which clearly separates the less abundant samples (below the range) from the more abundant samples (above the range).

To determine which range should be used, we tested five ranges, being each equal to μ ± l × μ, where μ is the mean and l is equal to 5%, 7.5%, 10%, 12.5%, and 15%.

As matrix D must be binary or discrete, for each protein (column) we assigned 0 and 1, respectively, to less abundant quantification level and to more abundant quantification level samples, while the samples with a level belonging to the above aforementioned ranges were considered with neutral abundance and with no value assignment. Therefore, our matrix D is slightly sparse.

Matrices A and B were computed as defined in the above formalization. Matrix A has 36 lines (6 × 6 = 36 pairs of samples of different types) and matrix B has 30 lines (combination of 6 taken 2 at a time, for control pairs and also for inoculated pairs, that is, $2\times \frac{6!}{2! 4!}=30$ pairs of samples of the same type), both with 80 columns (number of proteins).

Parameters definition: As mentioned above, the largest possible value for α is equal to the minimum number of edges linked to the nodes representing pairs of samples of different types. Analogously, β is at most equal to the minimum number of edges linked to the nodes representing pairs of samples of the same type. Those were the α and β values computed and used in our model.

Optimum solution: The algorithm to solve our α, β-k-FSP was implemented in MATLAB. The binary linear program optimum solution was obtained using the MATLAB intlinprog function, which was appropriately parametrized. The algorithm took an average execution time of 0.4866 ms, varying between 0.374 ms and 0.599 ms, on an Intel i5-4300U CPU (Mobile Haswell, 15 W TDP, 2C/4T, 1.90 GHz nominal and max turbo of 2.90 GHz).

3. Results and Discussion

3.1. (α, β)-k-Feature Set Problem

The results obtained for the optimum solution of the α, β-k-FSP, considering the mentioned five ranges for neutral quantification, are shown in

Table 3. Optimum values obtained for the α, β-k-feature set problem.

	R1	R2	R3	R4	R5
	μ ± 0.05 μ	μ ± 0.075 μ	μ ± 0.1 μ	μ ± 0.125 μ	μ ± 0.15 μ
Bio/statistical significative proteins	80	80	80	80	80
Alpha (α)	16	11	8	6	5
Beta (β)	18	15	9	5	4
Optimum number of proteins (k)	51	47	29	20	16

.

It can be seen that for ranges R1 and R2, a high number of proteins are obtained (k), while for R5 the obtained protein set seems rather small. We also observed that the set of proteins obtained, for each range, is included (except for no more than four proteins) in the next larger range, as depicted in Figure 3, where two Venn diagrams were used for a better visualization. Therefore, and considering also the alpha and beta obtained values, we chose the R3 range as the best compromise between the number of proteins and the legitimacy in finding similarities (beta) and dissimilarities (alpha) between samples. Therefore, we will focus our discussion on that optimum set of proteins (S-R3).

The approach carried out allowed the selection of 29 proteins (R3 range), which ensure with a high certainty the existing separation between the group of control plants and the group of plants inoculated with P. cinnamomi.

Table 4. Description of the optimum set of proteins (S) for R3 range (S-R3).

S	Protein Information a
R3
Arabidopsis UniProt Accession	Protein Name	Initial
P16181	40S ribosomal protein S11-1	RS111
P42798	40S ribosomal protein S15a-1	R15A1
Q9STX5	Endoplasmin homolog	ENPL
P10795	Ribulose bisphosphate carboxylase small chain 1A	RBS1A/RBCS1A
Q9SIM4	60S ribosomal protein L14-1	RL141
Q9LF37	Chaperone protein	ClpB3
F4J3Q8	P-loop containing nucleoside triphosphate hydrolases superfamily	F4J3Q8
Q9SII0	Probable histone H2A variant 2	H2AV2
Q9LRR9	(S)-2-hydroxy-acid oxidase GLO1	GLO1/GOX1
Q9LXG1	40S ribosomal protein S9-1	RS91
Q9SVR0	60S ribosomal protein L13a-3	R13A3
P56761	Photosystem II D2	PSBD
F4JYM8	Thiolase family protein	F4JYM8/AACT1
O04486	Ras-related protein	RABA2a
Q9FZ47	ACT domain-containing protein ACR11	ACR11
P38418	Lipoxygenase 2	LOX 2
Q9FLN4	50S ribosomal protein L27	RK27
Q9FGX1	ATP-citrate synthase beta chain protein 2	ACLB2
Q9SRV5	5-methyltetrahydropteroyltriglutamate-homocysteine methyltransferase 2	METE2
P27140	Beta carbonic anhydrase 1	BCA1
Q9SCW1	Beta-galactosidase 1	BGAL1
B3H4S6	Dicarboxylate transporter 1	DiT1
O49485	D-3-phosphoglycerate dehydrogenase 1	SERA1
Q9LF98	Fructose-bisphosphate aldolase 8	ALFC8/FBA8
P27323	Heat shock protein 90-1	HSP90-1
F4KDZ4	Malate dehydrogenase	FAKDZ4/PMDH2
P56778	Photosystem II CP43 reaction center protein	PSBC
A0A1P8B485	Protein translocase subunit Sec A	AGY1
O81644	Villin-2	VILI2

^a UniProt accession codes, protein name, and protein initials arise from the annotation using the Arabidopsis proteome database as a reference.

details the information regarding the selected group of proteins, including the UniProt accession codes, name, and protein initial.

When graphically representing the protein abundance for S-R3 in the 12 cork oak samples evaluated, a distinct pattern of colours is observed for the two plant groups, control and inoculated (Figure 4). The method allowed for the identification of the smallest number of proteins out of the group of 80 that has the potential to effectively separate, based on the quantification of the proteins identified, plants that are in an interaction with P. cinnamomi from plants that are not in an interaction with the pathogen. For the first 16 proteins on the x-axis, the profile of lower abundance and neutral is characteristic of the control samples, the opposite occurring in the following 13 proteins. Proteins with consistent profiles of more abundant/neutral or less abundant/neutral quantification levels are more frequent in control samples (21 proteins) than in the inoculated samples (10 proteins), showing that the diversity of quantification profiles is more homogeneous in the control samples. This means that the immune response of cork oaks to infection by P. cinnamomi falls within a certain pattern but displays a higher diversity of protein profiles. By using the normalised abundance values for each protein in each of the samples, it can be seen that the quantification profiles of the control samples are plotted in opposite quadrants to those of the inoculated samples, showing the dissimilarity in production levels of each of the proteins in the two types of samples (Figure 5). The heterogeneity of abundance patterns observed for each protein in each of the samples is also evident (column height of the stacked columns in Figure 5). For example, protein P56778 (first column) is always more abundant in control samples than in infected samples but the abundance values vary from sample to sample. In the case of protein P10795 (column fifteen), there is no consistent pattern of abundance in the control samples or in the inoculated samples. In nature, a slow form of cork oak decline predominates, with trees showing gradual degrees of severity of the phenomenon. It is possible that the different patterns observed in the infected samples for the selected proteins have correspondence with the phenotypic characteristics of the trees with decline symptoms. In view of the high degree of genetic diversity observed in Q. suber, it is expected that the homeostatic state of the trees emerging from the interaction with P. cinnamomi is diverse and complex. Thus, the defence response of cork oaks to P. cinnamomi infection observed in this study is heterogeneous and this optimal set of proteins may represent the diversity of protein patterns expected to be found in a natural population challenged by this oomycete. The protein biological features may help in understanding the role they play in trees with a healthy or declining profile.

3.2. Biological Relevance of the S-R3 Protein Set

The S-R3 protein set was used as input data on STRING v11 online platform (search tool for retrieval of interacting genes/proteins) [31,32] to build a protein–protein interaction (PPI) network among the selected proteins and screen for the associated KEGG metabolic pathways [33]. The default parameters have been maintained, with the confidence cut off value set at 0.400 and false discovery rate (FDR) stringency set at 0.05. The PPI network is based on known interactions (from curated databases and experimentally determined), predicted interactions (gene neighbourhood, gene fusions, and gene co-occurrence) and information gathered from text mining, co-expression, and protein homology. The links between proteins exhibit the joint contribution to a given function and allow us to infer about their biological relevance as a group. STRING analysis revealed that 24 of the 29 proteins of S-R3 are connected with each other by 36 edges, with a PPI enrichment p-value < 0.02 (Figure 6). Thus, it can be seen that there are proteins in this group that are partly biologically associated with each other, functioning as a group.

In order to understand which groups of proteins within the network could correspond to functional units, a hierarchical clustering of the string network was carried out using the k-means clustering method. Three main clusters were identified: the blue cluster includes 7 proteins functionally associated with photosynthesis, the red cluster comprises 16 proteins associated mainly with glyoxylate and dicarboxylate metabolism, carbon fixation in photosynthetic organisms, and the citrate cycle, and the green cluster contains 6 proteins associated with ribosomes (Figure 7).

In each of the clusters, the proteins Q9FLN4-50S ribosomal protein L27 (RK27), PMDH2-Malate dehydrogenase (FAKDZ4/PMDH2), P10795-ribulose bisphosphate carboxylase small chain 1A (RBS1A/RBCS1A), and A0A1P8B485-protein translocase subunit Sec A (AGY1) stand out because they interact with many proteins and/or are the link between the group to which they belong and the other groups. The number of edges in the ribosomal protein cluster indicates their functional dependence on each other, that is, alterations produced in one protein should influence the others. Therefore, it is expected that variations in the abundance of the RK27 protein will condition the whole group of ribosomal proteins, affecting the structures and processes of protein production essential for cellular reprogramming in situations of biotic stress. Plastid translation is required for embryo development in Arabidopsis with the engagement of at least nine nuclear-encoded plastid ribosomal proteins (PRPS5, S13, S20, L1, L4, L6, L21, L27, and L35), but the exclusive absence of PRPL27 or L1, L4, or L35 leads to alterations in cell division patterns, showing that L27 is required for basal ribosome activity [34]. This reinforces the image of a pivot protein on the PPI network. Changes in cellular metabolic processes in which the PMDH2, P10795, and AGY1 proteins participate will influence the RK27 protein also, and thus all the proteins from the ribosome cluster.

AGY1 encodes a chloroplast subunit of the protein translocase Sec A and the name originates from albino or glassy yellow phenotypes of Arabidopsis mutants, dysfunctional in the Sec A protein-transporting ATPase [35]. The biosynthesis of chloroplasts is compromised in these homozygous mutants due to the inability to import nuclear-encoded proteins destined to be translocated to the lumen of the thylakoids by the Sec pathway [35,36]. This inability gives rise to genetic reprogramming aimed at the production of proteins related to photosynthetic complexes, protein translocation systems for chloroplasts, and mitochondrial respiratory complexes. In this context, the connection found between AGY1 and RK27 makes sense, since an alteration in AGY function triggers the reprogramming of gene expression at the level of various organelles and concomitant protein synthesis. In cork oak, the interaction with P. cinnamomi induced changes in the amount of AGY1 produced in inoculated plants, reducing the amount when compared with that produced in control plants [12].

The peroxisomal malate dehydrogenase PMDH2 catalyses the re-oxidation of NADH during fatty acid beta-oxidation [37] and this protein is concomitantly implicated in several metabolic pathways in the peroxisomes, such as the glyoxylate and dicarboxylate metabolism, photorespiration, nitrogen metabolism, and plant hormone biosynthesis [38]. Its positioning in the red cluster interacting with several proteins participating in the above biological processes is fitting. Proteins from the red cluster that had connections to PMDH2 and had associated enzyme numbers (EC numbers) were submitted to KEGG [33], choosing the Arabidopsis organism as a reference, to identify the metabolic pathways in which they participate (Figure 8). Four red cluster proteins are assigned to the glyoxylate and dicarboxylate metabolism KEGG map, with malate dehydrogenase (EC:1.1.1.37; PMDH2) associated with the methylaspartate cycle, the thiolase family protein (EC:2.3.1.9; AACT1) associated with the ethylmalonyl pathway of terpenoid biosynthesis, and (S)-2-hydroxy-acid oxidase GLO1 (EC:1.1.3.15; Q9LRR9) and ribulose bisphosphate carboxylase small chain 1A (EC:4.1.1.39; P10795) linked to photorespiration (Figure 8).

Moreover, using the data made available in a previous publication [12], it can be seen that all the proteins mentioned in the previous text were less abundant in cork oak samples inoculated with P. cinnamomi when compared with the controls, showing negative values for Log2FC (fold change ratio logarithm of protein abundance in inoculated over control samples). It is possible that the strategy adopted by the plant to reduce the energy to be expended in reorganising the physiological state resulting from the interaction with the oomycete involves decreasing the abundance of proteins, which interferes with several metabolic processes simultaneously, some being housed exclusively in the peroxisome. In the context of biotic stress, the NO signaling molecule could eventually be responsible for the regulation of photorespiration and beta-oxidation in peroxisomes through S-nitrosylation of the H₂O₂-producing GLO1 proteins (EC:1.1.3.15; Figure 8) or the PMDH isoforms, in agreement with the hypothesis raised by Ortega-Galisteo and colleagues [39].

3.3. Harmonisation between the PPI Network and the Coverage Problem Approach

Looking to the coverage problem approach, based on the quantification level of the proteins, a correspondence between the biological processes enrolled in the clusters of the PPI network and some proteins of the optimal set S-R3 stand out. In

Table 5. Biological processes that stand out from the S-R3 solution in the cover problem approach.

Biological Process	Arabidopsis UniProt Accession	Number (N) and Relative Frequency (f) of Abundance Quantifications
			Control Samples			Inoculated Samples
			Less	More	Neutral	Less	More	Neutral
Protein synthesis	P16181; P42798 Q9SIM4; Q9SVR0 Q9LXG1; Q9FLN4	N f	21 0.58	3 0.08	12 0.33	7 0.19	22 0.61	7 0.19
Photosynthesis	P56778; P56761	N f	0 0.0	6 0.5	6 0.5	7 0.58	0 0.0	5 0.42
Glyoxylate and dicarboxylate metabolism	PMDH2; AACT1; Q9LRR9; P10795	N f	6 0.25	17 0.71	1 0.04	14 0.58	6 0.25	4 0.17

, we considered protein synthesis, photosynthesis, and glyoxylate and dicarboxylate metabolism that stand out in the PPI network. We also present their quantification levels for six, two, and four proteins, respectively, for the six control and six inoculated samples. Each of these biological processes includes 36, 12, and 24 measurements of quantification levels for each type of samples, given that each biological process includes p proteins, thus, it has $N$ quantification level measures, for the control and for the inoculated samples, where $N=6\times p$.

Regarding the quantification abundance level of the six proteins, with a full match to the green cluster in the PPI network, we see that the control and inoculated samples present opposite behaviours. For the control samples, 91% of quantifications have less (58%) or neutral (33%) abundance, while 80% of quantifications for the inoculated samples have more (61%) or neutral (19%) abundance.

The group of proteins belonging to photosynthesis is also very clear; 100% of quantifications for the control samples show more (50%) or neutral (50%) abundance, while 100% of quantifications for the inoculated samples have less (58%) or neutral (42%) abundance. Proteins P56761 (EC:1.10.3.9) and P56778 from the blue cluster in the PPI network are components of photosystem II (PSII) that uses light energy to abstract electrons from water, generating molecular oxygen and a proton gradient subsequently used for ATP formation. The altered production of these proteins is a sign that the infection of Q. suber roots by P. cinnamomi interfere with the photosynthetic process that occurs in the leaves and in cells distant from the site of infection. Necrosis of the host roots caused by P. cinnamomi infection and the disappearance of the thinner roots limit water uptake by the plant, and the inefficient functioning of photosystem II may be related to water availability. With regard to plant–pathogen interactions, there is also evidence of the disruption of photosystem II functioning during plant–bacteria interactions resulting from the action of effector molecules. According to Torres-Zabala et al. [40], effector molecules produced by Pseudomonas syringae reprogramme the expression of nuclear-encoded chloroplast genes and inhibit photosynthetic CO₂ assimilation through the disruption of photosystem II. This type of acting may well extend to oomycetes.

In the group of four proteins associated to glyoxylate and dicarboxylate metabolism included in the red cluster, 75% of quantifications have more (71%) or neutral (4%) abundance for the control samples, while 75% of quantifications for the inoculated samples have less (58%) or neutral (17%) abundance.

Therefore, we hypothesize that cork oak leaves point to a physiological state associated with an immune response resulting from the interaction with P. cinnamomi, with a focus on protein synthesis, photosynthesis, and glyoxylate and dicarboxylate metabolism biological processes when they contain:

–

A more or neutral abundant level of the proteins P16181, P42798, Q9SIM4, Q9SVR0, Q9LXG1, and Q9FLN4;

–

A less or neutral abundant level of the proteins P56761 and P56778;

–

A less or neutral abundant level of the proteins PMDH2, AACT1, Q9LRR9, and P10795.

4. Conclusions

The α, β-k-Feature Set Problem approach, applied to the group of proteins (80) previously identified as being associated with the cork oak immune response to inoculation by P. cinnamomi, allowed finding several subsets of proteins (solutions of the problem) that separate the control plants from the inoculated plants. These protein subsets are examples of molecular markers of the homeostatic state of cork oak plants that are in interaction with the oomycete. The best solution found (S-R3) includes 29 proteins highly connected between them in a network of protein–protein interactions, comprising three clusters uncovered through STRING analysis. These clusters are linked to each other through pivotal proteins that act as eventual regulators of the functioning of the whole group. Particularly noteworthy are the proteins Q9FLN4-50S ribosomal protein L27 (RK27), PMDH2-Malate dehydrogenase (FAKDZ4/PMDH2), P10795-ribulose bisphosphate carboxylase small chain 1A (RBS1A/RBCS1A), and A0A1P8B485-protein translocase subunit Sec A (AGY1). Interference with the production of these proteins may compromise protein synthesis, photosynthesis, chloroplast biogenesis, and glyoxylate and dicarboxylate metabolism.

An early diagnosis of the biochemical changes that occur in cork oaks in natural stands associated with infection by P. cinnamomi will contribute to the recognition of this pathogen as a real threat to cork oak ecosystems and will also allow the identification of trees with immunological profiles of resistance or tolerance. Monitoring the vegetative status of trees using protein molecular markers can greatly assist in the management of cork oak stands, helping producers to identify plants that interact with P. cinnamomi but do not die suddenly or show observable signs of decline (asymptomatic). We hypothesize that these 29 markers may be decisive for detecting trees that go into decline due to interactions with the pathogen, assisting the management of cork oak forest ecosystems.