We have developed a methodology for identifying further thermostabilizing mutations for an intrinsically thermostable membrane protein. The methodology comprises the following steps: (1) identifying thermostabilizing single mutations (TSSMs) for residues in the transmembrane region using our physics-based method; (2) identifying TSSMs for residues in the extracellular and intracellular regions, which are in aqueous environment, using an empirical force field FoldX; and (3) combining the TSSMs identified in steps (1) and (2) to construct multiple mutations. The methodology is illustrated for thermophilic rhodopsin whose apparent midpoint temperature of thermal denaturation Tm is ∼91.8 °C. The TSSMs previously identified in step (1) were F90K, F90R, and Y91I with ΔTm ∼5.6, ∼5.5, and ∼2.9 °C, respectively, and those in step (2) were V79K, T114D, A115P, and A116E with ΔTm ∼2.7, ∼4.2, ∼2.6, and ∼2.3 °C, respectively (ΔTm denotes the increase in Tm). In this study, we construct triple and quadruple mutants, F90K+Y91I+T114D and F90K+Y91I+V79K+T114D. The values of ΔTm for these multiple mutants are ∼11.4 and ∼13.5 °C, respectively. Tm of the quadruple mutant (∼105.3 °C) establishes a new record in a class of outward proton pumping rhodopsins. It is higher than Tm of Rubrobacter xylanophilus rhodopsin (∼100.8 °C) that was the most thermostable in the class before this study.