Let $M$ be a smooth manifold. Let $g$ be a Riemannian metric of regularity $C^k(M)$; that is, in any coordinates, $g$'s components' coordinate representations are $C^k$ functions in that chart. When one speaks of "coordinates" on a smooth manifold $M$, they are implicitly assumed to belong to the maximal smooth atlas $\mathscr{A}$ implicit in the definition of $M$.
Harmonic "coordinates" can be found using standard elliptic theory: starting from a chart $(U, x)$ on $M$ in which $g$ has $C^k$ components, one obtains $y : U \to \mathbb{R}^n$ solving $\Delta_g y^j = 0$ for $j = 1, \ldots, n$ ($n = \text{dim} M$). However, the resulting theory gives that the coordinates are $C^{k + 1}$ functions of $x$, and therefore are not guaranteed to be smooth local coordinates - hence not belong to the maximal atlas $\mathscr{A}$!
So in what sense are these coordinates? Does one have to relax the smooth structure of the manifold $M$?
You are right that harmonic coordinates of a $C^k$-metric are of class $C^{k+1}$ and thus, strictly speaking, do not belong to the underlined smooth atlas. But this is not a big deal: they are still $C^{k+1}$-compatible with any other coordinates of the said smooth atlas. They are "coordinates" in the sense that they give local parametrizations of the manifold, but only with $C^{k+1}$ maps. One can construct a maximal $C^{k+1}$-atlas out of the underlined smooth atlas simply by considering all the $C^{k+1}$ coordinates / charts that are compatible with the smooth atlas: harmonic coordinates belong to this $C^{k+1}$-atlas.
That said, harmonic coordinates are very nice in some contexts: since many geometric objects take a nice form in harmonic coordinates, they are often suitable for solving PDEs or computing the regularity of some other function. Afterwards, the result holds in the smooth atlas (with at least the same regularity).