A kernel density estimator for data on the polysphere $\mathbb{S}^{d_1}\times\cdots\times\mathbb{S}^{d_r}$, with $r,d_1,\ldots,d_r\geq 1$, is presented in this paper. We derive the main asymptotic properties of the estimator, including mean square error, normality, and optimal bandwidths. We address the kernel theory of the estimator beyond the von Mises-Fisher kernel, introducing new kernels that are more efficient and investigating normalizing constants, moments, and sampling methods thereof. Plug-in and cross-validated bandwidth selectors are also obtained. As a spin-off of the kernel density estimator, we propose a nonparametric $k$-sample test based on the Jensen-Shannon divergence. Numerical experiments illuminate the asymptotic theory of the kernel density estimator and demonstrate the superior performance of the $k$-sample test with respect to parametric alternatives in certain scenarios. Our smoothing methodology is applied to the analysis of the morphology of a sample of hippocampi of infants embedded on the high-dimensional polysphere $(\mathbb{S}^2)^{168}$ via skeletal representations ($s$-reps).