[proofplan]
We choose a small open rectangle around $z_0$ contained in $\Omega$ and construct $v$ by integrating the closed one-form $-u_y\,dx+u_x\,dy$ along a fixed rectangular path from $z_0$ to the point. The definition makes $\partial_y v=u_x$ directly, while differentiating the first integral and using $\Delta u=0$ gives $\partial_x v=-u_y$. These are exactly the Cauchy-Riemann equations for $u+iv$, so the Cauchy-Riemann holomorphicity criterion finishes the proof.
[/proofplan]
[step:Choose a coordinate rectangle contained in $\Omega$]
Write $z=x+iy$ and write $z_0=x_0+iy_0$, where $x,y,x_0,y_0 \in \mathbb R$. Since $\Omega \subset \mathbb C$ is open and $z_0 \in \Omega$, there exists $\rho>0$ such that the open square
\begin{align*}
U:=\{x+iy \in \mathbb C: |x-x_0|<\rho,\ |y-y_0|<\rho\}
\end{align*}
satisfies $U \subset \Omega$. Let
\begin{align*}
I&:=(x_0-\rho,x_0+\rho),&
J&:=(y_0-\rho,y_0+\rho).
\end{align*}
We regard $u|_U$ as a function
\begin{align*}
u: I\times J &\to \mathbb R \\
(x,y)&\mapsto u(x+iy).
\end{align*}
Because $u \in C^2(\Omega;\mathbb R)$ is harmonic, the partial derivatives $u_x,u_y,u_{xx},u_{yy}:I\times J\to \mathbb R$ are continuous and satisfy
\begin{align*}
u_{xx}(x,y)+u_{yy}(x,y)=0
\end{align*}
for every $(x,y)\in I\times J$.
[/step]
[step:Define the candidate conjugate by integrating the closed differential form]
For an oriented one-dimensional integral, we use the convention
\begin{align*}
\int_a^b g(t)\,d\mathcal L^1(t)=-\int_b^a g(t)\,d\mathcal L^1(t)
\end{align*}
when $b<a$. Define
\begin{align*}
v:I\times J&\to \mathbb R\\
(x,y)&\mapsto -\int_{x_0}^{x} u_y(s,y_0)\,d\mathcal L^1(s)+\int_{y_0}^{y}u_x(x,t)\,d\mathcal L^1(t).
\end{align*}
The integrands are continuous because $u\in C^2(\Omega;\mathbb R)$, so both one-dimensional integrals are well-defined for all $(x,y)\in I\times J$.
[guided]
The differential form suggested by the Cauchy-Riemann equations is
\begin{align*}
-u_y\,dx+u_x\,dy.
\end{align*}
If $v$ is a harmonic conjugate of $u$, then the Cauchy-Riemann equations should read
\begin{align*}
v_x=-u_y,\qquad v_y=u_x.
\end{align*}
Thus $dv$ should equal $-u_y\,dx+u_x\,dy$. Locally, we can build such a potential by integrating this form along a fixed path from the base point $(x_0,y_0)$ to $(x,y)$.
We choose the rectangular path that first moves horizontally from $(x_0,y_0)$ to $(x,y_0)$ and then vertically from $(x,y_0)$ to $(x,y)$. With the oriented integral convention
\begin{align*}
\int_a^b g(t)\,d\mathcal L^1(t)=-\int_b^a g(t)\,d\mathcal L^1(t)
\end{align*}
for $b<a$, define
\begin{align*}
v:I\times J&\to \mathbb R\\
(x,y)&\mapsto -\int_{x_0}^{x} u_y(s,y_0)\,d\mathcal L^1(s)+\int_{y_0}^{y}u_x(x,t)\,d\mathcal L^1(t).
\end{align*}
The first integral is the contribution of $-u_y\,dx$ along the horizontal segment, and the second integral is the contribution of $u_x\,dy$ along the vertical segment. Since $u\in C^2(\Omega;\mathbb R)$, the functions $(s,y_0)\mapsto u_y(s,y_0)$ and $(x,t)\mapsto u_x(x,t)$ are continuous on their respective intervals, so the one-dimensional integrals are well-defined.
[/guided]
[/step]
[step:Differentiate $v$ and use harmonicity to obtain the Cauchy-Riemann equations]
By the [Fundamental Theorem of Calculus](/theorems/632) applied to the continuous map $t\mapsto u_x(x,t)$ on $J$, for every $(x,y)\in I\times J$,
\begin{align*}
v_y(x,y)=u_x(x,y).
\end{align*}
For the $x$-derivative, differentiating the first integral by the [Fundamental Theorem of Calculus](/theorems/632) and differentiating the second integral under the integral sign, which is valid because $u_{xx}$ is continuous on $I\times J$, gives
\begin{align*}
v_x(x,y)
&=-u_y(x,y_0)+\int_{y_0}^{y}u_{xx}(x,t)\,d\mathcal L^1(t).
\end{align*}
Since $u$ is harmonic, $u_{xx}(x,t)=-u_{yy}(x,t)$ for every $(x,t)\in I\times J$. Hence
\begin{align*}
v_x(x,y)
&=-u_y(x,y_0)-\int_{y_0}^{y}u_{yy}(x,t)\,d\mathcal L^1(t)\\
&=-u_y(x,y_0)-\bigl(u_y(x,y)-u_y(x,y_0)\bigr)\\
&=-u_y(x,y),
\end{align*}
where the second equality is the [Fundamental Theorem of Calculus](/theorems/632) applied to the continuous map $t\mapsto u_{yy}(x,t)$. Therefore
\begin{align*}
u_x=v_y,\qquad u_y=-v_x
\end{align*}
on $I\times J$.
[/step]
[step:Apply the Cauchy-Riemann criterion to conclude holomorphicity]
Define
\begin{align*}
F:U&\to\mathbb C\\
x+iy&\mapsto u(x,y)+iv(x,y).
\end{align*}
The preceding step shows that $v$ has continuous first partial derivatives on $I\times J$ because $v_x=-u_y$ and $v_y=u_x$, and $u_x,u_y$ are continuous. Hence $F\in C^1(U;\mathbb C)$. The identities
\begin{align*}
u_x=v_y,\qquad u_y=-v_x
\end{align*}
are precisely the Cauchy-Riemann equations for $F=u+iv$. By the [Cauchy-Riemann Holomorphicity Criterion](/theorems/???), a $C^1$ complex-valued function on an open subset of $\mathbb C$ satisfying the Cauchy-Riemann equations is holomorphic. Therefore $F=u+iv$ is holomorphic on $U$, where $U\subset\Omega$ is an open neighbourhood of $z_0$. This proves the theorem.
[/step]