A different perspective for nonphotochemical quenching in plant antenna complexes

Light-harvesting complexes of plants exert a dual function of light-harvesting (LH) and photoprotection through processes collectively called nonphotochemical quenching (NPQ). While LH processes are relatively well characterized, those involved in NPQ are less understood. Here, we characterize the quenching mechanisms of CP29, a minor LHC of plants, through the integration of two complementary enhanced-sampling techniques, dimensionality reduction schemes, electronic calculations and the analysis of cryo-EM data in the light of the predicted conformational ensemble. Our study reveals that the switch between LH and quenching state is more complex than previously thought. Several conformations of the lumenal side of the protein occur and differently affect the pigments’ relative geometries and interactions. Moreover, we show that a quenching mechanism localized on a single chlorophyll-carotenoid pair is not sufficient but many chlorophylls are simultaneously involved. In such a diffuse mechanism, short-range interactions between each carotenoid and different chlorophylls combined with a protein-mediated tuning of the carotenoid excitation energies have to be considered in addition to the commonly suggested Coulomb interactions.


Parallel Tempering in the Well-Tempered Ensemble
The Parallel Tempering 14 simulation in the Well-Tempered Ensemble 15,16 (PT-WTE) was run with 15 replicas, each started from the last, equilibrated frame of the cMD CryoEM simulation and further equilibrated (10 ns) in the NPT ensemble at the corresponding temperature. The temperature for the i-th replica was determined as: with T 0 = 300 K (the temperature of the first replica) and k = 0.016. After equilibration, a parallel tempering well-tempered metadynamics simulation was run with the potential energy of the system as CV. Gaussians were deposited every 500 steps with a height of 5 kJ mol −1 , a bandwidth of 1000 kJ mol −1 and a bias factor of 20. This short metadynamics ensures a sufficient degree of overlap between the potential energy distributions of the replicas, thereby improving the acceptance probability of Monte-Carlo exchanges between replicas at different temperature during the production run. A 1 µs production run is then started for each replica, for a total of 15 µs of aggregate sampling. In order to preserve an optimal overlap between the potential energy distributions of neighboring replicas, a quasi-static regime for the underlying bias potential was necessary, where Gaussians were added rarely (every 50000 steps) during the course of the production simulation. Walls were included in order to prevent the dissociation of chlorophylls and the neoxanthin from the protein when exploring the phase space associated with higher temperatures. Figure 1: Replica exchanges in the PT-WTE Simulation. a Probability density estimates for the potential energy distributions of each replica in the PT-WTE simulation. Neighboring replicas do show overlapping potential energy distributions. b Replica exchanges along the simulation. c Diffusion in temperature space for each replica along the simulation.

Metadynamics
Metadynamics (MetaD) simulations have been performed with GROMACS 2018.6 17 simulation package patched with the PLUMED library 18 , version 2.5.1 19 . Conversion between AMBER and GROMACS has been handled with ParmEd 20 .
Equations of motion have been integrated using the leap-frog integrator with a 2 fs timestep. Holonomic constraints have been imposed on hydrogen bonds with the LINCS algorithm 21 . A cutoff of 10 Å has been applied for the electrostatic and van der Waals short-range interactions, combined with Particle Mesh Ewald (PME) for the treatment of the long-range electrostatics. The simulations are carried out in the NPT ensemble, with the employment of the velocity-rescaling thermostat from Bussi 22 and the Parrinello-Rahman pressure coupling method 23 . Periodic Bound-ary Conditions (PBCs) have been applied.

Multiple Walkers Well-Tempered Metadynamics
The preliminar multiple walkers well-tempered metadynamics 24,25 has been run using the P1 s and P1 l angles as collective variables (CVs). In this simulation, Gaussians are deposited with a time interval of 1 ps, a height of 1.0 kJ mol −1 , a standard deviation of 0.01 rad and a bias factor of 10. The bias has been stored on a grid to speed up the simulation. Harmonic walls have been applied to the sine of both CVs. Two walkers have been employed to rapidly explore a huge portion of conformational space.

Well-Tempered Metadynamics with Path CVs
An additional well-tempered metadynamics with path CVs 26 has been employed in order to improve on the previous CVs and increase the transitions between the cryo-EM and the Open structures. The path has been constructed from a starting steered MD simulation connecting the two basins, and selecting six landmarks along the MD so as to obtain an average interframe separation of 1.23 Å. The path was then parametrized according to Ref. 26 . We have used the RMSD as a metric to compare the structures from the metadynamics to those of the reference path landmarks. The RMSD is computed by first aligning on the Cα of the residues of helix B,C,A and D, the lutein, and the F211 and W226 residues. Then, the RMSD is computed on the previous selection with the exclusion of helix C. The λ parameter of the path CVs was set to 1.87 Å −2 .

