Consider a sample of size $n$ with values $(y_1, y_2, \ldots, y_n)$ drawn independently and randomly. The sample mean is
$$\bar y = \frac{1}{n}\sum_i y_i$$
and the sample variance is
$$s^2_y = \frac{1}{n-1} \left(\sum_i y_i^2 - \frac{1}{n}\left[\sum_j y_j\right]^2\right)$$
where all summation indexes range from $1$ through $n.$
Consequently, $\bar y s^2_y$ is a homogeneous cubic form in the $y_i.$ This means it is a linear combination of terms of the form $y_i^3,$ $y_iy_j^2,$ and $y_iy_jy_k$ where $i\ne j\ne k.$ Because the $y_i$ are identically distributed, the actual index values don't matter: we merely have to count how many of each kind of term there is and multiply their (common) coefficients by their counts.
Write $\mu_k$ for the raw $k^\text{th}$ moment $\mu_k = E[Y^k]$ where $Y$ is any random variable with this common distribution. Plug these formulas for $\bar y$ and $s^2_y$ into the expectations to find the coefficients and do the counting to obtain
$$E[\bar y s^2_y] = \frac{1}{n}(\mu_3 - (n-2)\mu_1^3 + (n-3)\mu_1\mu_2)$$
and
$$E[\bar y] = \mu_1;\ E[s^2_y] = \mu_2 - \mu_1^2$$
(as is well known).
For example, the coefficient of $y_1^3$ in $\bar y s^2_y$ is $1/n$ (from the coefficient of $y_1$ in $\bar y$) multiplied by $1/(n-1)$ times $1 - 1/n$ (from the coefficient of $y_1^2$ in $s^2_y$). Because there are $n$ such terms $y_1^3, y_2^3, \ldots, y_n^3,$ the coefficient of $\mu_3$ must be
$$n\left(\frac{1}{n}\right)\left(\frac{1}{n-1}\right) \left(1 - \frac{1}{n}\right) = \frac{1}{n}.$$
The other two coefficients in $E[\bar y s^2_y]$ are found similarly.
Thus, applying a tiny bit more algebra, we find
$$\operatorname{Cov}(\bar y, s^2_y) = E[\bar y s^2_y] - E[\bar y]E[s^2_y] = \frac{1}{n}(\mu_3 + 2\mu_1^3 - 3 \mu_1\mu_2).$$
In particular, in any symmetric distribution $\mu_1=\mu_3=0$ and the covariance must then be zero (but this doesn't imply $\bar y$ and $s^2_y$ are independent unless the distribution is Normal).
It is so easy to make little algebraic mistakes that an independent check is worthwhile. Here is R
code to draw repeated samples from any (finite) population and compute their means and variances. At the end it prints the covariance of the means and the variances followed by the value given by this formula. Here is an example:
Sample Formula
0.06382658 0.06380878
They agree closely, with small differences attributed to the randomness of the simulation.
#
# Create any finite population.
#
pop <- rexp(30)
#
# Create a sampling distribution of sample means and sample variances.
#
n <- 10 # Sample size: must be 2 or greater
sim <- replicate(5e4, {
y <- sample(pop, n, replace = TRUE)
c(mean(y), var(y))
})
#
# Compare the covariance of the simulation to the formula.
#
mu <- function(x, k = 1) mean(x^k) # Raw moments of `x`
v <- (mu(pop, 3) + 2 * mu(pop, 1)^3 - 3 * mu(pop, 1) * mu(pop, 2)) / n
c(Sample = cov(sim[1, ], sim[2, ]), Formula = v)