Assessing the seismic capacity of pier columns is a crucial element in the performance-based seismic design of bridges. Such assessment necessitates a probabilistic approach to accurately determine the marginal probability distributions of seismic capacities and to characterize the dependencies among these variables. In response to this need, this paper employs Multi-Output Gaussian Process Regression (MO-GPR), a probabilistic machine learning method, to jointly model the multiple seismic capacities of pier columns. We initially introduce a probabilistic seismic capacity model that utilizes MO-GPR for pier columns and validate its predictive accuracy in comparison to Bayesian linear regression and existing empirical methods. Subsequently, the methodology is augmented by the integration of hierarchical modeling within the MO-GPR framework, resulting in a Multi-Output Hierarchical Gaussian Process Regression (MO-HGPR) model that effectively estimates intraclass correlation coefficients for specific types of datasets. It is postulated that these correlation coefficients also reflect correlations associated with multiple components of the real structure. This study employs MO-HGPR and MO-GPR separately to investigate the potential correlations of seismic capacities among pier columns and distinct limit states. The results demonstrate that the MO-GPR model exhibits superior prediction accuracy and more effectively portrays uncertainty compared to existing empirical models. Moreover, the correlations of seismic capacities among piers and limit states are both robust and significantly impact the seismic fragility of bridges. This finding highlights the essential nature of considering capacities correlations in seismic fragility or risk assessment processes.