If the distribution of sample means is skewed, can I infer that the samples are not from identical populations? My understanding of the Central Limit Theorem is that the distribution of the sample mean will tend towards normal, for random sampling with replacement.
How reasonable is it for me to infer that if I have a bunch of samples but don't know if the underlying population is the same, a high skewness value is (perhaps circumstantial) evidence that the underlying population is in fact different (or the sampling wasn't random, etc)?
In my specific case, I am dealing with a problem where I have data samples from around 1000 different entities regarding compliance. Say it is wastewater samples with various amounts of lead, and there are ~1500 samples per entity. I have plotted the mean lead levels from each entity, and the distribution is noticeably right-skewed (R's 'moments' library gives a skewness of 0.36). The entities will claim that any high sample means are just due to the vagaries of the lead-pollution business, but does the skewness/lack of normality of data constitute (perhaps weak or circumstantial) evidence that in fact the underlying population is not identical, and some locations are systematically (rather than randomly) worse than others? My intuition was that if the lead level samples were normally distributed, it would fit with the claim that some entities are just unlucky (and in practice, be an argument against investing extra resources in investigating those entities compared with the other ones). But I feel like intuition is not a very reliable guide in the jungles of statistical parameter interpretation!
 A: NO
The central limit theorem is a limit theorem. In the limit, as $n\rightarrow\infty$, there is a convergence (and that convergence isn’t even quite of the sample mean). In finite samples, you don’t get the convergence, and some distributions have slow convergence.
For instance, a log-normal can produce a rather skewed distribution of sample means.
set.seed(2022)
N <- 1000
R <- 10000
xbars <- rep(NA, R)
for (i in 1:R){

    x <- rlnorm(N, 0, 2.5)
    xbars[i] <- mean(x)
}

hist(xbars, 25)
boxplot(xbars)



These log-normal distributions are $iid$ with finite variance, thus meeting the conditions of the central limit theorem, yet the distribution of sample means with a fairly large sample size of $1000$ is skewed.
Finally, note that the central limit theorem concerns the convergence of a transformation of the sample mean, not the sample mean itself. The distribution of the sample mean itself need not have support on the whole real line, such as this log-normal example that cannot take negative values and thus cannot have a sample mean less than zero.
$$
\text{Classical Central Limit Theorem}\\
\sqrt{n}\dfrac{
\bar X_n -\mu
}{
\sigma
}
\overset{d}{\rightarrow}N(0,1)\\
\bar X_n=\dfrac{X_1+\dots +X_n}{n}
$$
