Further thermo‐stabilization of thermophilic rhodopsin from Thermus thermophilus JL‐18 through engineering in extramembrane regions

Abstract It is known that a hyperthermostable protein tolerable at temperatures over 100°C can be designed from a soluble globular protein by introducing mutations. To expand the applicability of this technology to membrane proteins, here we report a further thermo‐stabilization of the thermophilic rhodopsin from Thermus thermophilus JL‐18 as a model membrane protein. Ten single mutations in the extramembrane regions were designed based on a computational prediction of folding free‐energy differences upon mutation. Experimental characterizations using the UV‐visible spectroscopy and the differential scanning calorimetry revealed that four of ten mutations were thermo‐stabilizing: V79K, T114D, A115P, and A116E. The mutation‐structure relationship of the TR constructs was analyzed using molecular dynamics simulations at 300 K and at 1800 K that aimed simulating structures in the native and in the random‐coil states, respectively. The native‐state simulation exhibited an ion‐pair formation of the stabilizing V79K mutant as it was designed, and suggested a mutation‐induced structural change of the most stabilizing T114D mutant. On the other hand, the random‐coil‐state simulation revealed a higher structural fluctuation of the destabilizing mutant S8D when compared to the wild type, suggesting that the higher entropy in the random‐coil state deteriorated the thermal stability. The present thermo‐stabilization design in the extramembrane regions based on the free‐energy calculation and the subsequent evaluation by the molecular dynamics may be useful to improve the production of membrane proteins for structural studies.


| INTRODUCTION
To prepare protein samples in structural biology, a thermostable mutant is often effective especially in the case of membrane proteins including G protein-coupled receptors, 1-3 AMPA receptor, 4 purine/H + symporter, 5 and integral membrane diacylglycerol kinase. 6 Currently, the establishment of such thermostable mutants requires a laborious screening in many cases. Therefore, a rational design of thermostable mutants is a desirable technology in the structural biology of membrane proteins. Recently, a computational method has been developed to identify possible thermo-stabilizing mutations in the intramembrane region of a membrane protein. 7,8 However, there is no report to date for the rational thermo-stabilization design in the extramembrane regions.
It is known that certain thermophilic organisms have "hyperthermostable" proteins tolerable at temperatures over 100 C. 9,10 Recently, we reported a rational hyperthermostabilization design of a soluble globular protein, CutA1, based on a computational prediction of folding free-energy differences upon a site-directed mutagenesis. 11 In the case of CutA1, the denaturation temperature was improved from 86 C to 137 C by introducing two hydrophobic and six ionic mutations; both the enthalpic effect of hydrophobic interactions and the entropic/enthalpic effect of electrostatic interactions contributed to the stabilization.
In this work, the same strategy as that used for CutA1 has been applied to a thermostable membrane protein, toward the design of hyperthermostable mutants. As a model membrane protein, we selected the thermophilic rhodopsin from Thermus thermophilus JL-18 (TR) that had been discovered as the first rhodopsin from thermophilic organisms. 12 The crystal structure of TR is available (PDB entry 5AZD), 13 indicating that the effect of mutation can be evaluated by the molecular dynamics (MD) simulation. Because the hyperthermostabilization strategy was originally developed on a soluble globular protein, we applied it in the extramembrane regions of TR. TR is composed of 260 amino-acid residues with seventransmembrane helices, and functions as a proton pump in response to the light where a covalently-bound retinal is used as a sensing cofactor. By utilizing its high thermostability, TR is expected as a useful tool for the optogenetics that controls an organism by the light. 13,14 Therefore, the hyperthermostable TR should expand its applicability as a thermostable biosensor in biotechnology. Here we report a rational design of further thermostable TR mutants through engineering in the extramembrane regions.

