One place where (a generalization of) the desired result is stated explicitly is in a 1973 paper by J. H. Sampson, On a theorem of Chern. Sampson gives a simplified proof of Chern's result in his 1957 paper On a generalization of Kähler geometry (Algebraic geometry and topology. A symposium in honor of S. Lefschetz, pp. 103–121. Princeton University Press, Princeton, N. J., 1957).
The form in which Sampson states it is as follows:
Theorem: If a Riemannian manifold $(M^n,g)$ has holonomy $H\subseteq\mathrm{O}(n)$, then the $g$-Laplacian commutes with all of the linear operators on $\Omega^*(M)$ constructed from the ring of $H$-equivariant linear maps $L:\Lambda^*(\mathbb{R}^n)\to\Lambda^*(\mathbb{R}^n)$.
(Chern had proved this in his 1957 paper for the ring of operators $L$ as above that are degree preserving, but, in fact, this implies the more general result. To be more precise, both Chern and Sampson phrase the result in terms of $H$-structures without torsion, but, of course, this is the same as Riemannian manifolds with holonomy contained in $H$.)
The desired result is a special case of this, since, if $\pi$ is a $g$-parallel $p$-form, then the holonomy $H$ of $g$ is contained in the stabilizer of $\pi$, and hence, by the above result, the operator $L(\alpha) = \alpha\wedge\pi$ commutes with the Laplacian of $g$. Of course, this implies the desired result since $\Delta(\alpha\wedge\pi) = \Delta\bigl(L(\alpha)\bigr) = L\bigl(\Delta(\alpha)\bigr) = 0$ when $\Delta(\alpha) = 0$, but is stronger.
There are earlier results in special cases by André Lichnerowicz, e.g. Généralisations de la géométrie kählérienne globale. (French) Colloque de Géométrie Différentielle, Louvain, 1951, pp. 99–122. Georges Thone, Liège; Masson & Cie., Paris, 1951, but neither his articles nor those of Chern or Sampson mention Bott.
There is also a very nice 1962 exposition (in French) of Chern's theorem by André Weil, Un théorème fondamental de Chern en géométrie riemannienne, in Séminaire N. Bourbaki, 1962, exp. no 239, p. 273-284.