The point of view you give of quantizing phase space is not so good. If you have a classical phase space, there will be more than one quantum analog, differing by commutator terms that vanish classically. I never found the idea that you start with a classical system particularly lucid--- the quantum system is defined by the path integral, and the answer must always be found in the path integral itself.
In 4d, as you said in the comments, the problem is sidestepped by using N=1 superspace plus the appropriate R symmetry, which rotate the different SUSY's into each other, ensuring that the full SUSY is there. In order for the theory to be SUSY, you just have to check that the theory is N=1 and that the R-symmetry is good. This is usually manifest, since the R symmmetry is simple global rotation of fields into each other. But this doesn't work in 10 or 11 dimensions.
A part of the answer suggested by Qmechanic, the Schwinger Dyson equation (Heisenberg equation of motion) is obeyed by the field operators, so the variation really does vanish as an operator, in a certain sense described precisely below. But this is not a complete answer because it isn't true that everything that vanishes classically when the equations of motion are satisfied also vanishes quantum mechanically, in the path integral. The reason is that the quantities in the path integral can be multiplied at the same point in space and time, and there are singularities of the coinciding operator products which make the naive equations of motion fail in exactly the type of expressions that occur in the action.
For a trivial example, consider the product of two X operators
$$\langle x(t) x(t')\rangle = 0 $$
In the trivial Euclidean free-particle (Brownian motion) 1-d (quantum mechanics, not field theory) path integral $S= \int \dot{x}^2$. The operator X obeys the equation of motion $\ddot{X}=0$, so the second derivative is zero, but of course it isn't at the coinciding point (or else propagator would vanish at all times with the reasonable boundary conditions). The correct second derivative of the correlation function is
$$ \langle \ddot{x}(t) x(t')\rangle = \delta(t-t')$$
The equation of motion holds except at coincident points. The same is true in field theories for the same reason--- the propagator is sourced at the point where the two operators coincide by a delta-function with unit amplitude (if the fields are canonically normalized).
Further, since $\ddot{X}=0$, you might think that the following transformation is a symmetry
$$ \dot x + \epsilon f(t)\ddot{X} $$
Since the variation of S to first order in $\epsilon$ is $\int \dot{x}\ddot{x}f(t)$. So naively, using the equation of motion $\ddot{x}=0$, you would naively find that this is a symmetry. But this is obvious nonsense--- if f is nonconstant (if f is constant, the transformation is a time translation) this is obviously not a symmetry of the action. since it is a time translation by a different amount at different times, and there is a unit of time defined by the diffusion.
The reason is that the $\ddot{x}$ is multiplied by a coinciding $\dot{x}$, and when there are coinciding quantities, the equations of motion are not necessarily satisfied. The actual variation in S is
$$ \int f(t) {d\over dt} {\dot{x}^2\over 2}$$
and you see that it is a perfect derivative if f is constant, assuming Stratonovich convention for time derivatives--- centered differences--- and this is the symmetry of time translation in this case), but it is not a symmetry if f is not constant.
When are equations of motion satisfied in a path integral?
When you have a path integral:
$$ \int e^{iS} D\phi $$
The integral is invariant under a shift of integration variables $\phi(x)+\delta\phi(x)$, where $\delta\phi$ is an arbitrary function of x, since each integral is separately translation invariant. There is no determinant for this transformation--- it's just a shift in the integration variable by a constant at each time.
The change in the integrand when you do the shift is, given by expanding the thing to first order.
$$ \int e^{iS(\phi+\delta\phi)} D\phi = \int e^{iS} D\phi + i\int (\int {\delta S\over \delta \phi(x)}\delta\phi(x) d^dx) e^{iS} D\phi $$
and from this, you learn that the equations of motion are satisfied. This is the correct demonstration of the equations of motion.
$$ {\delta S\over \delta\phi} = 0 $$
This is the equation of motion for the $\phi$ field only, only one field. Now suppose you have some other insertions (composite operators $O_k$):
$$ \int O_1(x_1) O_2(x_2) ... O_n(x_n) e^{iS} D\phi $$
Now doing the same shift, to first order, you find
$$ \int (O_1 +\delta O_1(x_1)) ... (O_n+\delta O_n(x_n)) e^{iS} (1+i\delta S) D\phi = \int e^{iS} D\phi $$
So you find that the equation of motion for the one field $\phi$ is satisfied even with composite operator insertions except at coinciding points:
$$ -i\, {\delta S\over \delta \phi(x)} \, O_{1}(x_1) ... O_{n}(x_n) = {\delta O_1 (x_1)\over \delta \phi(x)} O_2(x_2) ... O_n(x_n) + O_1(x_1){\delta O_2(x_2)\over \delta\phi(x)}...O_n(x_n) + ... + O_1(x_1)O_2(x_2) ... O_{n-1}(x_{n-1}){\delta O_n(x_n)\over \delta \phi(x)}.$$
This is the operator equation for the equation of motion multiplying arbitrary expressions of other fields. The expression is known as the Schwinger-Dyson equation.
It is not just true as a vacuum expectation value, because it is satisfied independent of the path integral boundary conditions, which don't have to be vacuum, and you don't need to do a path integral from time minus infinity to time infinity, the identity holds in an arbitrary state.
This Schwinger-Dyson argument for insertions of any type with lots of fields running around implies the following rules for using the equation of motion:
- The equations of motion derived from varying $\phi$ are satisfied as operator equations away from insertions.
- the equation of motion derived from varying $\phi$ are satisfied even at cross-insertions of other fundamental fields (independent path-integrations) other than $\phi$.
- the equations of motion fail only at insertions of composite or elementary operators which vary when you shift the field whose equation of motion you are talking about, and the failure is of the ward-identity type.
So it isn't true that anything that vanishes classically has no effect. But when you have a product of the equation of motion multiplying functions of fields which are different from the one you vary to get the equation of motion, these are still identically zero.
In the SUSY case, your variation involves products of fields with their SUSY partners, which are independent integration variables. The equations of motion are satisfied in the products that you reduce, so using the equations of motion in the SUSY closure is justified, but you check it theory by theory by hand, by looking to see which equation of motion you use, and what it is multiplying.
The interpretation is that when you do the SUSY transformation, you also have to slide the fields around a little bit in a measure preserving way to finish the transformation.
