Membrane Proteomics of Arabidopsis Glucosinolate Mutants cyp79B2/B3 and myb28/29.

Glucosinolates (Gls) constitute a major group of natural metabolites represented by three major classes (aliphatic, indolic and aromatic) of more than 120 chemical structures. In our previous work, soluble proteins and metabolites in Arabidopsis mutants deficient of aliphatic (myb28/29) and indolic Gls (cyp79B2B3) were analyzed. Here we focus on investigating the changes at the level of membrane proteins in these mutants. Our LC/MS-MS analyses of tandem mass tag (TMT) labeled peptides derived from the cyp79B2/B3 and myb28/29 relative to wild type resulted in the identification of 4,673 proteins, from which 2,171 are membrane proteins. Fold changes and statistical analysis showed 64 increased and 74 decreased in cyp79B2/B3, while 28 increased and 17 decreased in myb28/29. As to the shared protein changes between the mutants, one protein was increased and eight were decreased. Bioinformatics analysis of the changed proteins led to the discovery of three cytochromes in glucosinolate molecular network (GMN): cytochrome P450 86A7 (At1g63710), cytochrome P450 71B26 (At3g26290), and probable cytochrome c (At1g22840). CYP86A7 and CYP71B26 may play a role in hydroxyl-indolic Gls production. In addition, flavone 3′-O-methyltransferase 1 represents an interesting finding as it is likely to participate in the methylation process of the hydroxyl-indolic Gls to form methoxy-indolic Gls. The analysis also revealed additional new nodes in the GMN related to stress and defense activity, transport, photosynthesis, and translation processes. Gene expression and protein levels were found to be correlated in the cyp79B2/B3, but not in the myb28/29.


INTRODUCTION
Glucosinolates (Gls) as natural anticancer compounds are represented by three major classes of chemical structures (aliphatic, indolic, and aromatic; Yan and Chen, 2007;Sønderby et al., 2010). In addition to their anti-carcinogenic activities, they have a distinct role in plant defense against herbivores (Halkier and Gershenzon, 2006;Yan and Chen, 2007) and pathogens (Kissen et al., 2009). The activities are attributed to their hydrolysis products, such as isothiocyanates, thiocyanates, and nitriles (Halkier and Gershenzon, 2006). Gls biosynthesis starts from methionine, tryptophan or phenylalanine to produce aliphatic, indolic, or aromatic Gls, respectively (Yan and Chen, 2007;Sønderby et al., 2010). Briefly, the substrate amino acid is converted to aldoxime, then to aci-nitro compounds, thiohydroximate, and desulfoglucosinolate. After sulfation, the core Gls structure is formed. In aliphatic Gls biosynthesis, the methionine chainelongation and the core structure biosynthesis are under the control of three transcriptional factors MYB28, MYB29, and MYB76 (Yan and Chen, 2007;Frerigmann et al., 2012). In the core pathway, the formation of aldoximes is catalyzed by cytochrome P450s CYP79F1 and CYP79F2, and that of the aci-nitro compounds by CYP83A1 (Grubb and Abel, 2006). Then glutathione S-transferase U20 forms thiohydroximates, which are in turn rearranged to desulfoglucosinolate by UGT74B1 (Sønderby et al., 2010), followed by sulfation by SOT17 and SOT18 to produce intact Gls (Sønderby et al., 2010;Mostafa et al., 2016). Similar for indolic Gls, CYP79B2, CYP79B3, and CYP83B1 are responsible for aldoximes and aci-nitro compounds formation, followed by conversion to thiohydroximates by glutathione S-transferase F10, rearrangement to desulfoglucosinolates and sulfation to indolic Gls by SOT16 (Grubb and Abel, 2006;Mostafa et al., 2016). It is clear that the cytochrome P450s play a central role in the Gls biosynthesis, and these proteins are membrane localized (Neve and Ingelman-Sundberg, 2010).
Several studies have reported the relationship between the Gls biosynthetic pathway and other biological pathways in plants, e.g., amino acid and carbohydrate pathways using CYP79F1 RNAi lines , auxin biosynthesis using cyp79B2/B3 mutant (Zhao et al., 2002) and stress response pathways through environmental perturbation (Martínez-Ballesta et al., 2013). In our previous work, we used Arabidopsis double mutants (cyp79B2/B3 deficient in indolic Gls production and myb28/29 deficient in aliphatic Gls production), and discovered new nodes in the glucosinolate molecular network (GMN) that include stress and defense related proteins like glucan endo-1,3-beta-glucosidase, glutathione S-transferase F7 and glutathione S-transferase F2 and the electron carriers cytochrome B5 isoform C and cytochrome c oxidase subunit 5b-2 (Mostafa et al., 2016). To date, no studies have reported the glucosinolate molecular networks in the membrane proteome context.
Since many known glucosinolate proteins such as the cytochrome P450s are membrane or membrane associated proteins, here we investigated how perturbation of Gls metabolism using the aforementioned mutants affects the Arabidopsis membrane proteome using Tandem Mass Tag (TMT) labeling LC-MS/MS based quantitative proteomics. Analyses of protein interaction networks using STRING and functional enrichment of the identified proteins using agriGO allowed us to discover new nodes and edges in the GMN. With qRT-PCR, we were able to determine the correlation between gene transcripts and membrane proteins in the two mutants.
Together with our published soluble proteomics work (Mostafa et al., 2016), this study enables a comprehensive understanding of the Arabidopsis GMNs.

