Prediction of Thylakoid Lipid Binding Sites on Photosystem II.

The thylakoid membrane has a unique lipid composition, consisting mostly of galactolipids. These thylakoid lipids have important roles in photosynthesis. Here, we investigate to what extent these lipids bind specifically to the Photosystem II complex. To this end, we performed coarse-grain MD simulations of the Photosystem II complex embedded in a thylakoid membrane with realistic composition. Based on >85 μs simulation time, we find that monogalactosyldiacylglycerol and sulfoquinovosyldiacylglycerol lipids are enriched in the annular shell around the protein, and form distinct binding sites. From the analysis of residue contacts, we conclude that electrostatic interactions play an important role in stabilizing these binding sites. Furthermore, we find that chlorophyll a has a prevalent role in the coordination of the lipids. In addition, we observe lipids to diffuse in and out of the plastoquinone exchange cavities, allowing exchange of cocrystallized lipids with the bulk membrane and suggesting a more open nature of the plastoquinone exchange cavity. Together, our data provide a wealth of information on protein-lipid interactions for a key protein in photosynthesis.

Introduction

PSII is a major player of the photosynthetic machinery, the process by which light is converted into chemical energy in plants, algae, and cyanobacteria. PSII functions as a homodimer in vivo; each monomer consists of 27 subunits in plants and 20 in cyanobacteria, respectively. A large number of cofactors supplement PSII with its light harvesting and water splitting capabilities (1, 2).

PSII is located in the thylakoid membrane, together with photosystem I (PSI), cytochrome b6f complex, and ATP synthase, the other major proteins involved in photosynthesis. Over the last few decades, it has become apparent that there is a complex interplay between lipids and proteins. Lipids are vital for the functioning of proteins, and they have structural as well as functional roles (3). The thylakoid membrane is atypical as it has a rather unusual composition, containing mainly galactolipids instead of the more common phospholipids. The predominant galactolipids found in the thylakoid membranes of higher plants and cyanobacteria are digalactosyldiacylglycerol (DGDG), monogalactosyldiacylglycerol (MGDG), and sulfoquinovosyldiacylglycerol (SQDG), together with phosphatidylglycerol (PG) (4, 5). Although the role of these lipids is poorly understood, it is clear that they are essential for the photosynthetic process. First of all, the composition of the thylakoid membrane is conserved in evolution (5), indicating the importance of the galactolipids for photosynthesis. Secondly, experiments show that the lipid composition influences the dimerization state of the protein and also the attachment of the inner antenna proteins, which are important for the electron transport within PSII (6).

MGDG is the most abundant lipid in the thylakoid membrane and is present in the crystal structures of both PSI and PSII (7, 8, 9). It is a non-bilayer-forming lipid that prefers the inverted hexagonal phase when isolated (10, 11), which might be vital for the structural flexibility of the thylakoid membrane (12) and important for the functioning of the violaxanthin cycle (13, 14). MGDG is found at the dimer interface of PSII and it promotes the dimerization of PSII (6, 9). MGDG has additional roles in the functioning of the antenna complexes (15, 16), and is important for chloroplast development (17).

DGDG is also a major component of the thylakoid membrane and is found in the crystal structure of PSII, which contains five DGDG molecules (8). Like MGDG, DGDG has a major influence on structure and functioning of the antenna complexes (16, 18, 19). Absence of DGDG leads to the detachment of the external subunits from PSII, resulting in a decreased oxygen-evolving capacity of the photosystem (20). The existence of hydrogen bonds between DGDG and PSII tyrosines has been shown (21). DGDG has been further shown to influence the thylakoid membrane ultrastructure (22, 23, 24, 25) and lipid packing, the latter playing a role in heat resistance (26).

Although PG is only a minority component of the thylakoid membrane, it has a number of important roles. PG lipids sustain the electron transport activity of PSII, possibly through induced conformational changes at the QB site (6, 27, 28, 29, 30, 31). PG has also been implicated in PSII dimerization (32, 33) via stabilization of the attachment of the external subunits PsbO, PsbU, and PsbV, which are necessary for dimerization of PSII; a lack of PG could therefore indirectly lead to the monomerization of the PSII complex (34, 35). PG is furthermore involved in the stability of PSI and the antenna supercomplexes (16, 18, 19, 36, 37). During phosphate starvation, PG can be replaced by DGDG (38) or SQDG (39).

The anionic SQDG is also relevant for photosynthesis, with four molecules being found in the crystal structure of PSII from cyanobacteria (8). Absence of SQDG results in impaired electron flow from water to tyrosine Z, possibly due to conformational changes in PSII caused by the absence of SQDG (40, 41, 42). SQDG furthermore influences the assembly of the PSII subunits, including the antenna complexes (6, 16).

In addition to the headgroups, the lipid tail composition of the thylakoid membrane is highly regulated, in particular the number and position of double bonds. The saturation grade of the thylakoid membrane differs per organism (5). The presence of unsaturated lipid tails is especially important for photosynthesis (40, 43, 44). Simulation studies (45) show that the unsaturated fatty acids keep the thylakoid membrane of plants in a fluid phase at physiological temperatures. For cyanobacteria, which function at elevated temperatures, this is not required.

