We propose an effectively nonparametric approach to fitting hidden Markov models to time series of counts, where the state-dependent distributions are estimated in a completely data-driven way without the need to specify a parametric family of distributions. To avoid overfitting, a roughness penalty based on higher-order differences between adjacent count probabilities is added to the likelihood, which is demonstrated to produce smooth state-dependent probability mass functions. The feasibility of the suggested approach is assessed in simulation experiments, and further illustrated in two real-data applications, where we model the distributions of (i) major earthquake counts and (ii) acceleration counts of an oceanic whitetip shark (Carcharhinus longimanus) over time. The proposed methodology is implemented in the accompanying R package countHMM, which is available on CRAN.