| Design of mutants
The protein model of TR (UniProtKB H9ZSC3) used was prepared using the "Clean Protein" function of the software Discovery Studio v4.1 (BIOVIA; https://www.3dsbiovia.com/) from the chain A of the PDB entry 5AZD (amino acids 3-253) 13 : missing side chains at the residues 89-91, 178, and 253 were mended considering rotamer conformations; the residue 253 was modeled as a C-terminus by attaching the OXT atom. Thermo-stabilizing mutants of TR were designed using the program FoldX (http://foldxsuite.crg.eu/). 15 Before the design, the retinal group attached to Lys233 in the TR model was removed, and then the model was treated by the "Repair PDB" function of FoldX. Twenty-eight solvent-exposed residues in the extramembrane region were selected manually using Discovery Studio: Glu250(C), Glu251(C); characters in parentheses indicate the location of residues where N and C represent the N-terminal and the C-terminal sides of the protein, respectively. In the selection, a residue was excluded in case that: its main-chain dihedral angle was allowed only for limited residue types; it was a proline or an acidic residue for capping an α helix; it was a basic residue at the C-terminus of an α helix; it was an ionic residue making ion-pairs with other residues; a nativestate MD simulation predicted a burial of the residue in the protein molecule or in the membrane. The 28 residues were submitted to the calculation of the folding free-energy using the "Position Scan" function of FoldX. Twenty candidates of single mutations of high scores from the Position Scan were selected and were further analyzed using the  Table 1); a mutant was excluded in case that a substantial structural change upon mutation was suggested by the Build Model.

| Sample preparation
The gene of TR, which was optimized to Escherichia coli codon-usage, was inserted into the pET21c expression vector (Novagen, Madison, WI). The mutations were introduced by the site-directed polymerase chain reaction mutagenesis using the PrimeSTAR Max DNA polymerase (TaKaRa). For the protein expression, transformed E. coli C43 (DE3) cells harboring the plasmid were initially grown at 37 C in 10 mL of LB medium supplemented with ampicillin (the final concentration was 50 μg/mL) and were directly inoculated into 1 L of LB medium containing ampicillin. The cells were grown in a rotary shaker at 30 C after 2.5 hours incubation with 0.5 mM isopropyl β-D- for 20 minutes and at 80 C for 10 minutes so that TR trimers were dissociated into monomers as described previously. 16 The purified samples were diluted to be 0.5-0.8 mg/mL and were stored at 4 C until use.

| Measurement of thermostability
Unless otherwise noted, data was represented as mean ± SE of the mean (s.e.m.) from three independent experiments. To compare the thermostability of wild-type and the mutants of TR, the residual pigment after incubation at 90 C for 4 minutes of the partially purified samples was measured at maximum absorption wavelength of 530 nm by the UV-visible spectroscopy. Next, the residual pigment values at 530 nm of the purified samples in buffer-C were plotted against time for the incubation at 90 C as described previously. 12 The data was analyzed by a single exponential decay function to estimate a decoloration rate constant. We also determined the apparent midpoint temperature of denaturation (T m ) of the purified samples by the differential scanning calorimetry (DSC) as described elsewhere. 17 The DSC experiments in buffer-C were performed on a MicroCal PEAQ-DSC Automated (Malvern Panalytical, Malvern, United Kingdom) ramping from 30 C to 130 C at a scan rate of 1 C/min, for the wildtype and the mutants of TR. The data was processed employing the MicroCal PEAQ-DSC Analysis Software (Malvern Panalytical, Malvern, F I G U R E 1 Location of 10 residues for mutation. The TR model is shown as a ribbon representation with secondary-structure colorings: red, light blue, and green represent helix, sheet, and turn, respectively. The N-and C-termini are labeled. The model is embedded in a lipid bilayer membrane depicted as blue squares; the characters CP and EC mean cytoplasmic and extracellular sides of the cell, respectively. The 10 amino-acid residues are shown as space-filling models with labels and atom-type colorings. This figure was produced using the software Discovery Studio T A B L E 1 Ten selected candidates of thermo-stabilizing mutants Characters in parentheses indicate the location of residues where N and C represent the N-terminal and the C-terminal sides of the protein, respectively.
b Total folding-free-energy differences upon mutation from the "Build Model" function of FoldX.
c Suggestions from FoldX with their contributions to the folding freeenergy differences in parentheses.
United Kingdom). The signal of a buffer reference measured under identical conditions was subtracted from the signal of the protein samples, the data were normalized to the molar protein concentration, and an interpolated baseline was subtracted. The T m value was determined at the peak top of a heat capacity curve. For each construct, three independent experiments were performed to obtain the average T m value. To evaluate the statistical significance between a pair of average T m values from different constructs, a two-tailed t-test on averages was applied against the null hypothesis where the averages were the same; for each construct, the accordance with a normal distribution was confirmed at the significance level of 5% using the Kolmogorov-Smirnov test. The equality of valiances was examined at the significance level of 5% using the F-test. The Student's t test and the Welch's t test were adopted in the cases of equal and unequal variances, respectively.

| Molecular dynamics calculations
The MD simulation of TR was performed using the program GROMACS 18 v5.0.5, with the CHARMM36 force field. 19 Unless otherwise noted, the parameter was set to the default value in GROMACS.
The chain A of 5AZD was used to prepare the protein model for the wild-type TR. To prepare a mutant model, the residue to be mutated was replaced by a designated residue type with a rotamer conformation selected so as to provide favorable intramolecular interactions For the MD simulation of the TR system in the native state, first the system was equilibrated in an NPT ensemble at 300 K using GROMACS, according to the CHARMM-GUI procedure: a steepestdescents energy minimization with a tolerance of 1000 kJ /mol/nm; the first equilibration for 25 ps at 300 K with a random seed of initial velocity, a time step of 1 fs, the LINCS algorithm to constrain bonds involving hydrogen atoms, 22 the particle-mesh Ewald (PME) method for the electrostatic calculation, 23 and with a temperature coupling according to Berendsen 24 ; the second equilibration for 25 ps in the same way but with decreased geometric restraints; the third equilibration for 25 ps in the same way but with a pressure coupling at 1.0 bar according to Berendsen 24 and with decreased restraints; the fourth equilibration for 100 ps in the same way but with a time step of 2 fs and with decreased restraints; the fifth/sixth equilibration for 100 ps in the same way but with decreased restraints. Then, an NPT simulation at 300 K for 100 ns was performed according to the CHARMM-GUI protocol adopting a time step of 2 fs, the LINCS constraints, the PME electrostatics, the Nosé-Hoover thermostat, 25 and the Parrinello-Rahman barostat 26 at 1.0 bar; no geometric restraints was applied, the number of steps between the neighboring energy calculations (ie, nstcalcenergy) was set as 1 000 steps, the linear mode was selected as the center of mass motion removal at every 1 000 steps.
The calculation performance was 2.3 to 6.2 ns per day using a Linux PC with 11 to 30 threads. The native-state simulation was repeated for three or four times. For the MD simulation in the random-coil state, the TR system equilibrated at 300 K in the same way as that in the native-state simulation was submitted to an NVT simulation at 1 800 K for 1 ns with a time step of 1 fs, using the Berendsen thermostat, the LINCS constraints, and the PME electrostatics; no geometric restraints was applied, the nstcalcenergy was set as 200 steps; the linear mode was selected as the center of mass motion removal at every 200 steps. Since this high-temperature simulation could be executed within a half day using a Linux PC with 11 to 30 threads, 20 independent runs were performed for the statistical evaluation. The NPT equilibration at 300 K was performed with different initial velocity before each NVT run at 1 800 K, to avoid a complete reproduction.  3 | RESULTS

| Selection of candidates for thermo-stabilizing mutants
From 28 solvent-exposed residues in the extramembrane regions, we designed ten thermo-stabilizing mutants of TR based on the folding free-energy calculation using the program FoldX ( Figure 1 and Table 1). Of these single mutations, six and four sites were on the N-and C-terminal side of the protein, respectively. The that located at the N-termini of helices, the same stabilization in the mutants A177K and W210R that located at the C-termini of helices, and a stabilization due to the main-chain entropy in the mutants A115P and A215P that located at the N-termini of helices. The other suggestions were stabilization in the V79K mutant due to an ion-pair formation between the side chains of Glu6 and Lys79, and a stabilization in the T143K mutant due to a hydrogen-bond formation between the main-chain oxygen of Val141 and the ε-amino group of Lys143.

| Thermostability of TR
We purified the wild-type and 10 single mutants (S8D, V79K, T114D, A115P, A116E, T143K, A177K, W210R, G214D, and A215P) of TR. The UV-visible spectroscopy of these purified samples exhibited the maximum absorption at the same wavelength, 530 nm. It is generally accepted for the rhodopsin family that the maximum absorption is strongly correlated with the protein properties. [12][13][14] Therefore, to compare the thermostability of wild-type and the 10 mutants of TR, a thermal-decoloration of pigment at 530 nm after incubation at 90 C for 4 minutes was quantified ( Table 2). As a result, the residual pigment values of four mutants (V79K, T114D, A115P, and A116E) were significantly higher than that of the wild type, suggesting a mutationinduced thermo-stabilization, whereas the other six mutants were suggested to be destabilized. The four potentially stabilized mutants were characterized further. A kinetics study using the UV-visible spectroscopy at 530 nm showed that the thermal decoloration was significantly slower in the four mutants when compared to the wild type ( Figure 2A and Table 2). Another characterization by DSC showed the higher T m values of the four mutants when compared to the wild type (P < .01; Figure 2B and Table 2), although these are apparent values because the denaturation of TR was an irreversible/ kinetic process in the DSC experiment (Supplementary Figure S1). Residual pigment is calculated as "the maximum absorption at the wavelength 530 nm of TR after 4-minutes incubation at 90 C" divided by "that before incubation" and is represented in the percentage; the higher value indicates a thermo-stabilization. b For a TR construct, the relative amount of purified sample when compared to that of the wild type is shown. The expression level of T143K was too low to purify it for further characterizations. The difference of the residual pigment between T143K and WT was not also significant, thereby regarding it as a destabilizing mutant.  To evaluate the compactness of the TR molecule in the MD simulation at 1 800 K, the radius of gyration was calculated for the wildtype TR ( Figure 3D). When compared to the radius of gyration of the crystal structure of about 2.0 nm, that of the MD structure after 1 ns increased by 1.5-fold to be about 3.0 nm. This coincides with a reported experimental transition in the radius of gyration between the native and the random-coil states of apomyoglobin from a small-angle X-ray scattering study, 29 whereas, this MD structure is clearly different from that of a partially unfolded "molten-globule state" 30 with retaining secondary structures. Thus, the MD simulation of TR for 1 ns at 1800 K may provide a random-coil structure. Importantly, this random-coil TR molecule is highly reproducible in the simulation, and it does not be extended completely, indicating that long-range intramolecular interactions still present in this condition.

T A B L E 2 Experimental characterization of TR
According to the Boltzmann's principle, the entropy of a system increases by increasing the number of microscopic states of the system, indicating that an increment in the intramolecular structural fluc-

| DISCUSSION
Previously, we reported a hyperthermostabilization design of a soluble globular protein, CutA1, which was based on the computational prediction of folding free-energy differences upon mutation. 11 In the present study, we applied the same strategy to a thermostable membrane protein, TR. Here we designed 10 mutants in the extramembrane regions of TR, and four of them were found to be stabilized, indicating a hit-rate of 40%. This value is comparable to 39% in the case of CutA1. 31 Therefore, introducing a single mutation in the extramembrane regions based on the folding free-energy prediction may be useful to design a thermo-stabilized TR mutant. In the next step, the four successful single mutations obtained here will be used to design a further thermo-stabilized TR mutant with multiple mutations. It should be noted that the predicted changes in free energy for Model" function of FoldX may overestimate the free-energy changes.
Other MD-based design methods such as the free-energy perturbation calculation 32 may be examined in future to compare results. For practical applications such as optogenetics tools and biosensors, the functional activity of the mutants should be evaluated carefully.
In this study, we examined the possibility of an improvement of the thermo-stabilization design using the MD simulation. Although the MD simulation in the native state could not improve the hit rate of the design, it suggested a mutation-induced structural change of the T114D mutant that revealed the highest T m in the DSC experiment. Thus, the T114D mutation is thermo-stabilizing but may induce a structural change in the native state, although it is unclear whether the mutation-induced structural change contributes to the thermostabilization or not. This result indicates that the MD simulation in the native state may suggest a mutation-induced structural change that can affect the protein function. On the other hand, the MD simulation in the random-coil state could improve the hit rate by detecting an entropic destabilization of the S8D mutant that revealed the lower thermostability than that of the wild type. Thus, in certain cases, the MD simulation may suggest candidates with problems in the denatured state that is not considered explicitly in the folding free-energy prediction.
In conclusion, the present thermo-stabilization design in the extramembrane regions based on the free-energy calculation and the subsequent evaluation by the MD simulation may be useful to improve the production of membrane proteins. The engineering in the intramembrane region. 7,8 would be a complementary design method.
These two strategies in the extra-and intramembrane regions will be combined in future to design hyperthermostable mutants of TR tolerable at temperatures over 100 C. Further applications of the strategies to other membrane proteins may expand the availability of the rational design of proteins in general.

ACKNOWLEDGMENTS
The authors thank Dr. Hisashi Naitow contributed to the computational prediction. Takeshi Murata planned and supervised this study and co-wrote the manuscript.