Analysis of MD simulations
The analysis of the trajectories, the featurization, dimensionality reduction and clustering were performed with MDTraj 27 , PyEMMA 28 , Scikit-Learn 29 , and in-house Python scripts.

Analysis of the PT-WTE simulation
The PT-WTE simulation was reweighted prior to the analysis by employing the final bias on the potential energy as a static bias, from which the metadynamics weights have been computed. In addition to the bias on the potential energy, we have included the bias due to the walls on both Chl a612 and Chl a603 in the reweighting. Then, only those frames characterized by a high statistical weight have been employed for the subsequent analysis. In this way, 48.8 % of the frames of the 300 K replica was retained, accounting for the 99.8 % of the total weight of the trajectory. All the distributions reported for the PT-WTE simulation are reweighted ones, i.e., the probability density of an observable ϕ is estimated as: where K(ϕ; ϕ i , h) is a normalized kernel centered in ϕ i with parameter h, and where ω i is the statistical weight of the i-th frame. The probability density for non-periodic variables has been estimated using a Gaussian kernel with bandwidth h, whereas the density of periodic variables has been estimated with a von Mises kernel with a concentration parameter h.

dihedral PCA and clustering
The dimensionality reduction analysis was carried out employing residues of the main helices A and B, and residues on the lumenal side of the complex. More in detail, we have selected residues starting from R101 to D123, belonging to helix B and up to the beginning of helix E, and residues from A201 to T236, which belong to helix A up to the end of helix D. The analysis has been performed using the ϕ, ψ, χ 1 and χ 2 angles of each aminoacid. The cosine and sine of each angle have been employed as features for the following dimensionality reduction step, for a total of 386 features. Dimensionality reduction has been carried out with PCA. . In all cases, the conformational space sampled in the PT-WTE simulation is larger than the one sampled in the cMD CryoEM simulation.
The clustering was performed employing the first 30 principal components of the dPCA latent space. We have used an agglomerative clustering algorithm, employing the Ward's method for performing each merge. The best number of clusters was explored ranging from a minimum of two clusters up to a maximum of nine clusters, and monitoring the Bayesian Information Criterion (BIC), the Calinsky-Harabasz, and the mean silhouette scores. A final number of six clusters was determined from inflections in the scores curves. The final clustering and the clustering scores are shown in Figure 3.

Electronic analysis of quenching
A commonly accepted mechanism of chlorophyll quenching involves the excitation energy transfer (EET) from the Chl's lowest singlet excited state Q y and the dark S 1 state of the neighboring carotenoid, followed by rapid decay to the ground state with dissipation of the excitation energy as harmless heat. Following what is commonly done in the literature, [30][31][32] , we have modeled EET in the Förster weak coupling limit by assuming a constant spectral overlap between the Chl Q y emission and the Car S 1 absorption spectra. In this approximation, the LHC can tune the EET quenching through the modulation of the EET coupling V EET . The EET coupling can be written as a sum of a long-range Coulomb contribution and a short-range contribution:

Estimation of the Coulomb coupling
The Coulomb contribution to the total coupling is approximated as the Coulomb interaction between transition charges (TrEsp method 33 ) of the two interacting pigments: The transition charges for the Q y state of Chls were computed at the TDDFT B3LYP/6-31+G(d) level of theory; those for the S 1 state of Cars were obtained at the DFT-MRCI level of theory. 34 The transition charges are rescaled by a factor of 3.7, analogously to what has been done in Ref 1 , in order to be comparable with previous work. 1,31

Estimation of short-range effects
The short-range contribution V short was estimated for the S 2 /Q y coupling, as the S 2 state is accessible to standard excited-state calculations. V short was computed as the difference between the total coupling V EET and the Coulomb contribution V Coul . The total coupling V EET was obtained from an excited-state calculation of the Car-Chl pair at the TDA/ωB97X-D/6-31+G(d) level of theory followed by the application of the multi-FED-FCD diabatization scheme 35 . This scheme yields the energies of locally-excited (LE) and charge-transfer (CT) states, as well as all couplings between them, includig the EET coupling. Due to the double-excitation nature of the Car S 1 , these calculations cannot provide information about this state. The Coulomb coupling V Coul is obtained via direct integration of the transition densities of the two pigments 36 . These calculations were performed on 100 snapshots extracted from clusters 4,5, and 6 (See below). Short-range effects can be estimated geometrically by resorting to an approximation of the electronic wavefunction overlap in terms of the volume of the intersection between the van der Waals regions of the interacting pigments 37 . These van der Waals regions are defined as the union of interlocking spheres positioned on the conjugated atoms of the two pigments. In order to directly compare with our previous work, we use the same radii (1.4 times the van der Waals radius) defined there. This definition of "geometrical overlap" was sufficient to explain the variability of triplet energy transfer (TET) couplings among different Car-Chl pairs. The couplings involved in TET have vanishing Coulomb contribution due to the spin-forbidden nature of singlet-triplet transitions, 38 and only contain short-range terms. We leverage the relationship between overlap and TET couplings to estimate the variability of the short-range contribution. The TET couplings can be as large as 20 cm −1 when the overlap is above 20 Å 3 , and exponentially drop to ≲3 cm −1 for overlap below 10 Å 3 . 37 . Figure 4: Short-range coupling dependence on the overlap. Each point is colored according to the corresponding cluster (red for cluster 4, purple for cluster 5, brown for cluster 6). Empty circles denote the Lut-a612 pair, while filled circles denote the Vio-a603 pair.

