Consider an action that is gauge invariant. Do we obtain the same information from the following:
Find the equations of motion, and then fix the gauge?
Fix the gauge in the action, and then find the equations of motion?
Consider an action that is gauge invariant. Do we obtain the same information from the following:
Find the equations of motion, and then fix the gauge?
Fix the gauge in the action, and then find the equations of motion?
Here we will assume that we ultimately want to consider the full quantum theory, usually written in terms of a gauge-fixed path integral $$Z~=~\int \!{\cal D}\phi~ \exp\left(\frac{i}{\hbar}S_{\rm gf}[\phi]\right) \tag{1}$$ rather than just the classical action and the corresponding classical equations of motion (with or without gauge-fixing terms). If the gauge orbits have infinite volume (as is often the case), then we need to gauge-fix the path integral.
Of course, if we by brute force eliminate a field in the action via a gauge-fixing condition, then we can no longer carry out a variation of the action with respect to that field, and we would have lost information.
However, here we will only consider a 'softer' way to impose the gauge-fixing conditions via Lagrange multipliers, which may appear linearly or quadratically in the gauge-fixed action. The linear case leads directly to delta functions in the path integral, which impose the gauge-fixing conditions; while the quadratic case leads to Gaussian terms in the path integral, which suppress (but do not completely forbid) field configurations that violate the gauge-fixing condition. (Nevertheless, in a certain scaling limit, the Gaussian factors become delta functions.)
Together with the original fields, the (non-propagating, auxiliary) Lagrange multipliers are part of the fields $\phi$ in the path integral that is integrated over. In particular, the gauge-fixed action $S_{\rm gf}[\phi]$ can be varied with respect to these Lagrange multiplier fields as well.
These 'softly imposed' gauge-fixing conditions still do affect the variation of the action (as opposed to not imposing the gauge-fixing conditions). However, a more relevant question is:
Do the gauge-invariant physical observables of the theory depend on the specific gauge-fixing condition (e.g. Lorenz gauge, Coulomb gauge, etc)?
The answer is No, i.e. within the class of consistent gauge-fixing terms in the action, the specific form of the gauge-fixing terms in the corresponding equations of motion has no physical consequences.
For more general gauge theories, the equations of motion are e.g. not gauge-invariant, and it is better to encode the gauge symmetry via a (generalized) fermionic nilpotent BRST symmetry ${\bf s}$ that squares to zero $${\bf s}^2~=~0,\tag{2}$$ and preserves the original action $${\bf s} S_0~=~0.\tag{3}$$ The gauge-fixed action $$S_{\rm gf}=S_0+{\bf s}\psi\tag{4}$$ is the original action $S_0$ plus a BRST-exact term ${\bf s}\psi$ that depends on the so-called gauge-fixing fermion $\psi$, which encodes the gauge-fixing conditions.
A physical observable $F=F[\phi]$ in the theory is by definition required to be BRST-closed $${\bf s} F~=~0.\tag{5}$$
Statement. If the path integral measure is BRST-invariant, then the correlation function for a physical observable $$\langle F[\phi] \rangle ~=~ \frac{ \int \!{\cal D}\phi~ F[\phi]\exp\left(\frac{i}{\hbar}S_{\rm gf}[\phi]\right)} {\int \!{\cal D}\phi~ \exp\left(\frac{i}{\hbar}S_{\rm gf}[\phi]\right)} \tag{6}$$ is independent of the gauge-fixing fermion $\psi$.
Lemma 1. If $F$ is BRST-closed and ${\bf s} G$ is BRST-exact, then the product $F{\bf s} G$ is BRST-exact.
Lemma 2. The correlation function $$\langle {\bf s} F \rangle ~=~0\tag{7}$$ of a BRST-exact observable vanishes.
The lemma 2 follows by assuming ${\bf s}$ to be a Hermitian (or anti-Hermitian) operator.
Proof of statement: If $\widetilde{\psi}$ is another gauge-fixing fermion, then one may show via lemma 1 that the difference between the $F$-observables for different gauge-fixing fermions $\widetilde{\psi}$ and $\psi$ is $$ \left\{ \exp\left(\frac{i}{\hbar}{\bf s}(\widetilde{\psi}-\psi)\right)-1\right\}F~=~{\bf s}(\ldots), \tag{8}$$ which is BRST-exact, and therefore has a vanishing correlator function, cf. lemma 2. $\Box$
We should stress that the gauge-fixing fermion is required to satisfy certain rank conditions, and it can e.g. not be chosen to be identically zero.
Finally, let us mention, that an even bigger class of Lagrangian gauge theories can be treated with the help of the Batalin-Vilkovisky (BV) formalism.
No, it is not always consistent to first fix the gauge before deriving equations of motion. Consider electromagnetism coupled to matter. One can perform a gauge transformation to set $A_0=0$. However, if this is done in the action before deriving the equations of motion, then one will miss out on the $A_0$ equation of motion which guarantees the conservation of electric charge.
When attempting to derive sufficiently symmetric solutions, it is sometimes possible to begin by choosing an appropriate ansatz, even at the level of the action. But this is something of an art, and is not guaranteed to be consistent with the full equations of motion evaluated on that same ansatz.