Shear strength prediction for reinforced concrete (RC) beams is a complex non-linear problem impacted by many factors. Traditional Machine Learning (ML) models for shear strength prediction still need manual single-model selection and hyperparameter optimization, limiting the model deployment efficiency and prediction accuracy in practical engineering. This study proposes a novel shear strength prediction model based on Bayesian Optimization and multi-layer stacking (BO-Stacking). It trains and stacks several base ML models in a multi-layer structure to avoid manual model selection and boost shear strength prediction accuracy concurrently. Besides, Bayesian Optimization is developed to automatically find optimal hyperparameters of each base ML model, substantially reducing computational cost during training. With the developed techniques, the proposed model is able to automatically achieve high accuracy by simply specifying training time. Furthermore, SHapley Additive exPlanation is adopted to interpret the feature contributions in shear strength prediction. The proposed model is implemented on an experimental dataset consisting of 406 RC beams with stirrups and then compared with seven base ML models and three shear design codes. Comparative results show that BO-Stacking achieves the highest accuracy on the test set with R², RMSE and MAE of 0.9551, 39.98 and 25.99 respectively. The reliability and robustness of BO-Stacking is further verified through the trend analysis and Monte Carlo simulation. The proposed BO-Stacking model automates the tedious tasks of implementing an ML model for shear strength prediction and achieves stable and accurate results, facilitating the efficiency and reliability of model deployment in practice.