To start, I'm going to characterize angle- and orientation-preserving slightly differently. I will say that a function is angle preserving if
$$\frac{\langle(f\gamma_1)'(0),(f\gamma_2)'(0)\rangle}{\vert(f\gamma_1)'(0)\vert\vert(f\gamma_2)'(0)\vert}=\frac{\langle\gamma_1'(0),\gamma_2'(0)\rangle}{\vert\gamma_1'(0)\vert\vert\gamma_2'(0)\vert}.$$
It's what you get if you characterize the angle via the law of cosines, $\langle v,w\rangle=\vert v\vert\vert w\vert\cos\theta$, where $\theta$ is the (non-oriented) angle between $v$ and $w$. The above equation basically says that the cosine of the angle is preserved. And I call a map orientation preserving if for all $\gamma_1,\gamma_2$ whose derivatives at $0$ are linearly independent, the linear map characterized by $\gamma_1'(0)\mapsto(f\gamma_1)'(0)$ and $\gamma_2'(0)\mapsto(f\gamma_2)'(0)$ has positive determinant. This is because the determinant has the orientation of an ordered set of vectors baked into its sign by construction.
With these preliminaries out of the way, we can first notice that $f$ is real differentiable because it has continuous partial derivatives in an open set according to your premises. So $f$ has a differential $\mathrm Df$, and we can use the chain rule to obtain
$$(f\gamma_i(0))'=\mathrm (Df)(\gamma_i(0))\cdot\gamma_i'=\mathrm Df(z_0)\cdot\gamma_i'.$$
So we have
$$\frac{\langle \mathrm Df(z_0)\cdot\gamma_1'(0),\mathrm Df(z_0)\cdot\gamma_2'(0)\rangle}{\vert\mathrm Df(z_0)\cdot\gamma_1'(0)\vert\vert\mathrm Df(z_0)\cdot\gamma_2'(0)\vert}=\frac{\langle\gamma_1'(0),\gamma_2'(0)\rangle}{\vert\gamma_1'(0)\vert\vert\gamma_2'(0)\vert}.$$
And this holds for all smooth paths $\gamma_1,\gamma_2$, and $\gamma_{1/2}'(0)$ can be any vector except $0$. So $\mathrm Df(z_0)$ belongs to a class of linear maps $L:\mathbb R^2\to\mathbb R^2$ which satisfy
$$\frac{\langle Lv,Lw\rangle}{\vert Lv\vert\vert Lw\vert}=\frac{\langle v,w\rangle}{\vert v\vert\vert w\vert},\qquad\textrm{for all }v,w\in\mathbb R^2\backslash\{0\}.$$
Such linear maps are called conformal if they are orientation preserving and anticonformal if they are not. We will now find out how conformal maps look. First thing, we will show that $\vert Lv\vert$ is proportional to $\vert v\vert$. In particular, $\vert Lv\vert=\lambda\vert v\vert$ for some $\lambda\in\mathbb R$, independent of what $v$ actually is. For proof, consider two linearly independent vectors $v,w$ with $\vert v\vert=\vert w\vert$. We will use the fact that $v+w$ and $v-w$ are orthogonal to show that $\vert Lv\vert=\vert Lw\vert$. It holds that
$$\frac{\langle L(v+w),L(v-w)\rangle}{\vert L(v+w)\vert\vert L(v-w)\vert}=\frac{\langle v+w,v-w\rangle}{\vert v+w\vert\vert v-w\vert}.$$
Juggling a bit with the rules for inner products and linear maps we get
$$\frac{\langle Lv,Lv\rangle-\langle Lw,Lw\rangle}{\vert L(v+w)\vert\vert L(v-w)\vert}=\frac{\langle v,v\rangle-\langle w,w\rangle}{\vert v+w\vert\vert v-w\vert}.$$
Since $\vert v\vert^2:=\langle v,v\rangle$, and by assumption $\vert v\vert=\vert w\vert$, the right side is $0$. But then the left side is also $0$, meaning that
$$\begin{align*}\langle Lv,Lv\rangle&=\langle Lw,Lw\rangle\\
\vert Lv\vert^2&=\vert Lw\vert^2\\
\vert Lv\vert&=\vert Lw\vert.
\end{align*}$$
And since we chose $v,w$ in arbitrary direction, $\vert Lv\vert$ can thus only depend on $\vert v\vert$. For a linear map, this means $\vert Lv\vert=\lambda\vert v\vert$ for some $\lambda\geq0$. Also note that $\lambda\neq0$ because otherwise the defining equation could never be fulfilled, since we would have to divide by $0$. Now with this information, the defining equation becomes
$$\begin{align*}\frac{\langle Lv,Lw\rangle}{\lambda\vert v\vert\lambda\vert L\vert}&=\frac{\langle v,w\rangle}{\vert v\vert\vert w\vert}\\
\frac{\langle \lambda^{-1}Lv,\lambda^{-1}Lw\rangle}{\vert v\vert\vert w\vert}&=\frac{\langle v,w\rangle}{\vert v\vert\vert w\vert}\\
\langle\lambda^{-1}Lv,\lambda^{-1}Lw\rangle&=\langle v,w\rangle.
\end{align*}$$
This makes $\lambda^{-1}L$ an orthogonal transformation! And it's already known how orthogonal transformations on $\mathbb R^2$ look (the Wikipedia article on orthogonal matrices has descriptions). The ones with positive determinant (so the orientation preserving ones) look like this:
$$\begin{pmatrix}a&-b\\b&a\end{pmatrix},\qquad a^2+b^2=1.$$
And if $\lambda^{-1}L$ looks like that, then $L$ looks like this:
$$\begin{pmatrix}a&-b\\b&a\end{pmatrix},\qquad a^2+b^2=\lambda^2.$$
So that's how $\mathrm Df$ looks. You can immediately read out the Cauchy-Riemann equations from the components now, and also note that the condition $a^2+b^2=\lambda^2\neq0$ enforces that $\mathrm Df\neq0$, which translates to $f'\neq0$. Though by the way, the Cauchy-Riemann equations are just a means to verify that $\mathrm Df$ is of the above form, which is the actual condition needed for holomorphy. So we could also skip the CR equations and just directly say that $\mathrm Df$ is of the form required for holomorphy.
