For this question, I think two types of equations should be discussed together.
- Equations of motion for the fundamental fields in a Lagrangian theory. These are what comes up in the path integral formulation which tells you how probable a certain deviation should be.
- Conservation equations satisfied by polynomials in those fields that correspond to currents. Even though there could be anomalies, we should generically expect that these equations still have something to say about the quantum theory.
A few issues about quantization have already been brought up. For instance, a free quantum field is often expressed in terms of plane waves $e^{i p x}$ because these solve the Klein-Gordon equation. But exciting each plane wave individually is only something we can do because $\phi(x)$ is a simple sum of them and this is only true because the Klein-Gordon equation is linear. So in this sense, the fact that the Hilbert space is a Fock space is indeed a consequence of the equation of motion. It is therefore understandable that textbooks don't discuss interacting QFTs in the same way. If they did, they would essentially be limited to certain integrable QFTs that can be treated with a similar approach.
However, I think the biggest change in viewpoint comes about because we promote fields to operators. There is no way to measure an operator valued distribution directly. Experimentalists like S-matrix elements or Green's functions where we can plug in a position or a momentum and get a number. So we then have to ask whether equations of motion can be brought "outside the brackets" and taken to act on correlation functions themselves. The answer to this is yes but it comes with subtleties. For instance, the free field propagator is a Green's function of the Klein-Gordon equation so instead of
$$
(\partial_x^2 + m^2) \left < \phi(x) \phi(y) \right > = 0
$$
we have
$$
(\partial_x^2 + m^2) \left < \phi(x) \phi(y) \right > = \delta(x - y).
$$
This is what's known as an operator equation which turns into the classical equation if we only consider it at separated points. This distinction is related to the fact that two fields can be multiplied at will in a classical theory whereas composite operators are singular in a quantum theory and need to be renormalized. This happens for conservation laws too such as
$$
\partial_\mu \left < T^{\mu\nu}(x) O_1(x_1) \dots O_n(x_n) \right > = -\sum_i \delta(x - x_i) \partial^\nu_i \left < O_1(x_1) \dots O_n(x_n) \right >
$$
which is the operator equation for the stress tensor.
The thing is, once this leap to operator equations is made, they absolutely constrain the quantum theory in powerful ways. Ward identities in gauge theories, which relate different Green's functions, come about because of current conservation. Some scaling dimensions in conformal theories can be computed exactly because equations of motion are equivalent to shortening conditions. I.e. they specify that an operator must transform in a special representation of the conformal group where various descendants are guaranteed to vanish. This analysis is also important when we add interactions that break conformal invariance. We know that interactions lead to anomalous dimensions but this implies that short representations must become long. This again leads to constraints (sometimes powerful enough to obviate Feynman diagrams) because a generic Hilbert space would not provide the matter that the short multiplet needs to "eat". Finally, if one adopts superspace, all these short supermultiplets (with names like "half BPS" and such) can be regarded as solutions to equations of motion as well. And clever uses of these equations of motion is behind some modern methods of studying supersymmetric QFT by means of a "twist" which turns an otherwise complicated observable into something that is topological or holomorphic.
