To address the challenges of poor model interpretability and the limited prediction accuracy caused by strong multi-load coupling in integrated energy systems, an optimized XGBoost model integrating the Multi-Input Multi-Output (MIMO) strategy and the SHapley Additive exPlanations (SHAP) framework is proposed to achieve the joint forecasting of electricity, cooling, and heating loads. The results show that the mean absolute percentage errors for electricity, cooling, and heating load forecasting are 2.75%, 3.45%, and 3.44%, respectively, representing reductions of 12.8%, 15.3%, and 14.7% compared to the single-task XGBoost model. The global and local SHAP analysis quantifies the impact of input features, with a feature contribution ranking consistency reaching 0.8462, explicitly identifying temperature, historical loads, and solar radiation as the key factors affecting prediction accuracy. It is concluded that the proposed method effectively and explicitly captures the physical constraints and dynamic correlations among loads, significantly improving both the multi-load forecasting accuracy and the decision transparency of the model for integrated energy systems.