Plant Genotyping, Growth, and Sample Collection
Arabidopsis thaliana (L.) Heynh ecotype Columbia (Col-0) seeds were obtained from the Arabidopsis Biological Resource Center (Columbus, OH, USA). The seeds of cyp79B2/B3 and myb28/29 were kindly provided by Dr. John Celenza (Boston University, Boston, MA, USA) and Dr. Masami Hirai (RIKEN Plant Science Center, Yokohama, Japan), respectively. The mutant genotyping and chemotyping were reported in our previous study (Mostafa et al., 2016). Seed germination and seedling growth were conducted as previously described (Mostafa et al., 2016). Leaves from 5-week old wild type (WT), cyp79B2/B3 and myb28/29 were collected, frozen in liquid nitrogen and stored at −80 • C. Four replicates were included per genotype, and each replicate contains 2 g leaves pooled from 12 plants.

Protein Extraction and Peptide TMT Labeling
Protein was extracted according to Pang et al. (2010) by grinding the leaf tissues in liquid nitrogen and then homogenizing on ice in 10 mM Tris-HCl (pH 7.4), 10 mM KCl, 1.5 mM MgCl 2 , 10 mM dithiothritol (DTT), 0.5 M sucrose, and 10 mM phenylmethylsulfonyl fluoride (PMSF). The protein extracts were filtered through cheesecloth and centrifuged at 800 g for 10 min at 4 • C. The supernatant was transferred to ultracentrifuge tubes and centrifuged again at 100,000 g for 1.5 h at 4 • C. The formed microsomal membrane was washed with 100 mM sodium carbonate using a glass dounce homogenizer, followed by centrifugation at 100,000 g for 1.5 h at 4 • C. The microsome pellets were rinsed with 500 µl resuspension buffer containing 100 mM HEPES (pH 7), 1% triton X-100 and 0.5 M sucrose, and centrifuged at 800 g for 10 min at 4 • C. Protein was precipitated using 5 volumes ice cold 90% acetone overnight at −20 • C, followed by washing the pellets once with ice cold 90% acetone and twice with ice cold acetone before solubilizing in 7 M urea, 2 M thiourea, 4% CHAPS, and 0.25% Triton X-100. The protein amount was assayed using an EZQ assay kit (Invitrogen Inc., Eugene, OR, USA).
A total of 50 µg protein from each replicate was precipitated with ice cold 90% acetone at −20 • C overnight, followed by 20,000 g centrifugation at 4 • C for 15 min. After washing with ice cold 90% acetone, the pellets were solubilized, reduced, alkylated and digested with modified trypsin (Promega, Madison, WI, USA) at a 1:25 (w/w) ratio for 16 h at 37 • C, followed by TMT labeling according to the TMT 6-plex kit manual (Thermo Scientific Inc., San Jose, CA, USA). The WT replicates were labeled with 126 and 127 tags, cyp79B2/B3 replicates with 128 and 129 tags and myb28/29 replicates with 130 and 131 tags at room temperature for 2 h. After quenching with 8 µl 5% hydroxylamine for 30 min, the labeled samples were combined and lyophilized. Two independent experiments and four biological replicates each sample were performed.

