In which sense are Euler-Lagrange PDE's on fiber bundles quasi-linear? In what follows, all manifolds are smooth, Hausdorff, paracompact, connected and oriented, and all maps between any two of them are assumed to be smooth. Let $\pi:E\rightarrow M$ be a fiber bundle over the base manifold $M$, which we assume to have dimension $n$. Given a bundle map $$\mathscr{L}:J^kE\rightarrow\wedge^nT^*M$$ over $M$, where $J^kE$ is (the total space of) the jet bundle of order $k$ of $\pi$ - i.e. $\mathscr{L}$ is a Lagrangian density of order $k$ - we define the action functional $S_{\mathscr{L}}:C^\infty_c(M)\times\Gamma(\pi)\rightarrow\mathbb{R}$ associated to $L$ (here $C^\infty_c(M)$ is the space of smooth real-valued functions on $M$ with compact support, and $\Gamma(\pi)=\{\phi\in C^\infty(M,E)\ |\ \pi\circ\phi=\text{id}_M\}$ is the space of smooth sections of $\pi$) as $$S_{\mathscr{L}}(f,\phi)=\int_M f(j^k\phi)^*\mathscr{L}\ .$$ The Euler-Lagrange operator $E(\mathscr{L})$ associated to $\mathscr{L}$ is the (usually nonlinear) partial differential operator of order $\leq 2k$ defined by the formula $$\int_M E(\mathscr{L})[\phi](\vec{\phi})=\int_M\left.\frac{\partial}{\partial t}\right|_{t=0}f(j^k\phi_t)^*\mathscr{L}\ ,\quad\vec{\phi}\in\Gamma_c(\phi^*VE\rightarrow M)\ ,f|_{\text{supp}\vec{\phi}}\equiv 1\ ,$$ where $\phi_t=\Phi(t,\cdot)$ for any $\Phi\in C^\infty(I\times M,E)$ with $0\in I\subset\mathbb{R}$ an open interval, $\phi_0=\phi$, $\left.\frac{\partial\Phi}{\partial t}\right|_{t=0}=\vec{\phi}$ and $\phi_t\in\Gamma(\pi)$, $\phi_t|_{M\smallsetminus\text{supp}\vec{\phi}}=\phi|_{M\smallsetminus\text{supp}\vec{\phi}}$ for all $t\in I$. The above definition of $E(\mathscr{L})$ can be shown to be independent of the particular choice of $f$ and $\Phi$ under the above conditions. There are other possible ways to define $E(\mathscr{L})$ in the literature, all leading to the same object.
One can see from the above definition that $E(\mathscr{L})[\phi]\in\Gamma(\phi^*V^\circledast E\rightarrow M)$, where $\pi_{V^\circledast E}:V^\circledast E=\pi^*(\wedge^nT^*M)\otimes V^*E\rightarrow E$ is the so-called twisted dual of the vertical bundle $VE=\ker T\pi\rightarrow E$ of $E$. Since $$\phi^*V^\circledast E=\{(p,q)\in M\times V^\circledast E\ |\ \phi(p)=\pi_{V^\circledast E}(q)\}\subset M\times_M V^\circledast E$$ for all $\phi\in\Gamma(\pi)$, it is clear that $\phi^*V^\circledast E$ is a sub-bundle of $M\times_M V^\circledast E$ over $M$. As such, we can write $$E(\mathscr{L})=\rho(E(\mathscr{L}))\circ j^{2k}\ ,$$ where $$\rho(E(\mathscr{L})):J^{2k}(E)\rightarrow M\times_M V^\circledast E$$ is a bundle map covering $\text{id}_M$ ($M$ is seen above as a fiber bundle over itself with singleton fibers and $\text{id}_M$ as the projection map).
When one writes a local formula for $E(\mathscr{L})[\phi]$ using a local trivialization of $E$ and a local chart for the typical fiber $Q$ of $\pi$, one sees that $E(\mathscr{L})[\phi]$ is always an affine function of the highest-order derivatives of $\phi$ if the order of $E(\mathscr{L})$ happens to be nonzero (which we assume to be the case), so it would make sense to say that $E(\mathscr{L})$ is a "quasi-linear" partial differential operator. However, $M\times_M V^\circledast E$ is generally not a vector or affine bundle over $M$ (unless $E$ itself is), despite $\phi^*V^\circledast E$ being a (different) vector bundle over $M$ for each $\phi\in\Gamma(\pi)$, so $E(\mathscr{L})$ does not fit into the usual way to globally define quasi-linear partial differential operators (see e.g. the discussion in Section IX.2, pp. 393ff of the book Exterior Differential Systems by R.L. Bryant et al.), which requires the target bundle of the operator to be a vector (or at least an affine) bundle.

