[proofplan]
We prove the statement locally in holomorphic coordinate charts. In such charts, the canonical almost complex structures are represented by multiplication by $i$ on complex Euclidean space, so the displayed identity says exactly that the real differential of the coordinate expression is complex linear. Writing the coordinate expression in real and imaginary parts, this complex-linearity condition is equivalent to the several-variable [Cauchy-Riemann equations](/page/Cauchy-Riemann%20Equations). The forward direction follows from the Cauchy-Riemann equations for holomorphic coordinate functions, and the converse follows from the several-variable Cauchy-Riemann theorem.
[/proofplan]
[step:Translate the problem to holomorphic coordinate charts]
Fix a point $p \in X$. Let $n = \dim_{\mathbb{C}} X$ and $m = \dim_{\mathbb{C}} Y$. Choose a holomorphic chart $(U,\varphi)$ of $X$ with $p \in U$, where $\varphi: U \to \varphi(U) \subseteq \mathbb{C}^n$, and choose a holomorphic chart $(V,\psi)$ of $Y$ with $F(p) \in V$, where $\psi: V \to \psi(V) \subseteq \mathbb{C}^m$. After replacing $U$ by $U \cap F^{-1}(V)$, we may assume $F(U) \subset V$.
Define the coordinate representation
\begin{align*}
f: \varphi(U) \to \psi(V), \qquad f = \psi \circ F \circ \varphi^{-1}.
\end{align*}
Since $F$ is smooth and the coordinate charts are smooth maps, $f$ is a smooth map between open subsets of real Euclidean spaces after identifying $\mathbb{C}^n$ with $\mathbb{R}^{2n}$ and $\mathbb{C}^m$ with $\mathbb{R}^{2m}$.
Let $I_n: \mathbb{R}^{2n} \to \mathbb{R}^{2n}$ and $I_m: \mathbb{R}^{2m} \to \mathbb{R}^{2m}$ denote multiplication by $i$ under these identifications. Because $\varphi$ and $\psi$ are holomorphic coordinate charts, their tangent maps identify $J_{X,p}$ and $J_{Y,F(p)}$ with $I_n$ and $I_m$, respectively. Therefore
\begin{align*}
dF_p \circ J_{X,p} = J_{Y,F(p)} \circ dF_p
\end{align*}
holds if and only if
\begin{align*}
df_{\varphi(p)} \circ I_n = I_m \circ df_{\varphi(p)}.
\end{align*}
[guided]
We localize because both holomorphicity and the differential identity are local conditions. Fix $p \in X$, and choose holomorphic coordinate charts
\begin{align*}
\varphi: U \to \varphi(U) \subseteq \mathbb{C}^n
\end{align*}
around $p$ and
\begin{align*}
\psi: V \to \psi(V) \subseteq \mathbb{C}^m
\end{align*}
around $F(p)$. Replacing $U$ by $U \cap F^{-1}(V)$ ensures that $F(U) \subset V$, so the coordinate expression
\begin{align*}
f: \varphi(U) \to \psi(V), \qquad f = \psi \circ F \circ \varphi^{-1}
\end{align*}
is well-defined.
The point of using holomorphic charts is that the canonical almost complex structures become the standard complex structures on Euclidean space. Let $I_n: \mathbb{R}^{2n} \to \mathbb{R}^{2n}$ be multiplication by $i$ on $\mathbb{C}^n \cong \mathbb{R}^{2n}$, and let $I_m: \mathbb{R}^{2m} \to \mathbb{R}^{2m}$ be multiplication by $i$ on $\mathbb{C}^m \cong \mathbb{R}^{2m}$. The tangent map $d\varphi_p: T_pX \to \mathbb{R}^{2n}$ carries $J_{X,p}$ to $I_n$, and $d\psi_{F(p)}: T_{F(p)}Y \to \mathbb{R}^{2m}$ carries $J_{Y,F(p)}$ to $I_m$.
By the chain rule, the differential of the coordinate representation at $\varphi(p)$ is
\begin{align*}
df_{\varphi(p)} = d\psi_{F(p)} \circ dF_p \circ d(\varphi^{-1})_{\varphi(p)}.
\end{align*}
Thus the coordinate form of the identity
\begin{align*}
dF_p \circ J_{X,p} = J_{Y,F(p)} \circ dF_p
\end{align*}
is precisely
\begin{align*}
df_{\varphi(p)} \circ I_n = I_m \circ df_{\varphi(p)}.
\end{align*}
So the theorem reduces to the following Euclidean statement: a smooth map between open subsets of complex Euclidean spaces is holomorphic exactly when its real differential commutes with multiplication by $i$ at every point.
[/guided]
[/step]
[step:Identify complex linearity with the Cauchy-Riemann equations]
Let $z_j = x_j + i y_j$ be the standard complex coordinates on $\mathbb{C}^n$, where $x_j,y_j: \varphi(U) \to \mathbb{R}$ are the real coordinate functions for $1 \le j \le n$. Write the component functions of $f$ as
\begin{align*}
f_k: \varphi(U) \to \mathbb{C}, \qquad f_k = u_k + i v_k
\end{align*}
for $1 \le k \le m$, where $u_k,v_k: \varphi(U) \to \mathbb{R}$ are smooth real-valued functions.
At a point $a \in \varphi(U)$, the real differential $df_a: \mathbb{R}^{2n} \to \mathbb{R}^{2m}$ commutes with the standard complex structures,
\begin{align*}
df_a \circ I_n = I_m \circ df_a,
\end{align*}
if and only if, for every $1 \le j \le n$ and every $1 \le k \le m$,
\begin{align*}
\frac{\partial u_k}{\partial x_j}(a) = \frac{\partial v_k}{\partial y_j}(a)
\end{align*}
and
\begin{align*}
\frac{\partial u_k}{\partial y_j}(a) = -\frac{\partial v_k}{\partial x_j}(a).
\end{align*}
Indeed, $I_n(\partial_{x_j})=\partial_{y_j}$ and $I_n(\partial_{y_j})=-\partial_{x_j}$, while $I_m$ sends the target coordinate vector $\partial_{u_k}$ to $\partial_{v_k}$ and $\partial_{v_k}$ to $-\partial_{u_k}$. Comparing the $\partial_{u_k}$ and $\partial_{v_k}$ components of $df_a(I_n\partial_{x_j})=I_m df_a(\partial_{x_j})$ gives the displayed equations, and the equations from $\partial_{y_j}$ are the same system.
[/step]
[step:Derive the differential identity from holomorphicity]
Assume $F$ is holomorphic. Then for every choice of holomorphic charts as above, the coordinate representation
\begin{align*}
f: \varphi(U) \to \psi(V)
\end{align*}
is holomorphic. Hence each component $f_k = u_k + i v_k$ satisfies the several-variable Cauchy-Riemann equations at every $a \in \varphi(U)$. By the equivalence established in the previous step,
\begin{align*}
df_a \circ I_n = I_m \circ df_a
\end{align*}
for every $a \in \varphi(U)$. Taking $a = \varphi(p)$ and translating back through the tangent maps of the charts gives
\begin{align*}
dF_p \circ J_{X,p} = J_{Y,F(p)} \circ dF_p.
\end{align*}
Since $p \in X$ was arbitrary, the identity holds at every point of $X$.
[/step]
[step:Recover holomorphicity from complex linear differentials]
Assume that
\begin{align*}
dF_p \circ J_{X,p} = J_{Y,F(p)} \circ dF_p
\end{align*}
for every $p \in X$. Fix holomorphic charts $(U,\varphi)$ and $(V,\psi)$ with $F(U) \subset V$, and let
\begin{align*}
f: \varphi(U) \to \psi(V), \qquad f = \psi \circ F \circ \varphi^{-1}
\end{align*}
be the coordinate representation. The first step shows that
\begin{align*}
df_a \circ I_n = I_m \circ df_a
\end{align*}
for every $a \in \varphi(U)$. By the second step, the smooth real and imaginary component functions $u_k,v_k: \varphi(U) \to \mathbb{R}$ of each $f_k = u_k + i v_k$ satisfy the Cauchy-Riemann equations on $\varphi(U)$.
We now use the several-variable Cauchy-Riemann theorem: a smooth map from an open subset of $\mathbb{C}^n$ to $\mathbb{C}^m$ whose component functions satisfy the Cauchy-Riemann equations is holomorphic. The smoothness hypothesis is satisfied because $F$ and the coordinate charts are smooth. Therefore $f$ is holomorphic on $\varphi(U)$. Since this holds in holomorphic coordinate charts around every point $p \in X$, the map $F: X \to Y$ is holomorphic.
(citing a result not yet in the wiki: Several-variable Cauchy-Riemann theorem)
[/step]
[step:Conclude the equivalence]
The forward direction shows that every holomorphic smooth map $F: X \to Y$ satisfies
\begin{align*}
dF_p \circ J_{X,p} = J_{Y,F(p)} \circ dF_p
\end{align*}
for all $p \in X$. The converse shows that the same pointwise differential identity forces every holomorphic coordinate representation of $F$ to be holomorphic, hence $F$ is holomorphic as a map of complex manifolds. This proves the desired equivalence.
[/step]