Peptide Desalting, Strong Cation Exchange Fractionation, and LC-MS/MS Analysis
The TMT labeled peptides were desalted on Macrospin C-18 reverse phase mini-column (The Nestgroup Inc., Southborough, MA, USA) and fractionated using an Agilent HPLC 1260 strong cation exchange system as previously described (Mostafa et al., 2016). A total of 12 fractions were collected from each experiment. Each fraction was lyophilized, solubilized in solvent A (0.1% formic acid and 3% acetonitrile), and analyzed using an Easy-nLC 1000 system coupled to a Q-Exactive Orbitrap Plus MS (Thermo Fisher Scientific, Bremen, Germany) according to Mostafa et al. (2016) with minor modifications: The mobile phase gradient was ramped from 2 to 30% of solvent B (0.1% formic acid and 99.9% acetonitrile) in 57 min, then to 98% of solvent B in 6 min and maintained for 12 min. Mass analysis was performed in positive ion mode with high collision dissociation energy. The scan range was 400-2,000 m/z with full MS resolution of 70,000 and 200-2,000 m/z with MS 2 resolution of 17,500. The first mass was fixed at 115 m/z, and 445.12003 m/z (polysiloxane ion mass) was used for real-time mass calibration.

Protein Identification and Quantification
The MS data were searched using Proteome Discoverer 1.4 (Thermo Scientific, Bremen, Germany) against the Arabidopsis TAIR10 database with 35,386 entries. The searching parameters were set to 300 and 5,000 Da as minimum and maximum precursor mass filters, digestion with trypsin with two missed cleavages, Carbamidomethylation of cysteine was set as a static modification, and TMT6plex of N terminus, TMT6plex of lysine, phosphorylation of STY (serine, threonine, and tyrosine) and methionine oxidation were set as dynamic modifications. Precursor mass tolerance was 10 ppm, fragment mass tolerance was 0.01 Da, spectrum grouping maximum retention time difference was 1.1 and false discovery rate was 0.01 at the peptide level. Proteins quantification based on labeled unique peptides intensities and statistical analyses were performed as previously described Mostafa et al., 2016;Sun et al., 2017). The proteomics data were deposited to ProteomeXchange repository (accession number: PXD005781).

String Bioinformatics Analysis and Gene Ontology Enrichment
The relationship between the significantly changed proteins and Gls metabolic pathways (Chen et al., 2011;Mostafa et al., 2016) was analyzed using STRING bioinformatics tool (Baldrianová et al., 2015;Ji et al., 2016;Lim et al., 2017). The resulted networks were visualized in the confidence view relying on gene neighborhood, fusion, co-occurrence, co-expression, literature, and available data. To determine the enriched pathways, we performed Singular Enrichment Analysis (SEA) for the changed proteins and the results were compared using a cross comparison of SEA (SEACOMPARE) in the agriGO database (Silva-Sanchez et al., 2013).

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
To determine whether protein expression levels were correlated with transcript levels, we conducted qRT-PCR of 44 genes selected based on the proteomics data (32 for cyp79B2/B3 and 22 for myb28/28). This list of primers used in qRT-PCR is provided in Supplementary Table 1. Total RNA was extracted using a RNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA) and cDNA was synthesized with ProtoScript R II Reverse Transcriptase (New England BioLabs, Ipswich, MA, USA). qRT-PCR was performed with VeriQuest SyBr and a fluorescein kit (Affymetrix, Santa Clara, CA, USA) using CFX96 (Bio-Rad, Hercules, CA, USA) as described previously (Koh et al., 2012). For each reaction, three technical and three biological replicates were included. Relative expression of the target genes was calculated using the comparative C t method (Applied Biosystems, Framingham, USA). The differences in C t values ( C t ) between the target gene and two internal controls (AT4G34270 and AT5G44200) were calculated to normalize differences in the starting materials. The expression ratios of cyp79B2/B3 and myb28/29 to WT were calculated and compared to the ratios from the protein data using Pearson's r.

