Criteria of the holomorphic subbundle Note that $\mathcal{E} = (E, \bar{\partial}_{\mathcal{E}})$ is a holomorphic vector bundle over complex manifold $X$, where $\bar{\partial}_{\mathcal{E}}$ is a integrable Dolbeault operator on $E$.
Now, I consider a $h_0$-orthogonal projection $\pi \in C^{\infty}(End(E))$, that is $\pi^\ast = \pi = \pi^2$, where $h_0$ is a Hermitian metric on $\mathcal{E}$.
Why the following statement is true: if $\pi$ satisfies 
$$(Id_{\mathcal{E}} - \pi) \circ \bar{\partial}_{\mathcal{E}} \circ \pi = 0,$$
then $F := im(\pi)$ is a holomorphic subbundle on $\mathcal{E}$?
 A: Here's a sketch of the proof. Of course, it's a local question to see that $F$ has a holomorphic structure. Start with an adapted $h_0$-unitary frame $e_1,\dots,e_k,e_{k+1},\dots,e_m$ for $\mathcal E$ over an open subset $U\subset X$, where $e_1,\dots,e_k$ give a frame for $F$ and $e_{k+1},\dots,e_m$ give a frame for $F^\perp$. Note, first of all, that if $F$ is indeed a holomorphic subbundle, we can choose a holomorphic frame $Z_1,\dots,Z_k$ and write $e_\alpha = \sum a_\alpha^\beta Z_\beta$ for some smooth functions $a_\alpha^\beta$, $1\le \alpha,\beta\le k$. Then $\bar\partial_{\mathcal E} e_\alpha = \sum \bar\partial a_\alpha^\beta\otimes Z_\beta$ (since $\bar\partial_{\mathcal E} Z_\beta = 0$ by hypothesis). This means, in particular, that the projection of $\bar\partial_{\mathcal E}e_\alpha$ into $F^\perp$ is $0$, in agreement with the criterion.
To prove the desired result, we suppose that for $1\le \alpha \le k$ we have $\bar\partial_{\mathcal E} e_\alpha = 0 \pmod F$. This means that there are $(0,1)$-forms $\phi_\alpha^\beta$ so that
$$\bar\partial_{\mathcal E} e_\alpha = \sum \phi_\alpha^\beta e_\beta.$$
We want to show that we can find smooth functions $f_\alpha^\beta$ so that the $Z_\alpha = \sum f_\alpha^\beta e_\beta$ are holomorphic. That is, we want to solve
$$0 = \bar\partial_{\mathcal E} Z_\alpha = \sum_\beta\left(\bar\partial f_\alpha^\beta + \sum_\gamma f_\alpha^\gamma\phi_\gamma^\beta\right) e_\beta.$$
That is, given the $(0,1)$-forms $\phi_\alpha^\beta$, we want to solve $$\bar\partial f_\alpha^\beta + \sum_\gamma f_\alpha^\gamma\phi_\gamma^\beta = 0 \quad\text{for all } \alpha,\beta = 1,\dots,k.\tag{$\star$}$$
Now, as usually happens in differential geometry, the matrix $\phi = (\phi_\alpha^\beta)$ is not completely arbitrary. It must satisfy an integrability condition. Because $\bar\partial_{\mathcal E}^2 = 0$, it follows that $\bar\partial\phi - \phi\wedge\phi = 0$, and so 
$$d\phi - \phi\wedge\phi = \Phi$$
is a $k\times k$ matrix of $(1,1)$-forms. This, in turn, is the integrability condition needed to solve ($\star$).
With thanks to Robert Bryant for helping, here's roughly how this goes. By a clever trick, we can turn the $\bar\partial$-system of differential equations into a usual $d$ differential system. Let $(z^1,\dots,z^n)$ be holomorphic coordinates on $U\subset X$ and let $w^1,\dots,w^k$ be coordinates on $\Bbb C^k$.
Consider the ideal $\mathscr I$ of complex differential forms on $U\times\Bbb C^k$ generated by the $1$-forms $dz^1,\dots,dz^n$ and $\omega^\beta = dw^\beta + \sum w^\gamma\phi_\gamma^\beta$. Well, of course, $d(dz^j)=0$. Next, switching to matrix notation, thinking of $w=(w^1,\dots,w^k)$ as a row vector, we can rewrite this as $\omega = dw + w\phi$. Now, differentiating, we have
\begin{align*}
d\omega &= dw\wedge\phi + w\,d\phi = (\omega- w\phi)\wedge\phi + w(\phi\wedge\phi+\Phi) \\
&=\omega\wedge\phi + w\Phi = 0 \pmod{\mathscr I}
\end{align*}
inasmuch as $\Phi$, being a matrix of $(1,1)$-forms on $U$, is in the ideal generated by the $dz^j$. In summary, $d\mathscr I\subset\mathscr I$. We can now apply Nirenberg's complex Frobenius theorem: Since, moreover, $\mathscr I\cap\overline{\mathscr I} = \{0\}$, it follows that there are complex coordinates $x^1,\dots,x^{n+k}$ on a neighborhood $W\subset U\times\Bbb C^k$ of our starting point so that $\mathscr I$ is generated by $dx^1,\dots,dx^{n+k}$. This now endows $W$ with a (different) complex structure. Since projection $W\to U$ is a holomorphic submersion, we can suppose $x^j = z^j$ for $j=1,\dots,n$. Moreover, we can (shrinking open sets if necessary) choose $m$ independent holomorphic sections $s^\beta$ of this projection and write $s^\beta(z) = (z,f^\beta(z))$ for smooth functions $f^\beta\colon \tilde U\to\Bbb C^k$. 
Since the $\omega^\beta$ are in the ideal generated by the $dx^i$, it follows that $df^\beta + \sum f^\gamma\phi_\gamma^\beta$ are in the ideal generated by $dz^1,\dots, dz^n$, which tells us that, indeed, $\bar\partial f^\beta + \sum f^\gamma\phi_\gamma^\beta = 0$, as desired.
