Here is a proof by local calculation, but I am not sure whether it's right.
WLOG we may take $p = (0,\cdots,0,1)\in S^n\subset\mathbb{R}^{n+1}$, then we have a local coordinate $(U;z_1,\cdots,z_n)$ near $p$, where $U=\{x\in \mathbb{R}^{n+1}\mid x_{n+1}>0 \}$, and we have a homeomorphism(projection)
\begin{align*}
\varphi: U &\longrightarrow \mathbb{R}^n \\
(x_1,\cdots,x_{n+1})&\longmapsto (x_1,\cdots,x_n)
\end{align*}
The local expression of $Y$ is $Y\circ\varphi^{-1} =: \tilde{Y}$, hence
$$ dY(X) = \sum_{i=1}^{n+1}\left(\sum_{j=1}^n \partial_j\tilde{Y}^i\cdot {X}^j\right)\partial_i $$
where $\partial_i := \partial/\partial z_i,\ i=1,2,\cdots,n$, $\partial_{n+1}$ is defined to be the vector perpendicular to $T_pS^n = \text{Span}\{\partial_1,\cdots,\partial_n\}$ and preserves the orientation, $\tilde{Y} =\sum_{i=1}^{n+1} \tilde{Y}^i\partial_i$, $X = \sum_{j=1}^n X^j\partial_j$.
Similarly we have
$$ dX(Y) = \sum_{i=1}^{n+1}\left(\sum_{j=1}^n \partial_j\tilde{X}^i\cdot {Y}^j\right)\partial_i $$
Note that $\tilde{X}^i_p = X^i_p$, $\tilde{Y}^i_p = Y^i_p$ for $i=1,2,\cdots,n$ and $\tilde{X}^{n+1} = \tilde{Y}^{n+1} = 0$, so
$$ \left(dY(X) - dX(Y)\right)_p = \sum_{i=1}^n\sum_{j=1}^n\left( X^j\partial_jY^i - Y^j\partial_jX^i\right)\partial_i = [X, Y]_p $$
where the middle term takes value at $p$.