cyp79B2/B3 and myb28/29 Membrane Proteomes
Based on the MS/MS spectra of high confidence peptides derived from the WT, cyp79B2/B3 and myb28/29, we identified 4673 proteins in two independent experiments using Proteome Discover (Supplementary Table 2). Out of these proteins, 3,132 were identified in both experiments, while 1,076 and 465 were unique to experiments 1 and 2, respectively ( Figure 1A). A total of 4,655 proteins were available for quantification based on unique TMT labeled peptides, highlighting the high efficiency of labeling. PD enrichment analysis (based on TAIR and Uniprot annotations) of the identified proteins showed 2,171 to be membrane proteins ( Figure 1B and Supplementary Table 2). Comparative analysis of the protein expression changes between the mutants and WT at a fold change cutoff (>1.2 and <0.8), a p < 0.05 and transmembrane domain analysis revealed 93 proteins to be increased ( Figure 1C) and 99 to be decreased ( Figure 1D). Transmembrane domain analysis revealed that 175 out of the 192 differentially expressed proteins contained at least one transmembrane domains. The rest deemed to be membrane associated proteins (Supplementary Table 3). Correlating the changed proteins to those involved in Gls metabolism using STRING showed new nodes and edges (Figures 2, 3). The new nodes can be categorized according to their positions in the network as directly correlated or indirectly correlated to Gls metabolism. They can also be classified according to their biological roles as secondary (stress related) and tertiary (other biological process) connections (Detailed in next sections).  Nine membrane proteins showed common changes between the two mutants relative to WT, with only one protein increased while the other eight decreased ( Table 1). By STRING mapping of the significantly changed proteins (Figures 2,  3), we found seven of the nine proteins represented new connections with the glucosinolate metabolic network (GMN). The role of probable cytochrome c (CYC2) and plastocyanin minor isoform (PETE) in electron transport process (Pesaresi et al., 2009;Welchen et al., 2012) makes them biologically relevant tertiary connections in GMN in a way similar to cytochrome B5 isoform C and cytochrome c oxidase subunit 5b-2 (Mostafa et al., 2016). Photosystem I reaction center subunit IV B (PSAE2), 14-3-3-like protein GF14 nu (GRF7), adenine phosphoribosyltransferase 1 (APT1), alba DNA/RNAbinding protein (F28N24.7) and triose phosphate/phosphate translocator (APE2) form other tertiary nodes. Out of this group, APT1 was the only protein directly connected to the GMN (Figures 2, 3).

Specific Changes of cyp79B2/B3 Membrane Proteins
Sixty-four and 74 membrane proteins showed unique increases and decreases, respectively, in the cyp79B2/B3 mutant ( Table 2). Seventy-seven new nodes were discovered by the STRING mapping of these cyp79B2/B3 proteins to the GMN (Figure 2). It was obvious that perturbation of the indolic Gls metabolism affected a group of stressrelated membrane proteins forming new secondary nodes.
Frontiers in Plant Science | www.frontiersin.org FIGURE 3 | STRING analysis of myb28/29 changed proteins in relation to known proteins in Gls biosynthesis. Known Gls biosynthetic proteins are indicated by red balls, new proteins in GMN are indicated by gray balls, proteins changed in both mutants are indicated by italic labeling, and uniquely changed proteins in myb28/29 are indicated by non-italic labeling. Proteins involved in Gls biosynthesis, stress and defense, and other processes are labeled with green, brown, and violet labels, respectively. Connections strength are proportional to edges thickness as derived from neighborhood, gene fusion, co-occurrence, co-expression, previous experiments, and text-mining information at medium confidence score. Asterisk (*) indicates manual connections based on literature. Double asterisk (**) indicates known nodes in both mutants (Mostafa et al., 2016). Full names of the mapped proteins can be found in the abbreviation and protein name columns in Tables 1, 2. in this side network are bifunctional inhibitor/lipidtransfer protein (At2g45180; which has a proteolytic action) and a tetraspanin-18 (TOM2AH2) with unknown functions.

