Discussing with @jgon we got to the same conclusion, therefore I assume we have the right answer which I'm posting just in case someone else needs in the future.
In order to find the correct expression, it suffices to find $\sigma$ in:
$$\varepsilon(a_1, \ldots, a_{j_1-1}, b_{j_1}, a_{j_1+1}, \ldots, a_{j_i-1}, b_{j_i}, a_{j_i+1}, \ldots, a_p)=(-1)^\sigma \varepsilon(a_1, \ldots, \widehat{a_{j_1}}, \ldots, \widehat{a_{j_i}}, \ldots, a_p, b_{j_1}, \ldots, b_{j_i}).$$
Let us focus on $b_{j_1}$. For it to appear beside $a_p$ on the RHS of the above equality, it will have to go past $p-j_{1}$ positions. Therefore:
$$\varepsilon(a_1, \ldots, a_{j_1-1}, a_{j_1+1}, \ldots, a_{j_i-1}, b_{j_i}, a_{j_i+1}, \ldots, a_p)\\
= (-1)^{p-j_1} \varepsilon(a_1, \ldots, \widehat{a_{j_1}}, \ldots, a_{j_i-1}, b_{j_i}, a_{j_i+1}, \ldots, a_p, b_{j_1}).$$
Now we focus on $b_{j_2}$. It will have to go past $p-j_2+1$ places, therefore:
$$\varepsilon(a_1, \ldots, a_{j_1-1}, a_{j_1+1}, \ldots, a_{j_2-1}, b_{j_2}, a_{j_2+1}, \ldots, a_{j_i-1}, b_{j_i}, a_{j_i+1}, \ldots, a_p)\\= (-1)^{p-j_1} (-1)^{p-j_2+1} \varepsilon(a_1, \ldots, \widehat{a_{j_1}}, \ldots, \widehat{a_{j_2}}, \ldots, a_{j_i-1}, b_{j_i}, a_{j_i+1}, \ldots, a_p, b_{j_1}, b_{j_2}). $$
The last $b_{j_i}$ will have to go past $p-j_i+(1+\ldots+(i-1))$ places, therefore:
$$\varepsilon(a_1, \ldots, a_{j_1-1}, a_{j_1+1}, \ldots, a_{j_i-1},b_{j_i}, a_{j_i+1}, \ldots, a_p)=(-1)^{p-j_1}(-1)^{p-j_2+1}\cdots (-1)^{p-j_i+(i-1)}\varepsilon(a_1, \ldots, \widehat{a_{j_1}}, \ldots, \widehat{a_{j_i}}, \ldots, a_p, b_{j_1}, \ldots, b_{j_i}).$$ Now notice:
$$(-1)^{p-j_1}(-1)^{p-j_2+1}\cdots (-1)^{p-j_i+(i-1)}=(-1)^{pi-(j_1+\ldots+j_i)+ \frac{i(i-1)}{2}}$$ where I used the standard fact $$1+\ldots+n=\frac{n(n+1)}{2}$$ with $n=i-1$. Therefore:
$$\varepsilon(a_1+b_1, \ldots, a_p+b_p)= \varepsilon(a_1, \ldots, a_p)\\
+ \sum_{i=1}^{p-1} (-1)^{ip} (-1)^{i(i-1)/2} \sum_{1\leq j_1<\ldots<j_i\leq p} (-1)^{j_1+\ldots+j_i} \varepsilon(a_1, \ldots, \widehat{a_{j_1}}, \ldots, \widehat{a_{j_i}}, \ldots, a_p, b_{j_1}, \ldots, b_{j_p})\\
+\varepsilon(b_1, \ldots, b_p).$$ I know this can be written in a straightforward manner using shuffles, but I wanted to avoid that.