Taken together, the experimental data show that thylakoid lipids can exert their effect at different levels of the photosynthetic process. This includes influencing the dimerization of PSII monomers, stabilizing the interactions between the various subunits, facilitating the electron transport, and maintaining the structural integrity of the complex toward the reorganization of antenna complexes. However, due to the difficulty to directly probe lipid-protein interactions, it remains unclear how the thylakoid lipids affect the functioning of PSII (6). Computer simulations are a useful tool to complement experiments in this regard (46, 47, 48). Previous MD simulations have revealed preferential lipid-protein interactions and lipid binding sites on a large variety of membrane proteins (49, 50, 51, 52, 53, 54, 55, 56). PSII has also been subjected to computational modeling (57, 58, 59, 60, 61), but the interaction with the thylakoid lipids has not been investigated so far.

Here we fill this gap based on coarse-grain (CG) MD simulations of PSII embedded in a realistic model of the thylakoid membrane. By omitting atomistic detail, CG models allow us to probe the behavior of a system on much longer timescales (62). Our model is based on the cyanobacterial PSII structure, which was resolved by x-ray crystallography by Umena et al. (8). The PSII cores of plant and cyanobacteria are highly similar; any large differences are mainly found in the lumenal subunits and in the antenna systems (63). The thylakoid lipid composition also differs slightly in plants and cyanobacteria (5, 45). We previously reported on the overall dynamics of the cyanobacterial PSII complex including all of its cofactors (64) and on the exchange of the electron carriers plastoquinone and plastoquinol (65). Here we focus on the interaction of the PSII complex with the thylakoid lipids. We find an overall enrichment of MGDG and SQDG lipids around the complex, and identify a number of potential lipid binding sites on the membrane-exposed surface of PSII. We discuss the nature of these binding sites and their biological relevance in relation to the available experimental data.

Methods

System setup

The MD trajectory analyzed in this article is the same as described previously (64). (For a more detailed description of the system setup, we refer the reader to the Methods of van Eerden et al. (64).) In short, the simulations were based on the PSII x-ray structure of the thermophilic cyanobacterium Thermosynechococcus vulcanus, PDB: 3ARC (8). All cocrystallized lipids were included in the system, as well as all cofactors. This structure was embedded in a realistic representation of the thylakoid membrane consisting of 2686 lipids using the Insane script (66); see Fig. 1. The membrane is composed of 18:1 (9)–16:0 MGDG (40 mol %), 18:1 (9)–16:0 DGDG (25 mol %), 18:1 (9)–16:0 SQDG (15 mol %), di16:0 SQDG (10 mol %), and 18:1 (9)–16:0 PG (10 mol %) lipids, following the experimental data on thylakoid membrane composition of cyanobacteria (5). The system was solvated with 73,144 CG water molecules (representing four times as many water molecules), then 455 Na+ and 455 Cl− ions were added, corresponding to ∼100 mM Na+Cl−, plus 1032 Na+ counterions to neutralize the overall charge of the protein and the negative charge on the PG and SQDG headgroups.

Figure 1

System setup. (A) Given here is an overview of the simulation protocol. The left panel shows a snapshot of the simulation box of PSII embedded in the thylakoid membrane, with cofactors and solvent omitted for clarity. Lipid headgroups are colored as follows: PG, green; DGDG, blue; MGDG, red; SQDG, yellow. The center panel shows a side view of the equilibrated PSII dimer. When the lipid positions are averaged over the whole simulation (87 μs), lipid occupancies can be calculated and regions of high lipid occupancy become visible as lipid binding sites, as shown in the right panel. (B) Shown here are atomistic structures of the thylakoid lipids with the coarse-grained beads superimposed; aliphatic hydrogens are omitted for clarity. To see this figure in color, go online.