Gene Ontology Analysis of the Significantly Changed Membrane Proteins
AgriGO enrichment analysis of the changed proteins was conducted at the biological processes (BP), cellular components (CC), and molecular functions (MF) levels. By annotating 147 changed membrane proteins in the cyp79B2/B3 using SEA, we got 302 enriched GO terms for BP (Supplementary  Table 4). From this BP analysis, it was obvious that responses to stimuli including abiotic, chemical and stress were highly enriched in cyp79B2/B3 in addition to transport, photosynthesis and metabolic processes. In myb28/29, the most enriched BP terms were those related to translation process. This observation supported our results concerning the stimuli and translation-related proteins in the cyp79B2/B3 and myb28/29, respectively (Supplementary Table 4). On the level of CC, the high enrichment of membrane GO terms supported the effectiveness of our membrane preparation procedure (Supplemental Figures 2, 5).

Comparison of Protein Expression Data with Transcription Data
To determine whether protein level changes correlated with gene transcription changes, we examined the transcript levels of 32 and 22 genes from cyp79B2/B3 and myb28/28, respectively (Supplementary Table 5). The two mutants exhibited different patterns of correlation. In comparison of cyp79B2/B3 to WT, the genes investigated showed a positive correlation between transcript and protein levels in both direction and degree of expression (r = 0.6579, p = 4.269e −05 ; Supplementary Figure 7). However, in comparison of myb28/29 to WT, the genes did not show correlation between the transcript and protein levels (r = 0.0887, p = 0.6945; Supplementary Figure 7), only three out of the 22 genes showed similar regulation at both transcript and protein levels. For example, At4g13770 encoding cytochrome P450 83A1, exhibited down-regulation in myb28/29 compared to WT (Supplementary Table 5). The difference in the degree of correlation in these two mutants implies that different regulatory mechanisms are involved in the transcriptional and posttranscriptional processes in different genotypes (Marmagne et al., 2010;Koh et al., 2012).

DISCUSSION
As a result of Gls metabolism perturbation, many changes in the levels of soluble (Mostafa et al., 2016) and membrane proteins took place. It was interesting to discover new cytochromes to be involved in the GMN. In addition, several groups of stress and defense-related proteins as well as binding and transport activity proteins were related to the indolic and aliphatic GMNs, in addition to a group of ribosomal proteins in the myb28/29 mutant.

Three New Cytochromes in the Glucosinolate Molecular Network
Cytochromes play a key role in Gls biosynthesis. In aliphatic Gls biosynthesis, CYP79F1 and CYP79F2 catalyze the conversion of chain-elongated methionines to aldoximes, which are metabolized by another cytochrome (CYP83A1) to aci-nitro compounds, precursors of desulphoglucosinolates (Grubb and Abel, 2006). As to indolic Gls biosynthesis, CYP79B2 and CYP79B3 convert tryptophan to aldoximes, that are metabolized by CYP83B1 to form the aci-nitro compounds (Grubb and Abel, 2006). In addition, there is another CYP81F2 catalyzing the conversion of indolic-3-glucosinolate to 4-hydroxy-indolic-3-glucosinolate (Sønderby et al., 2010). Furthermore, CYP71A12 and CYP71A13 can metabolize indolic aldoximes to indole acetonitrile and subsequently indole acetic acid derivatives (Nafisi et al., 2007). In our previous study, we reported cytochrome B5 isoform C and cytochrome c oxidase subunit 5b-2 to be new nodes in the aliphatic and indolic GMNs, respectively (Mostafa et al., 2016). Here we discovered cytochrome P450 86A7 (CYP86A7) in redox reaction and metabolism of fatty    acids (Duan and Schuler, 2005), and cytochrome P450 71B26 (CYP71B26) as new nodes in the indolic GMN. Based on STRING analysis, CYP71B26 is connected to CYP81F2 through a direct edge, while CYP86A7 is connected indirectly to CYP81F2 through lectin family proteins (At5g03350 and At5g18470; Figures 2, 4). Given that their connection to a specific and key enzyme in indolic Gls biosynthetic pathway (CYP81F2) and their expression levels were decreased in the cyp79B2/B3 mutant ( Table 2), it is reasonable to hypothesize that CYP86A7 and CYP71B26 play specific roles in 4-hydroxy indolic-3glucosinolate production (Figure 4). Especially their precursor (indolic-3-glucosinolate) and the product were decreased in cyp79B2/B3 mutant as revealed in our previous study (Mostafa et al., 2016). Also by similarity, we can predict a role for the enzymes in hydroxy indolic-1-glucosinolate production ( Figure 5) as its synthesizing enzymes are not known (Sønderby et al., 2010). The third new cytochrome discovered in this study is a probable cytochrome c At1Gg22840 (CYC2), which plays a role in electron transport process (Welchen et al., 2012). CYC2 is in the shared decreased protein category, forming new connections with aliphatic GMN through ADP/ATP carrier protein 1 (AAC1) and 60S ribosomal protein L15-1 (RPL15A), which is connected to GSTF9, GSTF10 and GSTF11, and with indolic GMN through eukaryotic peptide chain release factor subunit 1-2 (ERF1-2), 60S ribosomal protein L28-1 (RPL28A) and adenine phosphoribosyltransferase 1 (APT1). APT1 is connected to GGP1 and SUR1. Although the CYC2 function awaits for further studies, it might play a role in the conversion of aci-nitro compounds to thiohydroximates.

