A good rule of thumb for complex analysis is that, if something has gone wrong, $0$ is almost certainly to blame.
Here is the difference quotient for $\overline{z}$ at $0$ taken along two different paths in the complex plane:
Tending to $0$ along the positive imaginary axis produces:
$\lim_{h\rightarrow0}\frac{\overline{\left(0+ih\right)}-\overline{0}}{ih}=\lim_{h\rightarrow0}\frac{-ih}{ih}=-1$
On the other hand, tending to $0$ along the positive real axis produces:
$\lim_{h\rightarrow0}\frac{\overline{\left(0+h\right)}-\overline{0}}{h}=\lim_{h\rightarrow0}\frac{h}{h}=1$
which is completely different. Since the limits do not agree, the limit of this difference quotient as $h$ tends to $0$ in $\mathbb{C}$ does not exist. Thus, the conjugate function is not complex differentiable at $0$.
Generally, when working "from first principles" as it were, one has to rely on the two-dimensionality of $\mathbb{C}$ (along with its algebraic/"arithmetic" properties) to show that things are different than the case of single real-variable analysis. Indeed, it is exactly because the condition of holomorphy requires these two-dimensional limits to exist that holomorphic functions are so marvelously well-behaved (in stark contrast to functions of multiple real variables).
Once you understand the basic trick of using Hahn-Banach to reduce everything to the situation where the target space is $\Bbb{C}$, I'm sure you can easily extend several results by yourself. What I did a few months back (and what I recommend you do) is pick your favorite complex analysis text (e.g I took Cartan), and I went through each of the basic theorem one-by one and extended them to the case where the target is an arbitrary Banach space, and the domain is an open subset of $\Bbb{C}^n$. Usually such extensions aren't too hard. Of course, from time to time I would refer to Mujica's text to see if I was missing out any details in the argument.
As an example of how Hahn-Banach can be used here, I shall prove Cauchy's integral formula:
Let $U\subset \Bbb{C}$ be an open set, $X$ a complex Banach space, and $f:U\to X$ a holomorphic mapping. Let $z\in U$ and let $\gamma:[a,b]\to U$ be a $C^1$ loop such that $z\notin \text{image}(\gamma)$. Then,
\begin{align}
I(\gamma,z)f(z)&=\frac{1}{2\pi i}\int_{\gamma}\frac{f(\zeta)}{\zeta-z}\,d\zeta
\end{align}
where $I(\gamma,z)$ denotes the index of $\gamma$ wrt $z$ (which is $+1$ for a circle containing $z$ "inside" of it).
Just so we're clear, by holomorphic on $U$, I mean once-complex differentiable at every point of $U$ ($X$ is a Banach space, so the difference quotient makes sense and limits make sense). Also, on the RHS, I'm referring to the Bochner integral of a continuous function $[a,b]\to X$, namely $\frac{1}{2\pi i}\int_a^b\frac{f(\gamma(t))}{\gamma(t)-z}\,\gamma'(t)\,dt$.
We prove this by letting $\lambda\in X^*$ be arbitrary (i.e a continuous linear map $X\to\Bbb{C}$). Then, note that $\lambda\circ f:U\to\Bbb{C}$ is holomorphic (chain rule), so by the one-dimensional Cauchy-integral formula,
\begin{align}
\lambda\bigg(I(\gamma,z)f(z)\bigg)&=I(\gamma,z)(\lambda\circ f)(z)\\
&=\frac{1}{2\pi i}\int_{\gamma}\frac{(\lambda\circ f)(\zeta)}{\zeta-z}\,d\zeta\tag{by 1D version}\\
&=\frac{1}{2\pi i}\int_{\gamma}\lambda\left(\frac{f(\zeta)}{\zeta-z}\right)\,d\zeta\tag{since $\lambda$ is linear}\\
&=\lambda\left(\frac{1}{2\pi i}\int_{\gamma}\frac{f(\zeta)}{\zeta-z}\,d\zeta\right)
\end{align}
In the last line, I used one of the basic properties of Bochner integrals, that for a Bochner-integrable function $\phi:S\to X$ from a measure space $S$ to a Banach space $X$, and $\lambda\in X^*$, we have $\lambda(\int_S\phi\,d\mu)=\int_S(\lambda\circ \phi)\,d\mu$ (in our case, you may wish to write out $\int_{\gamma}$ as an integral over $[a,b]$ to see that I'm using the theorem with $S=[a,b]$ and $\mu$ being Lebesgue measure).
Now, since $\lambda\in X^*$ was arbitrary, the Hahn-Banach theorem tells us the two things inside of $\lambda$ must be equal:
\begin{align}
I(\gamma,z)f(z)&=\frac{1}{2\pi i}\int_{\gamma}\frac{f(\zeta)}{\zeta-z}\,d\zeta
\end{align}
From here of course, the fact that "holomorphic on a disc $D_R(a)$ implies analytic on the disc $D_R(a)$" follows immediately by Taylor expanding the integrand as in the single variable case (I'm being slightly sloppy here, see this answer for the correct details). We fix $0<r<R$, then (by the just-established Banach-valued case of Cauchy's integral formula)
\begin{align}
f(z)&=\frac{1}{2\pi i}\int_{|\zeta-a|=r}\frac{f(\zeta)}{\zeta-z}\,d\zeta\\
&=\frac{1}{2\pi i}\int_{|\zeta-a|=r}\frac{f(\zeta)}{\zeta-a}\cdot\frac{1}{1-\frac{z-a}{\zeta-a}}\,d\zeta\\
&=\frac{1}{2\pi i}\int_{|\zeta-a|=r}\frac{f(\zeta)}{\zeta-a}\sum_{n=0}^{\infty}\left(\frac{z-a}{\zeta-a}\right)^n\,d\zeta\\
&=\sum_{n=0}^{\infty}\left(\frac{1}{2\pi i}\int_{|\zeta-a|=r}\frac{f(\zeta)}{(\zeta-a)^{n+1}}\,d\zeta\right)\cdot (z-a)^n,
\end{align}
where, as explained in the link, the exchange of series integral is possible due to uniform convergence.
This proves "holomorphic implies analytic". The fact that "analytic implies holomorphic" is much easier; simply take your favorite complex-analysis book and look at the proof there. The fact that $f$ is Banach-valued makes no difference (just replace absolute values by norms in appropriate places).
Now that you have the equivalence of holomorphic and analytic, I doubt you'll need Goursat's lemma (but you can still prove it if you want; just use Hahn-Banach to reduce to the one-dimensional case, and then invoke the already-known version of Goursat's lemma). Likewise, you can prove Morera's theorem, and establish Cauchy's inequalities, Louiville's theorem. With Morera's theorem in the Banach-case, you can once again prove that uniform limits of holomorphic functions are holomorphic.
Best Answer
Briefly, holomorphic functions are immensely useful because, on the one hand, they are surprisingly common (since any power series, for example, whose coefficients grow reasonably slowly defines a holomorphic function), and on the other hand one can prove very strong theorems about them. There is a web of results including Cauchy's integral formula and the identity theorem which assert that holomorphic functions are astonishingly rigid: given information about a holomorphic function in a very small part of its domain, one can extract information about the function's behavior in other a priori unrelated parts of their domain (and this is what allows things like contour integration to work).
For that reason, holomorphic functions are a powerful tool to apply to a problem when they do apply. For example, analytic number theorists frequently construct holomorphic or meromorphic functions which carry number-theoretic information, such as the Riemann zeta function, to prove theorems like the prime number theorem. Since you say you have a background in combinatorics, you might enjoy reading Flajolet and Sedgewick's Analytic Combinatorics, a thorough exposition of (among other things) ways to use complex analysis to provide asymptotics for combinatorial sequences.
Here is a simple example. Let $f_n$ denote the number of ways that $n$ horses can win a race, with ties. It turns out that this sequence has generating function
$$F(z) = \sum_{n \ge 0} \frac{f_n}{n!} z^n = \frac{1}{2 - e^z}.$$
This function is meromorphic with poles at $z = \log 2 + 2 \pi i k, k \in \mathbb{Z}$, each of which has residue $-\frac{1}{2}$. In fact, it turns out that $F(z)$ admits an infinite partial fraction decomposition
$$F(z) = \sum_{k \in \mathbb{Z}} \frac{1}{2(\log 2 + 2 \pi i k - z)}.$$
And by expanding the terms on the RHS in a geometric series, this gives an asymptotic expansion for $\frac{f_n}{n!}$ with leading term $\frac{1}{2 (\log 2)^{n+1}}$. In other words,
$$f_n \sim \frac{n!}{2 (\log 2)^{n+1}}.$$
The pole at $z = \log 2$ dominates the asymptotic expansion: the leading term in the error of the above expression is given by the other poles nearest the origin, which occur at $z = \log 2 \pm 2 \pi i$. Because these poles have nonzero imaginary part, if you plot the error in the above approximation you'll find that it oscillates. It is not so easy to explain why this should be the case without complex analysis.
A famous example is Hardy and Ramanujan's asymptotic formula for the partition function
$$p(n) \sim \frac{1}{4n \sqrt{3}} e^{ \pi \sqrt{ \frac{2n}{3} } }$$
which is proven using a much more sophisticated version of the above argument.
But really, there is too much to say about holomorphic functions, so again I suggest that you read a textbook. Besides Needham's book, I also personally enjoyed Stein and Shakarchi, which is very user-friendly and has good applications.