Symbolic regression is to search the space of mathematical expressions to find a model that best fits a given dataset. As genetic programming (GP) with the tree representation can represent solutions as expression trees, it is popularly-used for regression. However, GP tends to evolve unnecessarily large programs (known as bloat), causing excessive use of CPU time/memory and evolving solutions with poor generalization ability. Moreover, even though the importance of local search has been proved in augmenting the search ability of GP (termed as memetic algorithms), local search is underused in GP-based methods. This work aims to handle the above problems simultaneously. To control bloat, a multi-objective (MO) technique (NSGA-II, Non-dominant Sorting Genetic Algorithm) is selected to incorporate with GP, forming a multi-objective GP (MOGP). Moreover, three mutation-based local search operators are designed and incorporated with MOGP respectively to form three multi-objective memetic algorithms (MOMA), i.e. MOMA_MR (MOMA with Mutation-based Random search), MOMA_MF (MOMA with Mutation-based Function search) and MOMA_MC (MOMA with Mutation-based Constant search). The proposed methods are tested on both benchmark functions and real-world applications, and are compared with both GP-based (i.e. GP and MOGP) and nonGP-based symbolic regression methods. Compared with GP-based methods, the proposed methods can reduce the risk of bloat with the evolved solutions significantly smaller than GP solutions, and the local search strategies introduced in the proposed methods can improve their search ability with the evolved solutions dominating MOGP solutions. In addition, among the three proposed methods, MOMA_MR performs best in RMSE for testing, yet it consumes more training time than others. Moreover, compared with six reference nonGP-based symbolic regression methods, MOMA_MR generally performs better than or similar to them consistently.