Supplementary
As shown in Figure 4, the dependence on the geometrical overlap of the short-range singlet EET S 2 /Q y is very similar to the case of TET 37 . In addition, we note that the order of magnitude of the short-range couplings is comparable to that of Coulomb couplings, indicating that these short-range contributions might have a significant impact on the total EET coupling.

Determination of the Charge-Transfer states and Marcus analysis
The energy of the CT (Car + Chl − ) state and its coupling with the Q y state of the carotenoid's central chlorophyll were computed with the same protocol employed in Ref. 39 , in order to obtain comparable results. The snapshots employed were extracted from the clusters 4, 5, and 6, by first filtering out structures with a low statistical weight and then selecting points according to the Farthest Point Sampling (FPS) algorithm 40 , so as to obtain a set of geometries representative of each cluster. These snapshots (100 per cluster) were employed for performing QM calculations on both Lut-a612 and Vio-a603 pairs, for a total of 300 calculations per Car-Chl pair. The phytyl tail of the chlorophyll was cut after the first aliphatic carbon and kept within the MM region. The excited states are computed at the TDA/ωB97X-D/6-31+G(d) level of theory. The effect of the environment was included through MMPol. Our multi-FED-FCD diabatization scheme 35 has been employed in order to extract the CT energies and couplings. In order to compare with the results of Ref. 39 , where the energy of the Chl Q y state was computed at TD-DFT/ωB97X-D/6-31+G(d) level of theory, we have subtracted from the Q y energy determined via the diabatization scheme the difference in mean with the Q y energy obtained in LHCII. Excited-state calculations were performed with a locally modified version of the Gaussian package 41 .
The parameters employed for the determination of the charge-separation rates within the framework of Marcus theory were derived following Ref. 39 . In particular, the adiabatic energy of each state was computed as ∆G X = ⟨E X ⟩ − σ 2 X 2k B T , with X = Chl*/Car + Chl − . The reorganization energy λ for the Car + Chl − /Chl* transition was determined from the variance of the energy difference ∆E of the two states as

Carotenoid geometry optimization
A subset of frames from clusters 4, 5, and 6 was selected for the optimizations. Carotenoid geometries were optimized at the DFT/B3LYP/6-31G(d) level of theory within a QM/MM ONIOM scheme. The QM part consisted of the carotenoid molecule (Lut or Vio) only. The MM part was described with the same force field parameters used in the molecular dynamics. MM Residues within 6 Å from the carotenoid atoms were allowed to move during the optimization.

Semiempirical calculation of the S 1 energies
The carotenoid S 1 energy has been computed with semiempirical CI (SECI), following Ref. 42 , where the SECI calculations have proven accurate in predicting the excitation energy changes with geometry. Orbitals for the SECI calculations are determined through a SCF calculation on an open-shell singlet state with two singly occupied orbitals. The SCF calculation is based on the OM2 Hamiltonian 43,44 . The active space comprises all the π and π * orbitals of the conjugated main chain, selected according to the procedure proposed in Ref. 45 . All single excitations from two reference determinants, closed-shell HF singlet and the doubly excited HOMO 2 → LUMO 2 , were included in the multireference CI wavefunction. Definition of the dihedral angles shown in a, c. d Distributions of the first lumenal (d1 l ) and stromal (d1 s ) dihedral for both Lutein and Violaxanthin in each cluster. e Correlation of the dihedral angles (d1 l and d2 l for Lut, d1 s and d2 s for Vio) before (@MD) and after (@OPT) the carotenoid QM/MM optimization). The first dihedral of the conjugated chain, d1 X , X = s, l, is represented by filled circles, while the second dihedral d2 X , X = s, l, is represented by empty circles. Points are colored according to the respective cluster identity (red for cluster 4, purple for cluster 5, brown for cluster 6).