Machine learning (ML) is a rising and promising tool for Reynolds-Averaged Navier–Stokes (RANS) turbulence model developments, but its application to industrial flows is hindered by the lack of explainability of the ML model. In this paper, two types of methods to improve the explainability are presented, namely the intrinsic methods that reduce the model complexity and the post-hoc methods that explain the correlation between the model inputs and outputs. The investigated ML-assisted turbulence model framework aims to improve the prediction accuracy of the Spalart–Allmaras (SA) turbulence model in transonic bump flows. A random forest model is trained to construct a mapping between the input flow features and the output eddy viscosity difference. Results show that the intrinsic methods, including the hyperparameter study and the input feature selection, can reduce the model complexity at a limited cost of accuracy. The post-hoc Shapley additive explanations (SHAP) method not only provides a ranked list of input flow features based on their global significance, but also unveils the local causal link between the input flow features and the output eddy viscosity difference. Based on the SHAP analysis, the ML model is found to discover: (1) the well-known scaling between eddy viscosity and its source term, which was originally found from dimensional analysis; (2) the well-known rotation and shear effects on the eddy viscosity source term, which was explicitly written in the Reynolds stress transport equations; and (3) the pressure normal stress and normal shear stress effect on the eddy viscosity source term, which has not attracted much attention in previous research. The methods and the knowledge obtained from this work provide useful guidance for data-driven turbulence model developers, and they are transferable to future ML turbulence model developments.