[proofplan]
We prove the identity pointwise by evaluating both sides on arbitrary tangent vectors. The wedge product is expanded using the shuffle formula, and the shuffle sum is split according to whether the inserted vector $V_p$ is assigned to the $\alpha$-slot or to the $\beta$-slot. The first class of terms gives $(\iota_V\alpha)\wedge\beta$, while the second class gives $\alpha\wedge(\iota_V\beta)$ with the sign $(-1)^k$ coming from moving $V_p$ past the $k$ arguments assigned to $\alpha$.
[/proofplan]
[step:Reduce the identity to a pointwise alternating multilinear computation]
Fix a point $p\in M$. Let
\begin{align*}
X_1,\dots,X_{k+\ell-1}\in T_pM
\end{align*}
be arbitrary tangent vectors. Define tangent vectors
\begin{align*}
Y_0 &:= V_p, &
Y_j &:= X_j \quad \text{for } 1\leq j\leq k+\ell-1.
\end{align*}
For a form $\omega\in\Omega^r(M)$ with $r\geq 1$, the interior product is defined pointwise by
\begin{align*}
(\iota_V\omega)_p(Z_1,\dots,Z_{r-1})
=
\omega_p(V_p,Z_1,\dots,Z_{r-1}).
\end{align*}
Thus it is enough to prove
\begin{align*}
(\alpha\wedge\beta)_p(Y_0,Y_1,\dots,Y_{k+\ell-1})
&=
\bigl((\iota_V\alpha)\wedge\beta\bigr)_p(X_1,\dots,X_{k+\ell-1})\\
&\quad+
(-1)^k\bigl(\alpha\wedge(\iota_V\beta)\bigr)_p(X_1,\dots,X_{k+\ell-1}).
\end{align*}
Since the point $p$ and the tangent vectors $X_1,\dots,X_{k+\ell-1}$ are arbitrary, this pointwise identity proves the equality of forms.
[/step]
[step:Expand the wedge product by shuffles and isolate the terms where $V_p$ enters $\alpha$]
Let $\operatorname{Sh}(k,\ell)$ denote the set of $(k,\ell)$-shuffles of $\{0,1,\dots,k+\ell-1\}$, that is, permutations $\sigma$ such that
\begin{align*}
\sigma(1)<\cdots<\sigma(k),
\qquad
\sigma(k+1)<\cdots<\sigma(k+\ell).
\end{align*}
For a permutation $\sigma$ of this finite ordered set, let $\operatorname{sgn}(\sigma)\in\{1,-1\}$ denote its permutation sign, defined as $1$ for even permutations and $-1$ for odd permutations.
The shuffle formula gives
\begin{align*}
(\alpha\wedge\beta)_p(Y_0,\dots,Y_{k+\ell-1})
=
\sum_{\sigma\in\operatorname{Sh}(k,\ell)}
\operatorname{sgn}(\sigma)\,
\alpha_p(Y_{\sigma(1)},\dots,Y_{\sigma(k)})
\beta_p(Y_{\sigma(k+1)},\dots,Y_{\sigma(k+\ell)}).
\end{align*}
First consider the shuffles for which $0$ occurs among $\sigma(1),\dots,\sigma(k)$. Since these $k$ entries are increasing and $0$ is the least index, this means $\sigma(1)=0$. Removing the initial entry $0$ gives a bijection from this class of shuffles to $\operatorname{Sh}(k-1,\ell)$ on the index set $\{1,\dots,k+\ell-1\}$. Under this bijection the sign is unchanged, because $0$ already occupies the first position. Therefore this part of the sum equals
\begin{align*}
\sum_{\tau\in\operatorname{Sh}(k-1,\ell)}
\operatorname{sgn}(\tau)\,
\alpha_p(V_p,X_{\tau(1)},\dots,X_{\tau(k-1)})
\beta_p(X_{\tau(k)},\dots,X_{\tau(k+\ell-1)}).
\end{align*}
By the definition of $\iota_V\alpha$, this is precisely
\begin{align*}
\bigl((\iota_V\alpha)\wedge\beta\bigr)_p(X_1,\dots,X_{k+\ell-1}).
\end{align*}
[guided]
We now expand the left-hand side in the only place where a sign can enter: the antisymmetrised definition of the wedge product. Let $\operatorname{Sh}(k,\ell)$ be the set of $(k,\ell)$-shuffles of the ordered index set $\{0,1,\dots,k+\ell-1\}$. Thus $\sigma\in\operatorname{Sh}(k,\ell)$ means that the first $k$ selected indices remain increasing and the last $\ell$ selected indices remain increasing:
\begin{align*}
\sigma(1)<\cdots<\sigma(k),
\qquad
\sigma(k+1)<\cdots<\sigma(k+\ell).
\end{align*}
For a permutation $\sigma$ of this finite ordered set, let $\operatorname{sgn}(\sigma)\in\{1,-1\}$ denote its permutation sign: $\operatorname{sgn}(\sigma)=1$ when $\sigma$ is even and $\operatorname{sgn}(\sigma)=-1$ when $\sigma$ is odd.
The wedge product is therefore
\begin{align*}
(\alpha\wedge\beta)_p(Y_0,\dots,Y_{k+\ell-1})
=
\sum_{\sigma\in\operatorname{Sh}(k,\ell)}
\operatorname{sgn}(\sigma)\,
\alpha_p(Y_{\sigma(1)},\dots,Y_{\sigma(k)})
\beta_p(Y_{\sigma(k+1)},\dots,Y_{\sigma(k+\ell)}).
\end{align*}
We split this sum into two pieces. In the first piece, the inserted vector $Y_0=V_p$ is one of the arguments of $\alpha_p$. Since the $\alpha$-indices are increasing and $0$ is the smallest possible index, this forces $\sigma(1)=0$. After deleting this first entry, the remaining indices form a $(k-1,\ell)$-shuffle of $\{1,\dots,k+\ell-1\}$. The deletion does not change the sign, because the index $0$ was already in the first position. Hence this first piece is
\begin{align*}
\sum_{\tau\in\operatorname{Sh}(k-1,\ell)}
\operatorname{sgn}(\tau)\,
\alpha_p(V_p,X_{\tau(1)},\dots,X_{\tau(k-1)})
\beta_p(X_{\tau(k)},\dots,X_{\tau(k+\ell-1)}).
\end{align*}
By the definition of interior product,
\begin{align*}
\alpha_p(V_p,X_{\tau(1)},\dots,X_{\tau(k-1)})
=
(\iota_V\alpha)_p(X_{\tau(1)},\dots,X_{\tau(k-1)}).
\end{align*}
Substituting this into the preceding shuffle sum gives exactly the shuffle formula for
\begin{align*}
\bigl((\iota_V\alpha)\wedge\beta\bigr)_p(X_1,\dots,X_{k+\ell-1}).
\end{align*}
[/guided]
[/step]
[step:Identify the terms where $V_p$ enters $\beta$ and compute the sign]
Now consider the shuffles for which $0$ occurs among $\sigma(k+1),\dots,\sigma(k+\ell)$. Since those entries are increasing, this means $\sigma(k+1)=0$. To compare this with the shuffle formula for $\alpha\wedge(\iota_V\beta)$, move the entry $0$ from position $k+1$ to the first $\beta$-argument position after the $k$ arguments of $\alpha$ have been chosen from $X_1,\dots,X_{k+\ell-1}$. Equivalently, $0$ has crossed exactly the $k$ arguments assigned to $\alpha$. Therefore the sign changes by the factor $(-1)^k$.
More explicitly, deleting the entry $0$ gives a bijection from this class of shuffles to $\operatorname{Sh}(k,\ell-1)$ on $\{1,\dots,k+\ell-1\}$, and for the corresponding shuffle $\tau$ one has
\begin{align*}
\operatorname{sgn}(\sigma)=(-1)^k\operatorname{sgn}(\tau).
\end{align*}
Thus the second part of the shuffle sum equals
\begin{align*}
(-1)^k
\sum_{\tau\in\operatorname{Sh}(k,\ell-1)}
\operatorname{sgn}(\tau)\,
\alpha_p(X_{\tau(1)},\dots,X_{\tau(k)})
\beta_p(V_p,X_{\tau(k+1)},\dots,X_{\tau(k+\ell-1)}).
\end{align*}
By the definition of $\iota_V\beta$, this is
\begin{align*}
(-1)^k
\bigl(\alpha\wedge(\iota_V\beta)\bigr)_p(X_1,\dots,X_{k+\ell-1}).
\end{align*}
[guided]
The second piece of the shuffle sum consists of the terms in which the inserted vector $Y_0=V_p$ is assigned to $\beta_p$. Since the $\beta$-indices are increasing, the smallest index $0$ must appear first among the $\beta$-indices, so $\sigma(k+1)=0$.
We now compare this with the shuffle formula for $\alpha\wedge(\iota_V\beta)$. In that formula, the form $\alpha$ first receives $k$ of the vectors $X_1,\dots,X_{k+\ell-1}$, and then $\iota_V\beta$ receives the remaining $\ell-1$ vectors. Since
\begin{align*}
(\iota_V\beta)_p(Z_1,\dots,Z_{\ell-1})
=
\beta_p(V_p,Z_1,\dots,Z_{\ell-1}),
\end{align*}
the vector $V_p$ is placed before the remaining $\beta$-arguments.
The only issue is the sign. Starting from the shuffle term for $\alpha\wedge\beta$, the entry $0$ lies after the $k$ entries assigned to $\alpha$. To rewrite the term as one involving $\iota_V\beta$, we move this entry past exactly those $k$ entries. Each transposition changes the sign by $-1$, so the total sign factor is $(-1)^k$. Equivalently, deleting $0$ identifies this class of shuffles with $\operatorname{Sh}(k,\ell-1)$ on $\{1,\dots,k+\ell-1\}$, and the corresponding signs satisfy
\begin{align*}
\operatorname{sgn}(\sigma)=(-1)^k\operatorname{sgn}(\tau).
\end{align*}
Therefore this second class of terms is
\begin{align*}
(-1)^k
\sum_{\tau\in\operatorname{Sh}(k,\ell-1)}
\operatorname{sgn}(\tau)\,
\alpha_p(X_{\tau(1)},\dots,X_{\tau(k)})
\beta_p(V_p,X_{\tau(k+1)},\dots,X_{\tau(k+\ell-1)}).
\end{align*}
Using the definition of $\iota_V\beta$, this becomes
\begin{align*}
(-1)^k
\sum_{\tau\in\operatorname{Sh}(k,\ell-1)}
\operatorname{sgn}(\tau)\,
\alpha_p(X_{\tau(1)},\dots,X_{\tau(k)})
(\iota_V\beta)_p(X_{\tau(k+1)},\dots,X_{\tau(k+\ell-1)}),
\end{align*}
which is exactly
\begin{align*}
(-1)^k
\bigl(\alpha\wedge(\iota_V\beta)\bigr)_p(X_1,\dots,X_{k+\ell-1}).
\end{align*}
[/guided]
[/step]
[step:Combine the two shuffle classes and handle degree zero cases]
The two classes of shuffles above are disjoint and exhaustive: the index $0$ is assigned either to the $\alpha$-block or to the $\beta$-block. Adding their contributions gives
\begin{align*}
(\iota_V(\alpha\wedge\beta))_p(X_1,\dots,X_{k+\ell-1})
&=
\bigl((\iota_V\alpha)\wedge\beta\bigr)_p(X_1,\dots,X_{k+\ell-1})\\
&\quad+
(-1)^k\bigl(\alpha\wedge(\iota_V\beta)\bigr)_p(X_1,\dots,X_{k+\ell-1}).
\end{align*}
If $k=0$, then $\alpha$ is a smooth function and $\iota_V\alpha=0$ by convention; the first class of shuffles is empty and the same computation gives
\begin{align*}
\iota_V(\alpha\wedge\beta)=\alpha\wedge(\iota_V\beta),
\end{align*}
which is the stated identity because $(-1)^0=1$. If $\ell=0$, then $\iota_V\beta=0$ by convention; the second class is empty and the identity reduces to
\begin{align*}
\iota_V(\alpha\wedge\beta)=(\iota_V\alpha)\wedge\beta.
\end{align*}
Thus the formula holds for all $k,\ell\in\mathbb{N}_0$.
[/step]
