We use pseudospectral direct numerical simulations to solve the three-dimensional (3D) Hall-Vinen-Bekharevich-Khalatnikov (HVBK) model of superfluid helium. We then explore the statistical properties of inertial particles, in both coflow and counterflow superfluid turbulence (ST) in the 3D HVBK system; particle motion is governed by a generalization of the Maxey-Riley-Gatignol equations. We first characterize the anisotropy of counterflow ST by showing that there exist large vortical columns. The light particles show confined motion as they are attracted toward these columns, and they form large clusters; by contrast, heavy particles are expelled from these vortical regions. We characterize the statistics of such inertial particles in 3D HVBK ST: (1) The mean angle Θ (τ) between particle positions, separated by the time lag τ, exhibits two different scaling regions in (a) dissipation and (b) inertial ranges, for different values of the parameters in our model; in particular, the value of Θ (τ), at large τ, depends on the magnitude of U n s. (2) The irreversibility of 3D HVBK turbulence is quantified by computing the statistics of energy increments for inertial particles. (3) The probability distribution function (PDF) of energy increments is of direct relevance to recent experimental studies of irreversibility in superfluid turbulence; we find, in agreement with these experiments, that, for counterflow ST, the skewness of this PDF is less pronounced than its counterparts for coflow ST or for classical fluid turbulence. © 2023 Author(s).