The interactions were modeled with the Martini force field (67), using version 2.2 for the protein complex (68) in conjunction with the ElNeDyn elastic network to stabilize the tertiary structure of the protein subunits (69). The lipid parameters for the galactolipids were taken from (70) with the modifications as described in (45). The GROMACS 4.5.5 MD package (http://www.gromacs.org/) was used to perform the MD simulations with the common settings defined for the Martini force field (71). The system was simulated in the isothermal-isobaric (NpT) ensemble. The temperature was set to 328 K using the V-rescale thermostat with a coupling constant τ = 2.0 ps (72). The pressure was coupled semiisotropically to an external bath of p = 1 bar with a coupling constant of τ = 1.0 ps and a compressibility of χ = 3.0−4 bar−1 using the Berendsen barostat (73). Electrostatic interactions were calculated using a shifted potential with a cutoff of 1.2 nm in conjunction with a dielectric constant of 15. The van der Waals interactions were also calculated using a shifted potential, with a cutoff of 1.2 nm and a switch at 0.9 nm. The original simulation of 60 μs (64) was extended to 88.7 μs. The first microsecond was discarded as equilibration, yielding a total simulation time used for analysis of 87.7 μs.

Lipid density maps and definition of binding sites

To identify preferential enrichment of certain lipids around the protein, the average lateral density distributions of the lipids were calculated—the so-called “density maps”. The relative densities of the thylakoid lipids around the protein were calculated using the g_mydensity tool (74). To prevent artifacts from the periodic boundary conditions, jumps were removed from the simulation and the simulation box was quadrupled. The quadrupled trajectory was then centered on one of the four PSII subunits, and subsequently translationally and rotationally fitted on the centered PSII. The positions of the GL1 glycerol beads were used to represent the lipids. Both membrane lipids and cocrystallized lipids were included in the analysis. For each leaflet, the densities were calculated separately. The densities were normalized according to the average density of the lipid species in the entire membrane plane.

Potential binding sites were identified by calculating the time-averaged occupancies of headgroups throughout the simulation box. This was done using the Volmap tool of the VMD visualization package, which overlays a grid on the system and counts the occupation of each cell. The grid resolution was set to 0.2 nm and counting was performed after fitting the simulation on the protein backbone beads (75). The threshold for showing occupancies was based on the amount of MGDG binding sites and adjusted for the other lipids. The MGDG densities reveal clear spots in the simulation box with a threshold value of 36% (implying that at least in 36% of the frames a MGDG headgroup is present, which corresponds to roughly seven times the occupancy of MGDG in the bulk). We adjusted this value for the other lipids based on their relative abundance in the membrane and their headgroup size, which resulted in the following thresholds: 6% for PG, 45% for DGDG, and 30% for SQDG headgroups. Please note that not all grid points in the simulation box are occupied by a lipid in every single frame. Therefore, the bulk occupancy is not the same as the bulk concentration. Only binding sites that were observed in both monomers were selected. The last snapshot of Fig. 1 illustrates the iso-occupancy surfaces thus obtained.

Composition of binding sites

The amino acid composition of the binding sites was determined according to the procedure suggested by Arnarez et al. (55). At 10 ns intervals, for every lipid we created a so-called “phrase”, which is a list of protein residues found within 0.8 nm of the lipid headgroup at that time point. We only considered phrases that contain at least five residues. Subsequently the phrases were clustered with a GROMOS-like algorithm (55, 76). At the start of the clustering procedure all phrases are pooled. Phrases of which 70% of the residues are in common are clustered and the biggest cluster is removed from the pool. This is repeated until all phrases are clustered. Residues that are present in at least 80% of the phrases in a cluster are selected as part of the binding site. All binding sites were inspected manually. If it appeared that residues did not actually contribute to the site, but were merely a “bycatch” of the algorithm, they were removed from the binding site. Two examples are residues in the periphery of a binding site that only contribute marginally to the binding of lipids or residues that are part of a flexible loop, which occasionally diffuses toward the binding site, and thereby only temporarily becomes part of the binding site. An estimate of the uncertainty of the composition was calculated by considering three blocks of ∼29 μs as quasi-independent samples as well as each monomer, resulting in a sample size of six. Please note that this is not a robust error estimate; due to the long correlation times present in the system (exceeding the total simulation time), independent block averages could not be obtained. Consequently, we note that the composition of many binding sites is not completely identical over the two monomers. For the final binding site composition, we combined the compositions of the binding sites of the two monomers to improve our statistics, e.g., if the composition of site X in the one monomer was Arg, Lys, and Tyr, and in the other monomer Lys, Tyr, and Asn, the final composition became Arg, Lys, Tyr, and Asn.

We normalized the amino acid composition of the binding sites by dividing the result of the procedure above by the composition of the residues that come into contact with any lipid headgroup during the simulation. An estimate for the latter was obtained by counting all residues that are at some point of the simulation within 0.8 nm of any (PG, DGDG, MGDG, or SQDG) lipid headgroup. In the normalized distribution, an occurrence >1 corresponds to a preference for a certain residue in a binding site, whereas an occurrence <1 would mean the residue is avoided in the binding site.

Lipid residence times

To calculate the lipid residence times, we monitored, with a frequency of once per nanosecond, whether a lipid was in contact with a residue from a binding site. We used a double cutoff method in combination with a threshold. The double cutoff method is used to prevent the so-called “rattling in a cage” motion that would yield artificially high counts of short residence periods; we used an inner cutoff of 0.6 nm and an outer cutoff of 1.6 nm. A lipid was considered to be in contact with the binding site once it was within the inner cutoff of any binding site bead, until it moved further away than the outer cutoff of any bead in the binding site. We also used a threshold that the lipid headgroup has to be in contact with at least 20% of all the binding site beads to remove lipids that only bind peripherally to the binding site. To filter out membrane lipids that just bypass the binding site without actually binding, we only included residence times of at least 100 ns in the calculation of the average residence time and its standard deviation. Lipids that had been bound for >50 μs and were still bound at the end of the simulation were also excluded from the average value, because their residence time represents only a lower limit.

Results

We simulated the PSII dimer embedded in the thylakoid membrane over a time period of 87.7 μs. This timescale is long compared to the diffusion rate of individual lipids in the thylakoid membrane; a displacement of 1 nm takes, on average, ∼5 ns (45). The long timescale allowed us to analyze the preferential binding of thylakoid lipid types to the PSII complex. We first present results on the specific enrichment and depletion of these lipids based on density maps. We continue to identify specific lipid binding sites based on density iso-surfaces, and characterize the binding sites in terms of residue content and lipid residence times. Finally, a qualitative description is made of the behavior of the cocrystallized lipids.

Please note that we will use the stromal and lumenal nomenclature to refer to the inside and the outside, respectively, of the thylakoid membrane. Strictly speaking, this is only valid for plants, as the cyanobacterial thylakoid membranes are directly located in the cytoplasm and not, as in plants, inside the chloroplast.

MGDG and SQDG enriched around PSII

To probe the preference of PSII for specific lipids, we computed lipid density maps (see Methods). The density maps for the four classes of lipids (MGDG, DGDG, PG, and SQDG) are shown in Fig. 2. Remarkably, the density maps are quite different for each lipid. MGDG and SQDG are clearly enriched around the protein, showing several regions of high densities in the annular lipid shell. This is true for both the stromal and lumenal leaflets. PG and DGDG, on the other hand, are depleted from this region, with the exception of some high-density spots inside the protein that originate from cocrystallized lipids. In particular, PG appears to be excluded from the annular lipid shell, but also to have a somewhat increased density in the subsequent lipid layers. The densities of the cocrystallized lipids in the plastoquinone exchange cavities are also visible. Interestingly, for all lipids but PG the densities inside the cavities are connected with those in the membrane (see the arrows in Fig. 2), implying that lipids are able to exchange between the plastoquinone exchange cavity and bulk membrane.

Figure 2

Normalized density maps of the four different thylakoid lipids. A density of 1 corresponds to the average (bulk) density in the membrane plane. Densities have been capped at two times the bulk density. A cartoon representation of PSII is included to indicate its approximate position and orientation. The location of the plastoquinone exchange cavity is indicated by white lines. Arrows indicate connection points between the high lipid density regions in the plastoquinone exchange cavity and in the membrane; for clarity, the white lines are omitted here. To see this figure in color, go online.

MGDG and SQDG occupy distinct binding sites at membrane-exposed surface

To identify possible binding sites, we calculated iso-occupancy surfaces of the lipid headgroups around the PSII dimer complex. Using a threshold occupancy of five times the bulk density, we obtained clear spots of preferred interaction of lipids with the protein (see Methods). Spots that were present symmetrically, i.e., observed in both monomers, were classified as lipid binding sites. Using this procedure, we found 13 distinct binding sites at the membrane-exposed surface of PSII, in addition to the sites reflecting the cocrystallized lipids. All sites are indicated in Fig. 3 and listed in the Supporting Material (Tables S1 and S2).

Figure 3

Identification of lipid binding sites. Given here are stromal and lumenal views of the PSII backbone with lipid binding sites indicated. PG binding sites are depicted in green, DGDG in blue, MGDG in red, and SQDG in yellow. Binding sites resulting from cocrystallized lipids are numbered in black and refer to Table S2. Binding sites originating from membrane lipids, and which are present in both monomers, are labeled in blue. Labeling of the sites is given according to lipid type (first letter: M, MGDG; S, SQDG), membrane side (second letter: L, lumenal; S, stromal), and sequential numbering. To distinguish between the mirror sites in the two monomers, binding sites in the left monomer are marked with an “∗” and in the right monomer with a “#”. To see this figure in color, go online.

In line with the enrichment of MGDG and SQDG, and depletion of DGDG and PG around the protein (Fig. 2), peripheral binding sites are only found for MGDG and SQDG, with nine and four sites, respectively. There are a few spots where PG and DGDG bind, but they are never present in both monomers and therefore do not match our definition of a binding site. For SQDG, three out of four binding sites are located at the stromal side, as are all the SQDG cocrystallized lipids. MGDG has binding sites on both sides of the protein (five stromal versus four lumenal), although the majority of its cocrystallized lipids are located on the lumenal side.

The binding sites appear fairly homogeneously distributed around the entire membrane-exposed surface of the protein. However, with the exception of a SQDG lumenal binding site (SL1), there are no binding sites around subunits Cyt b 559α, PsbJ, ycf12, and PsbZ. The slightly higher mobility of these subunits (64) might prevent the formation of a well-defined binding site in this region.

Charged residues important for binding

To characterize the nature of the membrane-exposed binding sites, we analyzed the residues comprising those sites (see Methods). The full list of residues for each site is shown in Table S1. The statistical distribution of these residue types is shown in Fig. 4, as well as the normalized distribution corrected for the prevalence of the residues in the protein.

Figure 4

Binding site composition. Given here are bar plots with the amino acid and cofactor composition of the MGDG (red) and SQDG (yellow) binding sites. (A) Shown here is the relative occurrence in terms of percentage. In total, there are 81 residues comprising the MGDG sites and 39 comprising the SQDG sites. This means that an occurrence of respectively 1.2% for MGDG and 2.6% for SQDG sites corresponds to the presence of only one amino acid of that type across all the binding sites. (B and C) Given here is the normalized occurrence of amino acids and cofactors. The occurrence of heme, Ca2+, and Mg2+ is shown in a separate plot with a different scale to improve readability. Error bars represent an estimate of the uncertainty (see Methods). To see this figure in color, go online.

Charged residues appear to be important components of the binding sites; in fact, 12 out of 13 binding sites contain at least one charged residue. Consistently, the binding sites of SQDG (−1e charge) only contain cationic residues (Arg, Lys), whereas sites of MGDG (neutral) contain both negative and positively charged amino acids (Arg, Lys, Asp, Glu). When considering the normalized amino acid composition, in the MGDG binding sites the two positively charged amino acids (Arg and Lys) are equally present. However, this is not true for the negatively charged residues, where aspartic acid is more abundant than glutamic acid in MGDG sites. The importance of electrostatic interactions is further underlined by the occurrence of the cofactors chlorophyll a (CHL), heme, and cocrystallized ions in the binding sites. CHL does contain charged beads, but has a net charge of zero. This might explain why CHL is prevalent in the sites of the neutral MGDG rather than in those of the negatively charged SQDG lipid. Apart from this, some of the binding sites containing CHL, e.g., MS1 and ML1, are located toward the center of the bilayer preventing binding of charged lipids such as SQDG. Heme carries a net charge of −2e, and is not found in the SQDG binding sites. The SL1 sites, however, do contain a Ca2+ and a Mg2+ ion, which form a salt bridge with the SQDG lipids.

The polar threonine and the aromatic phenylalanine residues are enriched in both MGDG and SQDG binding sites, as well as tryptophan, serine, and glycine. Threonine and serine are both polar amino acids with clear hydrogen bonding capacities, and therefore interact favorably with the many hydroxyl groups of the MGDG and SQDG lipids. Although asparagine and glutamine are polar as well, asparagine only occurs in MGDG and glutamine in SQDG sites. Aromatic amino acids prevail at the interface between a polar membrane and aqueous solvent, which is also the location of the binding sites. One could therefore argue that aromatic residues are enriched in the binding sites, just because they share the same location in the protein. However, in the normalized distributions, the aromatic residues stand out, especially tryptophan, and in the SQDG sites, tyrosine as well. It is unlikely that the hydrophobic effect is solely responsible for the enrichment of aromatic residues in binding sites, because there are many other hydrophobic residues present at the membrane interface. It remains therefore unclear what drives the favorable interaction between these aromatic residues and the thylakoid lipids.

Residues like leucine and glycine are, in absolute number, a big contributor to the binding sites. After normalizing for their occurrence in the protein, the contribution of leucine and glycine to the binding sites seems less prominent and more based on their ubiquity in the membrane-exposed surface of the protein than due to specific interactions with the thylakoid lipids.

Lipid residence times vary by orders of magnitude

To characterize the kinetics of the lipids in the binding sites, we analyzed the lipid residence times. To quantify the residence times, we counted how long lipids stay within 1.6 nm of the binding site from the moment the lipids entered within a 0.6-nm cutoff radius of the binding site (see Methods). The distribution of residence times of the MGDG and SQDG lipids is shown in Fig. 5 A. A breakup per individual binding site can be found in Fig. S1.

Figure 5

Distribution of lipid residence times. Logarithmically binned histograms of lipid residence times are shown for the MGDG (left) and SQDG (right) binding sites. Gray bars indicate bins containing lipids that were still bound at the end of our simulation and had been bound for at least 50 μs. In (A), the total count of binding events of a certain duration are shown. In (B), the value of each bar represents the summed time of the binding events in the respective bin (total bound time, i.e., the count multiplied with the residence time), giving proportional prominence to longer residence-time events, which are fewer in number but can represent most of the simulation time.

We observed a broad range of residence times, varying over at least three orders of magnitude from 100 ns to at least 86 μs (i.e., limited only by the length of our trajectory). The vast majority of MGDG and SQDG lipids have residence times below 1 μs. Using a threshold residence time of 100 ns to filter out transiently passing membrane lipids, we calculate an average residence time of 0.5 μs for both MGDG and SQDG, with SDs of 2.5 and 1.9 μs, respectively. There are only a few lipids that remained in their binding pocket longer than 30 μs, and most of them are MGDG lipids. One MGDG lipid stayed bound during the entire trajectory, implying a residence time >86 μs. For SQDG, the longest binding time observed is 74 μs. Considering the total bound time of the lipids (obtained by multiplying the number of bound lipids in a given residence time interval with the residence time; Fig. 5 B), it is clear that a significant fraction of the total bound lipid time can be attributed to lipids with long residence times.

Most of the binding sites bind multiple lipids for a relatively short time, e.g., MS3, MS4, MS5, ML2, ML3, ML4, SS1, and SS2. There are, however, two binding sites that bind a single or a few lipids for an elongated time, sites MS1 and ML1, depicted in Fig. 6. These two sites are both located at the dimer interface. In both of these sites, CHLs play an important role in the binding of the lipids. Whereas MS1 mostly coordinates a single lipid, ML1 is occupied by two or three lipids. The sites MS2, SS3, and SL1 form an intermediate category that bind lipids for a significant amount of time, but at the same time bind many lipids for a short time. SL1 is of particular interest, because it is located at plastoquinone (PLQ)/plastoquinol (PLQol) exchange channels I and III (9, 65) and might therefore have an influence of the PLQ/PLQol exchange (Fig. 6). It is also the only binding site to contain ions. For a more detailed description of every binding site, we refer the reader to Fig. S2.

Figure 6

Visualization of key binding sites. Snapshots of binding sites MS1 (A), ML1 (B), and SL1 (C) are shown. The MS1 site consists of a MGDG lipid coordinated by chlorophylls (green) and various amino acids, the ML1 sites contain three MGDG lipids coordinated by amino acids and a single chlorophyll, and the SL1 site contains three SQDG lipids coordinated by amino acids and ions. MGDG headgroups are colored red and SQDG yellow; lipid tails are colored dark gray. The lipids are shown as sticks and the proteins are shown in white cartoon representation, with the binding-site amino acids as pink sticks, the CHLs in ball-and-stick representation in green, and the ions as orange spheres. To see this figure in color, go online.

Dynamics of cocrystallized lipids

Our simulation data also allow us to look at the mobility of the lipids that are found in the crystal structure. Based on visual inspection, we determined for all the cocrystallized lipids whether they stay in their binding site or diffuse into the membrane. We distinguish between lipids that show no or very limited movement with respect to their position in the x-ray structure (“none”, “limited”), and between lipids that show larger displacements (“large”) or even escape from their binding site (“escape”). In some cases, escaped lipids get replaced by other lipids (“exchange”). The result of this qualitative assessment is shown in Table S2 for each of the 60 cocrystallized lipids. The average behavior is summarized in Table S3.

From the 60 cocrystallized lipids, we found 16 lipids to diffuse away from their initial binding pocket during the entire simulation. Especially MGDG and DGDG lipids were found leaving their binding site (∼30%), which seems caused by the fact that many MGDG and DGDG cocrystallized binding sites are located at the protein periphery. In six cases, the escaped lipid gets replaced by a similar lipid. In the other 10 cases, the binding site is simply lost. The mobility of the lipids located at the surface of the protein is greatly dependent on how well they are shielded from the membrane. The diffusion of PsbX (64) results in the exposure of lipid sites 18 and 46 to the membrane, resulting in their escape. The mobility of the lipids in the protein interior and at the dimer interface is very limited, and none of these interfacial lipids leaves the dimer interface.

Of particular interest is the behavior of the cocrystallized lipids in the PLQ exchange cavity. Some of these lipids have very limited mobility (lipids 11, 36, 37, 38, 40, 55, 56), whereas others are very mobile (lipids 12, 21, 22, 23, 24); see Tables S2 and S3 and Fig. 3. A few cocrystallized lipids also unbind, and are capable of leaving the PLQ cavity (lipids 35 and 39). In fact, we find multiple lipids to diffuse in and out of the PLQ exchange cavity; see Table S4 and Fig. 7. This is in line with the continuous densities between the membrane and the PLQ exchange cavity shown in Fig. 2. There are two pathways that the lipids use to enter and leave the PLQ exchange cavity. The first pathway (channel I) is the PLQ channel formed by Cyt b 559α and PsbJ (77). The second pathway is a new channel that we recently described as a potential exchange channel for PLQ and PLQol (65). We denoted the new pathway “channel III”, as a second PLQ channel was already described by Guskov et al. (9). Channel III emerges when PsbJ moves toward Cyt b 559α and PsbK and ycf12 move in the opposite direction, creating a tunnel between PsbJ on the one side and PsbK and ycf12 on the other side. The sizes of channels I and III can vary (65), and during this simulation, channel III is larger in size than channel I, allowing more lipids to diffuse in and out of the protein. In total (counting both cocrystallized lipids and membrane lipids), channel I is used four times whereas channel III is used 12 times (Table S4). Out of the in-total 16 exchanges, nine were lipid entrances into the PLQ exchange cavity, whereas the remaining seven were lipid exits out of the PLQ exchange cavity. This corresponds to a flux of, respectively, 51 and 40 entrance and exit events per monomer ms−1. MGDG enters and leaves the protein most often, followed by SQDG and DGDG. Remarkably, PG is not seen entering or leaving the PLQ exchange cavity during the simulation at all.

Figure 7

Lipid dynamics in the PLQ exchange cavity. Given here is a time series of snapshots of two diffusion events of a single MGDG lipid (left) and a SQDG lipid (right) into and out of the PLQ exchange cavity. The upper panels show a stromal view and the lower panels, a side view. The MGDG lipid originates from the stromal leaflet and diffuses through channel I, whereas the SQDG lipid originates from the lumenal leaflet and diffuses through channel III. The three PLQ channels are labeled, as well as the stromal and lumenal side of the protein. The color of the lipid is indicative of the relative time with respect to the first snapshot; the snapshots in the left and right panel are from a timespan of, respectively, ∼10 and ∼30 μs. In the upper-right panel, channel II is not visible and in the lower-right panel, all the three channels are not visible at the given viewpoints. To see this figure in color, go online.

The exchange of lipids results in a change in the number of lipids in the cavity, and alters the composition. At the start of the simulation, each cavity contains three PG, two DGDG, two MGDG, and two SQDG lipids, whereas during the last 20 μs the average composition in each cavity is 3.0 PG, 1.5 DGDG, 3.2 MGDG, and 2.5 SQDG lipids. There is thus an increase of 1.2 lipids in the exchange cavity per monomer during the simulation.

Discussion

Preferential interaction of PSII with MGDG and SQDG lipids

Lipids are of great importance for the functioning of membrane proteins, as has become clear in the ever-evolving picture of protein-lipid interactions (3). Our simulations provide a detailed view on how the complex PSII system interacts with the lipids constituting the thylakoid membrane and thereby help us to understand the functioning of this important protein. The major findings from the simulations are 1) the enrichment of MGDG and SQDG in the annular lipid shell around the protein, and 2) the occurrence of multiple MGDG and SQDG binding spots at the membrane-exposed surface. Some of these sites bind lipids for elongated times, in particular sites MS1, ML1, and SL1, which are occupied by at least one lipid most of the time (see Fig. 6). Other binding spots are in fact more transient binding regions that bind multiple lipids shortly, e.g., sites MS3, MS4, ML2, ML4, and SS1. Combining the peripheral cocrystallized sites, which are mostly stable in our simulations, with the new binding spots, we predict an average number of 31.0 MGDG (28.0 lipids from binding sites; see Fig. S2, plus 3.0 cocrystallized lipids, #41–43; see Table S2) and 13.8 SQDG lipids (12.8 lipids from binding sites, see Fig. S2; plus 1.0 cocrystallized, #45, see Table S2) bound to the membrane-exposed surface of the PSII dimer.

An intriguing question is, what are the driving forces leading to the binding site enrichment of MGDG and SQDG instead of DGDG and PG, two lipid types that are also present in large amounts in the thylakoid membrane? We attribute the relative depletion of DGDG to a combination of steric and entropic factors arising from its large and flexible headgroup, consisting of two sugar rings. MGDG, with only one sugar ring, can more easily approach the protein surface. In particular, chlorophyll has an important role in the coordination of MGDG lipids in some of the binding sites, e.g., MS1 and ML1. The outcompeting of PG by SQDG, when both carry a negative charge, is more puzzling. Positively charged residues are important components of the SQDG binding sites, but potentially PG lipids could also bind. A possible driving force is the general affinity of sugar residues for proteins, as is observed in, e.g., all-atom simulations of proteins in sugar solutions (78) and in CG simulations of proteins interacting with ganglioside lipids (79, 80). Additionally, the lipid tails could play a role. SQDG is the only lipid in our simulated mixture that also occurs with fully saturated tails instead of one saturated and one unsaturated tail. However, further analysis of the SQDG binding sites does not reveal a strong preference for either the fully saturated or singly unsaturated variant (data not shown).

Putative biological role of binding sites

The observation of specific lipid binding sites warrants a discussion on their putative biological role. One possibility is that the lipid binding sites have a supporting or catalyzing function in the creation of PSII superstructures, by functioning as a bridge or glue between different complexes, including the associated antenna complexes as recently shown in the formation of megacomplexes composed of PSI, PSII and a phycobilisome antenna complex (81). Additionally, CG MD simulations have shown the structural importance of lipids in supercomplex formation of respiratory chain complexes in mitochondrial membranes (55, 82).

Stabilization of the dimer structure is another form of steering supercomplex formation. Of the two MGDG binding sites at the dimer interface, the ML1 site has contributions from residues originating from the two different monomers. The lipids in this site might therefore contribute to the dimerization. This is in agreement with experimental studies that have shown the importance of MGDG in PSII dimerization (6). The low mobility of the interfacial cocrystallized MGDG lipids 31 and 32 is also in agreement with the importance of interfacial lipids for dimerization. This might be also true for the low mobility of the interfacial DGDG lipids 1–6. Note that these lipids were detergents in the crystal structure that were modeled as DGDG molecules, but in vivo might actually be MGDGs.

Another possibility is that the lipids bound to the surface of PSII are important to sustain the electron transport in PSII, either by participating directly in the electron chain or by stabilization of other residues and molecules that take part in the electron transport. SQDG deficiency in Chlamydomonas reinhardtii results in a decreased efficiency of electron donation from water to Tyr Z (D1-TYR161) (42). It has been proposed that these changes in the electron donation could perhaps be the indirect result of changed interactions of the extrinsic proteins PsbO, PsbU, and PsbV (35). Although the residues of PsbV are not part of any binding site, during the simulation some PsbV residues interact with the SQDG lipids that reside at the SL1 site. Among the PsbV residues interacting with SQDG lipids are the charged residues PsbV-LYS24 and PsbV-ARG31. These residues thus might play a role in the attachment of PsbV, in an SQDG-mediated manner.

Furthermore, some of the bound lipids are near the PLQ exchange cavity, suggesting a role in steering the PLQ/PLQol exchange pathways. For instance, the SL1 binding site is located at the lumenal side of the PLQ exchange cavity, and might thereby influence the flux of PLQ and PLQol in and out of the cavity. Experiments show that SQDG deprivation results in changes around the QB binding site found inside the exchange cavity (35, 41, 42, 83). The low mobility of the cocrystallized PG lipids 51 and 55, and 52 and 56, which are also located in and around the PLQ exchange cavity, supports a functional role, in agreement with the experiments of Itoh et al. (31).

Lipid exchange in the PLQ cavity

With the exception of the binding sites discussed above, we observed lipids in the exchange cavity to be rather mobile in general, and able to exchange with lipids from the bulk thylakoid membrane. The fluxes of lipids through the PLQ channels are of the order of a few tens of lipids ms−1 per monomer, similar to the exchange rates estimated for PLQ/PLQol in our previous study (65). The diffusion of lipids in and out of the cavity appears therefore coupled to the exchange of PLQ/PLQol in and out of PSII. For PLQ and PLQol, we previously showed that in addition to the known channels I and II, a third channel is used (65). Here we observed lipids to also exchange using channel III, in fact as the dominant route in addition to channel I. Channel II, being the narrowest channel, did not allow lipid flow in our simulations, but cannot be ruled out as an additional lipid channel as the opening and closing of the channels is a stochastic event (65). The high fluxes together with the existence of multiple channels suggest a much more open nature of the PLQ exchange cavity than the word “cavity” implies. PLQ exchange “bay” might be a more appropriate description. In this model, the PLQ exchange bay has a relatively open nature, where helices Cyt b 559α + β and PsbJ prevent the bay from being blocked by larger molecules and proteins. The promiscuous nature of the exchange may prevent high entropic energy barriers that are usually associated with ligand binding.

Limitations and outlook

Here we briefly discuss a number of limitations of our study that are important in the interpretation of our results, and provide an outlook for future studies.

The first limitation is the absence of lipid asymmetry in our simulations. In reality, lipid membranes typically have an asymmetric distribution of lipids between the two leaflets. In PSII, all cocrystallized MGDG lipids are located lumenally and all SQDG lipids stromally (8), indicative of the existence of membrane asymmetry. However, a structure of PSI with LHCI reveals cocrystallized DGDG and MGDG in the stromal leaflet (84). This result would make such an asymmetrical distribution of lipids over the two leaflets more unlikely. Due to lack of more quantitative information on lipid asymmetry in thylakoid membranes, we opted for a symmetric membrane model. The newly predicted binding sites of both MGDG and SQDG lipids are located on either side of the protein. Only when more conclusive data on lipid asymmetry become available can the relevance of these sites be reassessed.

The second limitation is that our simulations have not fully converged. By coarse-graining a system it becomes possible to sample much longer timescales, important to achieve statistically meaningful results on protein-lipid binding. Yet, even though we simulated for 87 μs, which is substantial for current practices, the simulation has not converged. The differences between the two monomers are the proof of this.

The third limitation pertains to the use of a CG model, implying a loss of detail as a tradeoff for the increased sampling speed. For instance, consider hydrogen bond capabilities; although implicitly included in the parameterization of the various bead subtypes, their directionality is lost, potentially impacting the relative binding strength of the glycolipids. For an extensive discussion on limitations of the Martini model, we refer the reader to Marrink and Tieleman (85). To verify our CG model, one could simulate (parts) of the system on an all-atom level; for example, to study the interactions between some of the lipids and the protein in more detail. A complete all-atom model of PSII, including all cofactors and thylakoid lipids, could be constructed using efficient backmapping tools (86) and recently parameterized all-atom force fields of the cofactors and the thylakoid lipids (45, 87, 88, 89).

Furthermore, one should be aware that in this study the PSII crystal structure of Synechococcus vulcanus was simulated, which means that care should be taken when interpreting these results for plants. The PSII cores of plants and cyanobacteria are very similar, but their antenna complexes are completely different. The cyanobacterial phycobilisome binds stromally to PSII, whereas the LHCII complexes in plants bind laterally to PSII, which is also the location of the lipid binding sites. One could therefore expect differences between the lipid binding sites in plant and cyanobacterial PSII. A second difference between plants and cyanobacteria that affects the lipid binding sites is the composition of their thylakoid membranes, which mainly differs in fatty acid composition.

Ultimately, our predictions require experimental validation. The experimental identification of lipid binding sites on membrane proteins is a rapidly growing area, thanks to an increasing number of techniques (90). Traditionally, x-ray techniques may reveal tightly bound, cocrystallized lipids, but advanced mutagenesis studies or chemical cross-linking techniques are also used to probe weaker bound lipids that are not washed away under the harsh crystallization conditions. New techniques such as the use of lipid nanodisks to isolate membrane proteins with their native lipid environment (91) and mass spectrometry (92) hold a lot of promise to further this development. Studies in which photosynthetic complexes are embedded in model membranes, that contain one or more of the thylakoid lipids, could further elucidate the role of specific lipids and their role in electron transport and supercomplex formation.

Conclusions

Taken together, we revealed specific lipid binding sites on a large multisubunit protein complex, PSII, embedded in a realistic lipid environment, the thylakoid membrane. Surprisingly, only MGDG and SQDG lipids form strong contacts with the PSII system, whereas DGDG and PG lipids are largely depleted. We find that charged amino acids and chlorophylls are important components of the lipid binding sites. Furthermore, we showed that the PLQ exchange cavity is much more dynamic than previously assumed, allowing the exchange of lipids along with the electron carriers. In general, our large-scale coarse-grain simulations open the way to improve our understanding of protein-lipid interactions in the complex setting of native membranes.

Author Contributions

F.J.V.E. and S.J.M. designed the research. F.J.V.E. and M.N.M. performed the research. All authors analyzed the data and wrote the manuscript.