Write the dephasing channel, acting on arbitrary $d$-dimensional states, as
$$D(\rho)\equiv D_p(\rho) = (1-p)\rho + p \sum_k\operatorname{Tr}(E_{kk}\rho)E_{kk}.$$
Using the notation $E_{ij}\equiv|i\rangle\!\langle j|$.
(1) Compute the Choi representation
It follows that $D(E_{ij})=(1-p)E_{ij}$ for $i\neq j$, and $D(E_{ii})=
E_{ii}$.
Denote the Choi representation of $D$ as $J(D)$. One way to define the Choi of a channel is as
$$J(D)\equiv \sum_{ij} D(E_{ij})\otimes E_{ij},$$
which using the above rules immediately becomes
$$J(D) = (1-p)\sum_{i\neq j} E_{ij}\otimes E_{ij}
+ \sum_i E_{ii}\otimes E_{ii} \\
= \sum_i E_{ii}\otimes E_{ii} + (1-p)\sum_{i\neq j}E_{ij}\otimes E_{ij} \\
\equiv \sum_i |ii\rangle\!\langle ii| + (1-p)\sum_{i\neq j}|ii\rangle\!\langle jj|.$$
The same result could also equivalently be obtained observing that $$D_p=(1-p)\operatorname{Id} + p \operatorname{FullyDephasingCh},$$
that the map $D\mapsto J(D)$ is linear, and that $J(\operatorname{Id})=|m\rangle\!\langle m|$ with $|m\rangle\equiv\sum_i |ii\rangle$, and
$$
J(\operatorname{FullyDephasingCh}) = \sum_i |ii\rangle\!\langle ii|.
$$
(2) Eigendecomposition of the Choi
A straightforward way to get a Kraus decomposition is via the eigendecomposition of the Choi.
To this end, start by observing that support and image of $J(D)$ are spanned by the vectors $\{|ii\rangle : i=1,...,d\}$. We can thus make a formal simplification replacing $|ii\rangle\to|i\rangle$. Upon this replacement, $J(D_p)$ gets a nice matrix representation:
$$J(D_p)=I + (1-p)(\boldsymbol1-I),$$
where $I$ is the identity matrix, and $\boldsymbol1$ the constant matrix with $(\boldsymbol1)_{ij}=1$ for all $i,j$.
For concreteness, in the $d=2$ case the above amounts to
$$J(D_p) = \begin{pmatrix}1 & 1-p \\ 1-p & 1\end{pmatrix}.$$
It is not too hard to verify that the eigenvalues of $J(D_p)$ are a nondegenerate $\lambda_1=d(1-p)+p$, and a $(d-1)$-fold degenerate $\lambda_2=p$. A corresponding set of eigenvectors is
$$|\lambda_1\rangle = \frac{1}{\sqrt d}\sum_i |i\rangle, \\
|\lambda_2^j\rangle = \frac{1}{\sqrt2} (|j\rangle - |1\rangle),
\qquad j=2,...,d.$$
The $|\lambda_2^j\rangle$ here are eigenvectors corresponding to $\lambda_2$, but note that they are not orthogonal.
An orthonormal system of eigenvectors can be written in the form
$$|\lambda_2^j\rangle =
\frac{1}{\sqrt{j(j+1)}}\left(\sum_{i=1}^j|i\rangle - j|j+1\rangle\right),
\qquad j=1,...,d-1.
$$
These correspond to the following spectral decomposition for the Choi:
$$J(D_p) = ((1-p)d+p)\frac{\boldsymbol1}{d} + p\left(I-\frac{\boldsymbol1}{d}\right),$$
where we observe that $\frac{\boldsymbol1}{d}$ is a unit-trace orthoprojection, and $\left(I-\frac{\boldsymbol1}{d}\right)$ an orthogonal projection with trace $d-1$, as expected.
(3) Kraus operators from Choi's eigendecomposition
Given a map $\Phi\in\mathrm T(\mathcal X,\mathcal Y)$ we have $$\Phi(X)=\sum_{a}A_a X B_a^\dagger \iff
J(\Phi) = \sum_a \operatorname{vec}(A_a)\operatorname{vec}(B_b)^\dagger,
$$
for some linear operators $A_a,B_b:\mathcal X\to\mathcal Y$. Also, one can always find such a representation, with operators such that $\langle A_a,A_b\rangle=\langle B_a,B_b\rangle=0$ for any $a\neq b$.
In particular, if the map $\Phi$ is CP, and thus its Choi is positive semidefinite, we get the "standard" Kraus decomposition, and
$$\Phi(X)=\sum_a A_a X A_a^\dagger
\iff
J(\Phi) = \sum_a \operatorname{vec}(A_a)\operatorname{vec}(A_a)^\dagger.
$$
This means that the eigenvectors of the Choi are tightly related to (one choice of) Kraus operators. More specifically, if we have an eigendecomposition of the form
$$J(\Phi) = \sum_a \lambda_a \mathbb P(v_a),$$
with $v_a\in\mathcal Y\otimes\mathcal X$ such that $\|v_a\|=1$, and $\mathbb P(v_a)$ denoting the corresponding rank-1 projection, then a valid choice of Kraus operators is
$$A_a = \sqrt{\lambda_a} \operatorname{unvec}(v_a) \in\mathrm{Lin}(\mathcal X,\mathcal Y).$$
Going back to the case at hand, using the eigenvalues and orthonormal eigenvectors previously derived, we get the Kraus operators
$$
A_1 = \frac{\sqrt{d+p(1-d)}}{\sqrt d} \sum_i E_{ii} \equiv \frac{1}{\sqrt d}I,
\qquad
A_2^k = \frac{\sqrt{p}}{\sqrt{k(k+1)}} \left(
\sum_{j=1}^k E_{jj} - k E_{k+1}
\right),
$$
for $k=1,...,d-1$.
(4) Examples
For $d=2$, the above amount to
$$A_1 = \frac{\sqrt{2-p}}{\sqrt2} I,
\qquad
A_2 = \frac{\sqrt{p}}{\sqrt2} Z.
$$
For $d=3$, we get instead
$$
A_1 = \frac{\sqrt{3-2p}}{\sqrt3} I,
\quad
A_2 = \frac{\sqrt{p}}{\sqrt2} (E_{11} - E_{22}),
\quad
A_3 = \frac{\sqrt{p}}{\sqrt6} (E_{11} + E_{22} - 2 E_{33}).
$$
