I think this has to do with what you treat as the dynamic variables. In your first formulation, $\rho$ was given and fixed, and you wanted to calculate $V$ by varying only $V$, not $\rho$. In the second formulation, $V$ is just a quantity derived from $\rho$, and this is the energy suitable for obtaining $\rho$ by varying $\rho$, with $V$ just a shorthand for a certain linear transformation of $\rho$. A third possibility is to consider $V$ as given and find the motion of a charge distribution in this potential -- here, again, there would not be a factor of $\frac{1}{2}$. This is the case, for instance, if we calculate the motion of the electron in the potential of a hydrogen nucleus, or of the Earth in the gravitational potential of the Sun.
Mathematically speaking, we need a factor of $\frac{1}{2}$ in front of quadratic terms (where $V\rho$ is quadratic in $\rho$ if $V$ is considered as a derived quantity derived from the dynamical variable $\rho$) but not in front of linear terms. Physically speaking, the factor of $\frac{1}{2}$ avoids double-counting when the energy is regarded as the interaction energy of a charge distribution with itself.
Now you may ask: But surely there is some well-defined energy and we can't just pick and choose the factors according to expedience? That's true, but consider again the hydrogen atom: If we consider the potential of the nucleus as given, then we're not counting the potential energy of the nucleus in the field of the electron in the integral over $\rho V$. But the entire energy is there and has to be accounted for in the integral, hence the factor $1$. On the other hand, in treating a helium atom, where both electrons are treated as dynamic and their interaction energy enters into the Hamiltonian, we do include a factor of $\frac{1}{2}$ to avoid double-counting the energy in the double integral over $\rho_1\rho_2/r_{12}$.
Best Answer
We will show that one need not calculate the energy stored in the electrostatic field by integration. To that end, we proceed.
We start with the electrostatic energy $W$ stored as given by
$$W=\frac{\epsilon_0}{2}\int_V \vec E\cdot \vec E \,dV \tag 1$$
where $V$ is all of space.
We now assume that the field is induced by a charge distribution residing on a conductor of volume $V$ bounded by surface $S$. Inasmuch as the electric field is zero inside the conductor, the volume of integration in $(1)$ can be taken as all of space outside of the conductor.
We can change the form of $(1)$ using $\vec E=-\nabla \Phi$ along with the vector identity
$$\nabla \cdot \Phi\,\vec E = \Phi \nabla \cdot \vec E+\nabla \Phi \cdot \vec E\tag2$$
Then, we have for any conducting surface $S$
$$\begin{align} W&=\frac{\epsilon_0}{2}\int_V \vec E\cdot \vec E \,dV \\\\ &=-\frac{\epsilon_0}{2}\int_V \nabla \Phi \cdot \vec E \,dV \tag 3\\\\ &=\frac{\epsilon_0}{2}\int_V \Phi \nabla \cdot \vec E-\nabla \cdot (\Phi \vec E) \,dV \tag 4\\\\ &=\frac{\epsilon_0}{2}\oint_S\Phi (\hat n \cdot \vec E)dS \tag 5\\\\ &=\frac{\epsilon_0}{2}\left. \Phi(\vec r)\right|_{\vec r\in S}\frac{Q}{\epsilon_0} \tag 6\\\\ &=\frac12Q\left. \Phi(\vec r)\right|_{\vec r\in S} \tag 7 \end{align}$$
NOTES:
In arriving at $(3)$, we used $\vec E = -\nabla \Phi$.
In going from $(3)$ to $(4)$, we used the vector identity $(2)$.
In going from $(4)$ to $(5)$, we used the Divergence Theorem, recognizing that the unit normal points out of $V$ and therefore into the surface of the conductor; we then absorbed the minus sign wherein the final result, the unit normal points out of the conductor's surface. We also used that Gauss's Law, $\nabla \cdot \vec E=\rho/\epsilon_0=0$ in $V$.
In going from $(5)$ to $(6)$, we used the fact that the potential is constant on the conducting surface and the integral form of Gauss's Law, $\oint \vec E\cdot \hat n dS=Q/\epsilon_0$.
Thus, the electrostatic energy stored is given by
$$\bbox[5px,border:2px solid #C0A000]{W=\frac12 Q\left. \Phi(\vec r)\right|_{\vec r\in S}} \tag 8$$
For a conducting sphere, the potential on the surface is $\Phi = \frac{Q}{4\pi \epsilon_0R}$ which from $(8)$, immediately gives the result
$$W=\frac{Q^2}{8\pi\epsilon_0R}$$
as expected!
So, we see that one does not need to integrate $\vec E\cdot \vec E$ to determine the stored energy! Rather, we can simply use $(8)$ directly.