We introduce the Wick integral $\int_s^t p(X_u) \Diamond \mathrm{d} X_u$ for a class of stochastic processes $X$ which are not necessarily Gaussian, in the regime of bounded $2> q$-variation. The integral is defined for polynomial integrands, and has the property of being centred if $X$ is such. In the case of $1/2 < H$-fractional Brownian motion, the Wick integral agrees with the divergence operator in Malliavin calculus. It satisfies a correction formula with the Young integral $\int p(X)\mathrm{d} X$ and an Itô formula which have arbitrarily many correction terms (only limited by the degree of $p$), given by integration against the cumulant functions of $X$, and reduce to familiar identities in the Gaussian case. These results are obtained by first developing diagram formulae for Appell polynomials. Our theory applies to a range of processes taking values in bounded Wiener chaos, such as the Rosenblatt process.