The channel representation of a quantum gate $U$ essentially corresponds to the action of a gate on the density matrix $\rho$, as in $\mathcal{S}(U) = U \rho U^{\dagger}$, instead of its action on the pure state wave function $|\psi \rangle$, as in $U |\psi \rangle$. The conventional description of quantum circuits, particularly in the absence of mid-circuit measurements, is the latter, as we can regard the state of the quantum register as a pure state $|\psi\rangle$ acted upon by unitary operators $\{U_i\}$ (i.e., quantum gates), which maintain the purity of the state. The channel representation is the appropriate description of a circuit where some of its operations are non-unitary (e.g., mid-circuit measurements), in which case the state of the quantum register can only be fully described in terms of a density operator $\rho$ acted upon by superoperators $\mathcal{S}(U)$.
On page 3 of the paper, Mitarai and Fujii explain why they adopted this superoperator decomposition instead of the standard operator decomposition.
It is noteworthy that as we perform decompositions of a superoperator rather than an operator such as $U$ itself, the method becomes friendly for a realistic quantum device. A direct decomposition of $U$ into some simple operators {$V_i$}, i.e. $U = \sum_i c_i V_i$, can also be utilized for the same task; however, as expectation values are calculated as $\langle 0|U^{\dagger} O U|0\rangle$ where $|0\rangle$ is an initial state, this approach requires us to evaluate $\sum_{i,j} c_i c^{*}_j \langle 0|V_{j}^{\dagger} O V_i|0\rangle$, which are rather hard for the NISQ devices. This fact demonstrates the advantage of using the above formalism. The tensor network representation of the superoperator formalism allows us to graphically understand the decompositions.
To rephrase the quote, it is advantageous to use the superoperator formalism in this context because it allows us to express the expectation value that we wish to compute as a linear combination of expectation values with the symmetric form $\langle 0|P^{\dagger} O P|0\rangle$, while the operator formalism would yield some expectation values $\langle 0|Q^{\dagger} O P|0\rangle$ with $P \neq Q$. The former can be straightforwardly computed, while the latter typically require more sophisticated methods that roughly double the circuit depth (cf., e.g., Fig. 9 in R. Somma et al., Phys. Rev. A 65, 042323 (2002))
The point about this paper by Mitarai and Fujii is not to express a two-qubit gate as a linear combination of tensor products of single-qubit gates. In fact, for the sort of two-qubit gate they consider, $e^{i \theta A_1 \otimes A_2}$, with $A_1^2 = \mathbb{1}$ and $A_2^2 = \mathbb{1}$, there is one such trivial decomposition: $e^{i \theta A_1 \otimes A_2} = \cos \theta \; \mathbb{1} \otimes \mathbb{1} + i \sin \theta \; A_1 \otimes A_2$. Hence, as you rightly pointed out in your edit, the equalities in Fig. 1 of the paper do not refer to the $4 \times 4$ matrix representations of the two-qubit operations but rather to the expectation values computed with such two-qubit operations. Or, equivalently, in the words of one of the authors, the equality holds when we consider the operations as quantum channels.
Regarding your ultimate goal of obtaining a quasiprobabilistic decomposition of the CNOT gate, it is straightforward to adapt the expression provided by Mitarai and Fujii for the CZ gate using the identity $CNOT = (\mathbb{1} \otimes H) CZ (\mathbb{1} \otimes H)$, where $H$ is the Hadamard gate. The analogue of Fig. 2 for the CNOT gate is as follows.
Basically, all the Z operations are replaced by X operations in the target qubit only. You can easily check that this linear combination yields exactly the same expectation values as the original CNOT gate. However, its $4 \times 4$ matrix representation is
$\begin{pmatrix}
1 & -i & 0 & 0 \\
-i & 1 & 0 & 0 \\
0 & 0 & 0 & 1+i \\
0 & 0 & 1+i & 0
\end{pmatrix}$,
which is clearly not the CNOT gate.
