We study the $L^q$-regularity of the density of barycenters of $N$ probability measures on $\mathbb{R}^d$ with respect to the $p$-Wasserstein metric ($1<p<\infty$). According to a previous result by the first author and collaborators, if one marginal is absolutely continuous, so is the $W_p$-barycenter. The next natural question is whether the $L^q$- regularity on the marginals is also preserved for any $q > 1$, as in the classical case ($p=2$) of Agueh--Carlier, or for $W_p$-geodesics ($N=2$). Here we prove that this is the case if one marginal belongs to $L^q$ and the supports of all the marginals satisfy suitable geometric assumptions. However, we show that, as soon as $N>2$, it is possible to find examples of $W_p$-barycenters which are not $q$-integrable, even if one marginal is compactly supported and bounded, thus highlighting the role played by the geometry of the supports. Furthermore, we provide a general estimate of the $L^q$-norm, including a detailed study of the sources of singularities, and a characterization of the $W_p$-barycenters à la Agueh--Carlier in terms of the associated Kantorovich potentials. Finally, we explicitly compute the $W_p$-barycenters of measures obtained as push-forward of special affine transformations. In this case, regularity holds without any additional requirement on the supports.