The problem is equivalent to finding $\min_{e_i, \lambda} \sum (w(e_i) - \lambda)^2$, where $\lambda$ is a free real parameter subject to optimization as well as the edges of a spanning tree. Indeed, for a chosen set of edges $e_i$ we have that $\sum (w(e_i) - \lambda)^2$ is minimized when $\lambda$ is the mean of $w(e_i)$.
If $\lambda$ is fixed, we can use, say, Kruskal's algorithm which considers edges by increasing of $|w(e) - \lambda|$. Note that if we gradually increase $\lambda$, the order of edges only changes at $\lambda = (w(e) + w(e')) / 2$ for some edges $e, e'$, therefore there are polynomially many different orders. Run Kruskal's algorithm with respect to each of these orders and choose the best answer.
The complexity is roughly $O(|E|^3)$ if done straightforwardly. UPD: This approach can be optimized to $O(|E|^2 \log |V|)$, roughly as follows: when we process a "two edges $e_1$, $e_2$ switch their order to $e_2$, $e_1$" event, the result of the Kruskal algorithm changes iff all edges on the cycle of the current MST induced by $e_2$ are $\geq e_1$, then we have to delete $e_1$ and insert $e_2$. This can be done efficiently with link-cut trees.