This study presents a multi-output Gaussian process regression (GPR) surrogate model for seismic-fragility analysis of bridge structural systems. Multi-output GPR can model the correlations among multiple outputs, accuracy and stability are achieved with fewer training data, which reduces the computational cost of fragility analysis. Furthermore, the explainability of the constructed surrogate model is implemented by adopting an automatic relevance-determination (ARD) kernel in the GPR. The estimated hyperparameters can provide the contribution of the uncertainty of each input parameter to the outputs. The fragility analysis using the multi-output GPR surrogate model was verified by applying it to a seismic isolation highway bridge with multiple spans and a curved geometry. The effectiveness of the multi-output GPR was demonstrated by the construction of an accurate and stable surrogate model with 46 inputs and 28 outputs. The relative contributions of the uncertainties to the structural properties and input earthquake loads could also be understood. The fragility curves, at both the component and system levels, were appropriately obtained using a sufficient number of samples in a Monte Carlo calculation. Furthermore, the failure modes were evaluated, identifying which structural components contributed to the system failure. This enabled discussions on structural system failures from the viewpoint of the structural dynamic characteristics of the bridge and earthquake-load properties.