Question: in which global (i.e. coordinate-independent) way one may define quasi-linear partial differential operators in order to encompass Euler-Lagrange operators $E(\mathscr{L})$ acting on smooth sections of general fiber bundles?

 A: Let $\mathcal{E} \subset J^{2k}E$ be the submanifold (provided that this subset is a submanifold) of all $2k$-jets sitting in the zero-level set of $\rho(E(\mathscr{L}))$, the PDE submanifold. This is Definition 8 in my [arXiv:1211.1914]:

A PDE submanifold $\mathcal{E} \subset J^{l}E$ is called linear if $E\to M$ is a vector bundle and $\mathcal{E} \to M$ is a vector sub-bundle of the $J^{l}E \to M$ vector bundle. The PDE system is called quasilinear if $\mathcal{E} \to J^{l-1}E$ is an affine sub-bundle of the affine bundle $J^{l}E \to J^{l-1}E$.

In my notation all the bundle projections (as well as their restrictions to sub-bundles) like $E\to M$ and the naturally induced $J^{l}E\to M$ or $J^{l}E\to J^{l-1}E$ are anonymous. One also needs to recall that the projection $J^{l}E \to J^{l-1}E$ is naturally an affine bundle.
That definition is immediately followed by the obvious Lemma 1:

A quasilinear PDE submanifold $\mathcal{E} \subset J^{l}E$ can always be represented as the kernel of a morphism of affine bundles. Namely, there exists a vector bundle $F\to M$ (naturally interpreted also as an affine bundle) and a morphism $f$ such that the following diagram of affine bundle morphisms is fiber-wise exact:
  $$\require{AMScd}
\begin{CD}
  \mathcal{E} @>>> J^{l}E @>\xrightarrow{~0~}>f> F \\
  @VVV @VVV @VVV \\
  J^{l-1}E @= J^{l-1}E @>>> M
\end{CD}$$

Now, I don't think that I had found the same exact definition anywhere else (that's why I don't give any reference), but I think that it's the one that makes the most sense.
A: Igor gave a coordinate independent definition of quasilinear equations  which the Euler-Lagrange equations satisfy. Still missing is a definition of quasilinear differential operator, which the Euler-Lagrange operator satisfies. 
You already observed that

$M\times_M V^\circledast E$ is generally not a vector or affine bundle over $M$.

hence the definition of quasilinear operator found in Bryant, Chern et. al. on p. 397 (which agrees with the definition found in Krasil'shchik, Lychagin, Vinogradov, Geometry of jet spaces and nonlinear partial differential equations p.160) cannot be applied directly, since they assume vector bundles over $M$. But $V^\circledast E$ is a vector bundle over $E$, and the fix might be to expand the definition of nonlinear DO:  
Definition: Let $V\to E$ be a vector bundle over $E$ and $\psi:J^k E \to V$ a morphism of fiber bundles over $E$ (not neccessarily linear, since we don't assume a linear structure on $E\to M$), then the map $\psi \circ j^k$ is called a nonlinear DO. Call such an operator quasilinear if $\psi$ is affine, when restricted to any fiber of $J^k E \to J^{k-1} E$.
I haven't checked, but suspect that the E-L-operators satisfies this definition of quasilinear DO (with $V=V^\circledast E$). Observe that the zero set of a quasilinear $\psi$ is a quasilinear PDE $\mathcal{E}\subset J^k E$, as in Igors answer. Observe also that a nonlinear DO in the more restrictive sense of: a map of fiber bundles $\rho: J^k E \to W$ over $M$ is a special case of the previous definition by pulling back $W$ to $E$.
Note also that such a generalization of DO is unavoidable in cases like minimal surfaces, where we don't have a fiber bundle structure on the ambient space $E$. (In fact, in the case of jets of submanifolds we probably need to allow vector bundles $V$ over $J^1 E$ instead of over $E$, since the vertical bundle $VE$ doesn't exist and is replaced by a "normal" bundle.)
