The zearalenone hydrolase (ZHD101) derived from Clonostachys rosea can effectively degrade the mycotoxin zearalenone (ZEN) present in grain by-products and feed. However, the low thermal stability of ZHD101 hampers its applications. High throughput screening of variants using spectrophotometer is challenging because the reaction of hydrolyzing ZEN does not change absorbance. In this study, we used ZHD101 as a model enzyme to perform computation-aided design followed by experimental verification. By comparing the molecular dynamics simulation trajectories of ZHD101 at different temperatures, 32 flexible sites were selected. 608 saturated mutations were introduced into the 32 flexible sites virtually, from which 12 virtual mutants were screened according to the position specific score and enzyme conformation free energy calculation. Three of the mutants N156F, S194T and T259F showed an increase in thermal melting temperature (ΔTm>4 °C), and their enzyme activities were similar to or even higher than that of the wild type (relative enzyme activity 95.8%, 131.6% and 169.0%, respectively). Molecular dynamics simulation analysis showed that the possible mechanisms leading to the improved thermal stability were NH-π force, salt bridge rearrangement, and hole filling on the molecular surface. The three mutants were combined iteratively, and the combination of N156F/S194T showed the highest thermal stability (ΔTm=6.7 °C). This work demonstrated the feasibility of engineering the flexible region to improve enzyme performance by combining virtual computational mutations with experimental verification.