In Nielsen and Chuang (page:379), it is shown that the operator sum representation of a depolarizing channel $\mathcal{E}(\rho) = \frac{pI}{2} + (1-p)\rho$ is easily seen by substituting the identity matrix with
$$\frac{\mathbb{I}}{2} = \frac{\rho + X\rho X + Y\rho Y +Z\rho Z}{4}.$$
What is the more systematic way to see this result? Particularly, for the higher dimensional analogue, I cannot see how to proceed.
This really depends where you want to start from. For instance, you can construct the Choi state of $\mathcal E$, i.e., $$ \sigma = (\mathcal E \otimes \mathbb I)(|\Omega\rangle\langle\Omega|)\ , $$ with $\Omega = \tfrac{1}{\sqrt{D}}\sum_{i=1}^D |i,i\rangle$, and then extract the Kraus operators of $\mathcal E(\rho)=\sum M_i\rho M_i^\dagger$ by taking any decomposition $$ \sigma = \sum |\psi_i\rangle\langle\psi_i|\ ,\tag{*} $$ and writing $|\psi_i\rangle = (M_i\otimes\mathbb I)|\Omega\rangle$ (which is always possible).
Note that the decomposition $(*)$ is highly non-unique (any $|\phi_j\rangle = \sum V_{ij} |\psi_i\rangle$, with $V$ an isometry, is also a valid decomposition), which relates to the fact that the Kraus decomposition is equally non-unique. Obviously, the eigenvalue decomposition is a simple choice (which, moreover, minimizes the number of Kraus operators).
Let's look at your example in a bit more detail. Here, $D=2$. You have that $$ \mathcal E(X)=p\mathrm{tr}(X)\,\frac{\mathbb I}{2}+(1-p)X $$ for any $X$ (due to linearity) -- the $\mathrm{tr}(X)$ is required to make this trace-preserving for general $X$.
We now have that \begin{align} \sigma &= (\mathcal E \otimes \mathbb I)(|\Omega\rangle\langle \Omega|) \\ & = \tfrac1D \sum_{ij} \mathcal E(|i\rangle\langle j|)\otimes |i\rangle\langle j|\ \end{align} inserting the definition of $|\Omega\rangle$ and using linearity.
This yields $$ \sigma = \frac{p}{2D}\mathbb I\otimes \sum_{i}|i\rangle\langle i| + (1-p)\frac1D \sum_{ij}|i\rangle\langle j|\otimes |i\rangle\langle j|\ . $$ The second term is just $(1-p)|\Omega\rangle\langle\Omega|$, and the first term is $\frac{p}{2D}\mathbb I\otimes\mathbb I$.
You can now see that one possible eigenvalue decomposition of $\sigma$ is given by the four Bell states (I leave it to you to work out the weights), and it is well known and easy to check that that the four Bell states can be written as $$ (\sigma_k\otimes \mathbb I)|\Omega\rangle\ , $$ where $\sigma_k$ are the three Pauli matrices or the identity.
Thus, you get that the $M_i$ in the Kraus representation are the Paulis and the identity, with the weight given by the eigenvalue decomposition of $\sigma$.
While the procedure in the existing answer, based on channel-state duality, applies to general channels, there's a more direct way to obtain Kraus operators for this particular case of the depolarising channel, $\mathcal{E}(\rho) = (1-p)\rho + p \frac{I}{d}$, where $d$ is the Hilbert space dimension. Since $\mathrm{tr}(\rho) = 1$ for any density operator, we can write $\frac{I}{d} = \frac{I}{d}\mathrm{tr}(\rho)$. For any orthonormal basis $\{|1\rangle, \dots, |d\rangle\}$, $I = \sum_i |i\rangle \langle i|$ and $\mathrm{tr}(\rho) = \sum_i \langle i|\rho |i\rangle$, so $$ \frac{I}{d} = \frac{1}{d}\sum_i|i\rangle\langle i| \sum_j \langle j|\rho|j\rangle = \frac{1}{d} \sum_{i,j} |i\rangle \langle j|\rho|j \rangle\langle i|,$$ from which we see that a possible set of Kraus operators for $\mathcal{E}$ are $\{|i\rangle \langle j|/\sqrt{d}\}$, with $i,j \in \{1,\dots, d\}$. (You can check that these satisfy the completeness relation.)
For the one-qubit case, note that $\{I/2,X/2,Y/2,Z/2\}$ are related to $\{|0\rangle\langle 0|/\sqrt{2}, |0\rangle\langle 1|/\sqrt{2}, |1\rangle\langle 0|/\sqrt{2},|1\rangle\langle 1|/\sqrt{2}\}$ by a unitary, so by using the unitary freedom in the Kraus representation, we can get from the above representation to the standard $\frac{I}{2} = \frac{\rho + X\rho X + Y\rho Y + Z\rho Z}{4}$.
