Overview
As in the question, suppose data points are sampled i.i.d. from a Gaussian distribution, but we haven't observed these values. Instead, we've observed a set of sample means for different subsets of the data points. I'll assume these subsets don't intersect, but you could relax this assumption and follow a similar approach. Below, I'll find maximum likelihood estimators for the mean and variance of the distribution, given the observed sample means.
As one might expect, the maximum likelihood estimator for the mean is a weighted mean of the observed sample means, with weights proportional to the number of data points contributing to each observation. But, the maximum likelihood estimator for the variance is not the analogously weighted variance of the observations (as suggested in the question). Rather, it's a scaled version of this. Scaling compensates for the fact that averaging reduces the variability of the observed sample means compared to the original data.
The averaging used to produce the observations discards information about the original data points. So, even if there are many data points, estimating the variance as described here can't be expected to work well if there are too few observed sample means.
Formulating the problem
Let $X_1, \dots, X_n$ be random variables representing the original $n$ data points, which are assumed i.i.d. Gaussian with mean $\mu$ and variance $\sigma^2$. We can represent the entire dataset with a random vector $X = [X_1, \dots, X_n]^T$, which has a multivariate Gaussian distribution:
$$p(x \mid \mu, \sigma^2) = \mathcal{N}(x \mid \mu \vec{1}_n, \sigma^2 I_n)$$
where $\vec{1}_n$ is a $n \times 1$ vector of ones and $I_n$ is the $n \times n$ identity matrix.
Unfortunately, we haven't directly observed the data points, so we treat them as latent variables. Instead, we have $k \le n$ observations represented by the random vector $M = [M_1, \dots, M_k]^T$. Each observation $M_i$ is the sample mean of a specified subset of data points. This subset contains $n_i$ data points, whose indices are given in set $S_i$:
$$M_i = \frac{1}{n_i} \sum_{j \in S_i} x_j$$
This can be written more succinctly as:
$$M = W X$$
$W$ is a $k \times n$ weight matrix, where $W_{ij} = \frac{1}{n_i}$ if $j \in S_i$ (i.e. point $j$ contributes to observation $i$), otherwise $0$. This implies that each row of $W$ sums to one. I will also assume that each data point contributes to exactly one observation (so none of the $S_i$ intersect, and their union is $\{1, \dots, n\}$). This implies that each column of $W$ has exactly one nonzero entry.
The goal is to find values of the parameters $\mu, \sigma^2$ that maximize the marginal likelihood of the observed measurements. This is equivalent to minimizing the negative log likelihood, which will be more convenient:
$$\min_{\mu, \sigma^2} \ -\log p(m \mid \mu, \sigma^2)$$
Note that $M$ is a linear transformation of $X$. Because $X$ is Gaussian, $M$ is also Gaussian with convenient forms for the mean and covariance matrix:
$$p(m \mid \mu, \sigma^2) = \mathcal{N}(m \mid \mu W \vec{1}_n, \sigma^2 W W^T)$$
Furthermore, because of the special structure of $W$ above:
$$W \vec{1}_n = \vec{1}_k \quad
W W^T = \text{diag}(n_1^{-1}, \dots, n_k^{-1})$$
Therefore, the problem is:
$$\min_{\mu, \sigma^2} \ -\log \mathcal{N} \Big( m \left| \
\mu \vec{1}_k, \
\sigma^2 \text{diag}(n_1^{-1}, \dots, n_k^{-1})
\right. \Big)$$
Solution
Plug the expression for a Gaussian density function into the optimization problem above. After a little algebra to simplify things (taking advantage of the diagonal structure of the covariance matrix), this yields:
$$\min_{\mu, \sigma^2} \
\frac{k}{2} \log(2 \pi)
+ \frac{k}{2} \log \sigma^2
- \frac{1}{2} \sum_{i=1}^k \log n_i
+ \frac{1}{2 \sigma^2} \sum_{i=1}^k n_i (m_i - \mu)^2$$
Terms that don't depend on $\mu$ or $\sigma^2$ don't affect the solution and can be dropped:
$$\min_{\mu, \sigma^2} \
\frac{k}{2} \log \sigma^2
+ \frac{1}{2 \sigma^2} \sum_{i=1}^k n_i (m_i - \mu)^2$$
Differentiate the objective function with respect to $\mu$, set the derivative equal to zero, and solve. This yields the maximum likelihood estimate for the mean:
$$\mu_{ML} = \frac{1}{n} \sum_{i=1}^k n_i m_i$$
Notice that $\mu_{ML}$ is a weighted mean of the observations, where the weights are proportional to the number of data points contributing to each observation (since the $n_i$ sum to $n$ under our assumptions above).
Similarly, differentiate the objective function with respect to $\sigma^2$, set the derivative equal to zero, and solve. This yields the maximum likelihood estimate for the variance:
$$\sigma^2_{ML} = \frac{1}{k} \sum_{i=1}^k n_i (m_i - \mu_{ML})^2$$
This is not equivalent to the weighted variance of the observations (with weights proportional to $n_i$). In that case, we'd have:
$$\sigma^2_{W} = \frac{1}{n} \sum_{i=1}^k n_i (m_i - \mu_{ML})^2$$
But, notice that $\sigma^2_{ML} = \frac{n}{k} \sigma^2_W$. So, the maximum likelihood estimate of the variance is a scaled version of the weighted variance. The scaling factor will be greater than one if there are fewer observations than data points. This makes sense because the averaging used to produce the observations reduces their variability compared to the original data points. You could think of the scaling as compensating for this effect.
If the true mean is known (as stated in the question), it can be substituted into the expression for $\sigma^2_{ML}$ in place of $\mu_{ML}$.
