Biological evidence demonstrates that fish can tune their body stiffness to improve thrust and efficiency during swimming locomotion. However, the stiffness-tuning strategies that maximize swimming speed or efficiency are still unclear. In the present study, a musculo-skeletal model of anguilliform fish is developed to study the properties of variable stiffness, in which the planar serial-parallel mechanism is used to model the body structure. The calcium ion model is adopted to simulate muscular activities and generate muscle force. Further, the relations among the forward speed, the swimming efficiency, and Young’s modulus of the fish body are investigated. The results show that for certain body stiffness, the swimming speed and efficiency are increased with the tail-beat frequency until reaching the maximum value and then decreased. The peak speed and efficiency are also increased with the amplitude of muscle actuation. Anguilliform fish tend to vary their body stiffness to improve the swimming speed and efficiency at a high tail-beat frequency or small amplitude of muscle actuation. Furthermore, the midline motions of anguilliform fish are analyzed by the complex orthogonal decomposition (COD) method, and the discussions of fish motions associated with the variable body stiffness and the tail-beat frequency are also presented. Overall, the optimal swimming performance of anguilliform fish benefits from the matching relationships among the muscle actuation, the body stiffness, and the tail-beat frequency.