Stress Related Membrane Protein Changes as a Secondary Result of Glucosinolate Metabolism Perturbation
Plant Gls metabolism is responsive to stress conditions, e.g., temperature and light stress (Martínez-Ballesta et al., 2013), water stress (Khan et al., 2010), salt stress (Guo et al., 2013), and microbial stress (Clay et al., 2009). In our previous study, glucan endo-1,3-beta-glucosidase, glutathione S-transferase F2 and glutathione S-transferase F7 in addition to others as stress-related proteins were found to connect to the Gls pathway (Mostafa et al., 2016). Here we found the levels of 51 stress-related proteins changed significantly in the cyp79B2/B3 mutant and six with changes in the myb28/29 mutant. In the cyp79B2/B3 membrane proteome, a group of general stimuli response-related proteins exhibited significant changes compared to WT (  Figures 2, 4). It is known that GAPC2 participates in the oxidation of glyceraldehydes-3-phophate to glycerate FIGURE 4 | Predicted positions of the directly connected nodes and connected cytochrome nodes on the glucosinolate metabolic pathway. Italic indicates proteins changed in both mutants, * indicates proteins changed in cyp79B2/B3, and ** indicates proteins changed in myb28/29. Red color means increased proteins, green color means decreased proteins. Full names of the directly connected proteins can be found in the abbreviation and protein name columns in Tables 1-3. from which pyruvate is formed. The pyruvate can be converted to acetylCoA for methionine chain-elongation in aliphatic Gls biosynthesis or for synthesis of tryptophan in indolic Gls pathway (Mann, 1987). Both glucosinolate classes were decreased in the cyp79B2/B3 mutant in our previous study (Mostafa et al., 2016) together with GAPC2 in this study. Therefore, the connection between GAPC2 and MYBs in the STRING maps reflects functional relationship and does not necessarily indicate direct physical interaction. Another stress related group showing expression level changes was the salt stress and water deficiency group represented by chloride channel protein CLC-c (Jossier et al., 2010), aquaporin PIP2-2 (Javot, 2003;Tournaire-Roux et al., 2003), annexin D1 (ANN1; Gorecka et al., 2005;Jia et al., 2015; formed edge with GSTF9), early-responsive to dehydration stress protein (Rai et al., 2016), probable aquaporin PIP1-5 (Weig et al., 1997), aquaporin PIP2-3 (Daniels et al., 1994), probable aquaporin PIP1-4 , plasma membrane intrinsic protein 1B (Alexandersson et al., 2005), aquaporin PIP1-3 (Kammerloher et al., 1994), probable aquaporin PIP2-6 (Alexandersson et al., 2010), and aquaporin PIP2-7 (Weig et al., 1997 ; Figures 2, 4). The decreased expression of this group of aquaporins ( Table 2) confirms crosstalk between indolic Gls production and water deficiency enzymes (Khan et al., 2010).
The mechanism underlying such crosstalk is intriguing. The reduction in aquaporins potentiates our observation of retarded growth of Gls mutants (Mostafa et al., 2016). The decreased Gls production resulted in stress status, which led to decreased water uptake and decreased expression of aquaporins, and thus growth retardation.

