Energy functional in Poisson's equation: what physical interpretation?

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}$.

In the Lagrangian formulation of classical mechanics, the functional that we extremise to derive the Euler-Lagrange equation of motion is the Lagrangian, defined as

$$L[u] = T[u] - V[u]$$

where $T[u]$ is the kinetic part and $V[u]$ is the potential part. Contrast this with the Hamiltonian $H=T+V$ from which we derive Hamilton's equations of motion (these are equivalent to the Euler-Lagrange equations, and the two are related by Legendre transform). That should explain the minus sign discrepancy.

I am not sure about the factor of 1/2. Is it possible that we are simply defining $\rho_{\mathrm{Feynman}} = 2\rho_{\mathrm{Evans}}$?

The energy functional has to agree with the total energy and there can't be any ambiguity or "freedom of conventions" that would add the factors of two, among other things. After all, the total energy/mass is a subject to conservation laws and it determines both inertial of the system as well as the strength of its gravitational field.

If you integrated Feynman's two $\int E^2/2$ and $\int \rho V/2$ contributions over the whole 3-space, you would surely be double-counting the energy. The correct energy is just one of them, as given by the Evans formula. In the localized case, the two expressions may be shown to be equal.

You may also use Feynman's formula but you must be careful to use the $\int E^2/2$ exclusively for the part of the electric fields that are "unrelated" to localized charged sources.