Effects of Glucosinolate Metabolism Perturbation on Other Processes and Nodes
Gls biosynthetic pathway is organelle specific and involves transport starting from methionine chain-elongation, sulfate transport, and ending with Gls storage in the seeds (Sønderby et al., 2010;Gigolashvili and Kopriva, 2014;Jørgensen et al., 2015). Here we report a decrease in ABC transporter B family member 19 (Lin and Wang, 2005) in both mutants ( Table 1). In addition to their role in sulfate transport, ABC transporters are involved in transporting Gls hydrolysis products (Kang   Figures 2,  4). How this lectin family protein function is not known.
Another biological process affected by the Gls perturbation is photosynthesis as revealed by the increase of photosystem I reaction center subunit IV B in both mutants (Table 1), and increases in cyp79B2/B3 photosystem II stability/assembly factor HCF136 (Meurer et al., 1998), protein curvature thylakoid 1B, NAD(P)H-quinone oxidoreductase subunit H, light-harvesting complex I chlorophyll a/b binding protein 1 and light-harvesting chlorophyll protein complex II subunit B1 ( Table 2). The increased activity in the photosynthetic process could be a strategy to compensate for the internal stress in the mutants as indicated by changes of many stress-related proteins (Tables 2, 3; Mostafa et al., 2016). It was obvious that aliphatic Gls metabolism perturbation activated the ribosomal protein expression as reflected by the increased levels of 18 ribosomal proteins in the myb28/29 ( Table 3). The biological implication of this change is not known although we can correlate it to the regulation of aliphatic Gls biosynthetic pathway by MYB28 and MYB29 . In both mutants, adenine phosphoribosyltransferase 1 (APT1) acting on adenine phosphorylation (Allen et al., 2002) showed connections with GGP1 and SUR1, so it might have a role in thiohydroximate formation (Figures 2-4). Its decrease in levels may be a feedback of the decreased Gls production in the mutants. In cyp79B2/B3, FtsZ homolog 1 (FTSZ1) involved in chloroplast division and protein binding (Osteryoung et al., 1998) was found to connect with BCAT3 and GSTF9, suggesting it may affect methionine chain-elongation and thiohydroximate synthesis. Interestingly, another FtsZ homolog 2-2 (FTSZ2-2; McAndrew et al., 2008) was also connected with GSTF9 (Figures 2, 4). Isoform 3 of dihydrolipoyllysine-residue succinyltransferase component of 2-oxoglutarate dehydrogenase complex 2 (At4g26910) is a member of tricarboxylic acid cycle and can affect methionine biosynthesis and its coupling to acetylCoA in the chain elongation process. Interestingly, it was found to form multiple connections with GMN via BAT5, BCAT3, IMD1, IMD2, IMD3, GSTF9, and SUR1 (Figures 2,  4). In addition, ATP sulfurylase 1 (APS1), a hydrogen sulfide biosynthesis enzyme, formed edges with GGP1 and SUR1, suggesting its potential role in thiohydroximate synthesis (Figures 2, 4). The increased levels of the aforementioned proteins may reflect a feedback mechanism to compensate for reduced Gls levels in the cyp79B2/B3. Flavone 3'-Omethyltransferase 1 (OMT1) in flavonoid metabolism (Muzac et al., 2000) was connected with FMO1, so it could participate in sulfinyl Gls formation (Figures 2, 4). This finding provides another line of evidence for the pathway interaction between phenylpropanoids and glucosinolates. Previously, methionine derived aldoximes were shown to directly or indirectly inhibit caffeic acid O-methyltransferase (COMT) and caffeoyl-CoA O-methyltransferase CCoAOMT), leading to low levels of phenylpropanoid metabolites (Hemm et al., 2003). Here the decreased levels of OMT1 in cyp79B2/B3 may contribute to the decreased production of sulfinyl Gls in the mutant. The data support our metabolomics finding concerning the decreased shikimate level (Mostafa et al., 2016). Another possibility of the OMT1 activity is methylation of hydroxy-indolyl Gls to form methylated indolic Gls (unknown before, Sønderby et al., 2010) in a way similar to methylation of quercetin into isorhamnetin (Figure 5). In myb28/29, 60S ribosomal proteins L13-1 (BBC1) and L15-1 (RPL15A) might be a component in thiohydroximate synthesis through the connections with GSTF9, GSTF10 and/or GSTF11. Both proteins were increased, presumably to compensate for the deficiency of aliphatic Gls in the mutant (Mostafa et al., 2016).

The Proteome and Transcriptome Correlation
In the cyp79B2/B3, the defense and stress-related genes calreticulin 3 (At1g08450; Sun et al., 2014), calmodulin (At2g41100; Cazzonelli et al., 2014), lectin (At5g03350; Armijo et al., 2013), and SNAP25 (At5g61210; Eschen-Lippold et al., 2012) showed significant upregulation in the transcriptome and increases in the proteome. Malate dehydrogenase 2 expression was decreased at both the transcript and protein levels, and it is known to be involved in bacterial defense (Jones et al., 2006). These data have provided additional evidence for the relationship between indolic glucosinolates and stress responses. The overall positive correlation between protein and gene expression levels in the cyp79B2/B3 indicates transcriptional regulation of indole glucosinolates. In myb28/29, although there was no overall correlation between transcript and protein levels, isoform 2 of LysM (At1g21880; Willmann et al., 2011) and AIG2 (avirulence induced gene, At5g39730) exhibited similar downregulation patterns as their corresponding proteins. Both genes are involved in cellular stress responses (Jiang et al., 2007;Willmann et al., 2011). Post-transcriptional and post-translational regulations may contribute to the non-correlation between the expression of some of the genes and their encoded proteins in myb28/29.

CONCLUSIONS
Glucosinolate biosynthetic process is controlled by several cytochrome proteins known to be localized to the membrane, but little is known about how Gls metabolism would affect the membrane proteome. In this study, we aim to address this important question utilizing the TMT labeling based quantitative proteomics of two genetic mutants, i.e., cyp79B2/B3 as the indolic Gls mutant and myb28/29 as the aliphatic Gls mutant. We identified 4,673 proteins, out of which 2,171 were membrane proteins. From these membrane proteins and after transmembrane domain analysis, 192 exhibited different levels relative to WT, with cytochrome P450 86A7, cytochrome P450 71B26 and probable cytochrome c representing new cytochromes potentially involved in GMN. Based on our analyses, the first two might play a role in hydroxyl-indolic Gls production. In addition, a flavone 3 ′ -O-methyltransferase 1 is hypothesized to participate in the methylation process of the hydroxyl-indolic Gls to form methoxy-indolic Gls. GO functional enrichment revealed important processes related to stress response, transport activities and photosynthesis in the cyp79B2/B3 and those related to protein translation in the myb28/29. A transcription profiling of both mutants showed a strong correlation between transcript and protein levels in cyp79B2/B3, and no significant correlation in myb28/29. Overall, the new nodes and edges discovered in the GMNs are useful resources for future hypothesis-testing experiments and ultimately toward engineering and breeding of Gls profiles with positive impacts on human health and plant defense.

AUTHOR CONTRIBUTIONS
IM performed the experiments, data analysis and paper drafting; MY performed qRT-PCR experiment and data analysis; NZ participated in protein extraction and peptides labeling; SG conducted the statistical analysis; CD contributed in LC/MS analysis of peptides; MA and ME provided supervision and advice, and SC designed the experiments, supervised the work and finalized the manuscript.

ACKNOWLEDGMENTS
We would like to thank Chen laboratory members for their support and co-operation. The US National Science Foundation (NSF CAREER 0845162), University of Florida, and the Egyptian Government represented by the Egyptian Cultural and Educational Bureau at Washington DC are acknowledged for funding this project.

SUPPLEMENTARY MATERIAL
The