No, you diagonalize in the truncated degeneracy subspace D, so your eigenkets are definitely not eigenkets of the full hamiltonian!
Most texts are needlessly formal and lose you in a thicket of abstract formalism, and you never thought to illustrate them by a daft minimal example. Here is one.
Take
$$
H_0+\lambda V= \begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 2
\end{pmatrix} +\lambda \begin{pmatrix}
0 & 1 & 2 \\
1 & 0 & 3\\
2 & 3 & 7
\end{pmatrix} .$$
Consider
$$|1\rangle=\begin{pmatrix}
1 \\
0 \\
0
\end{pmatrix} ,\qquad |2\rangle=\begin{pmatrix}
0 \\
1 \\
0
\end{pmatrix} , \qquad |3\rangle=\begin{pmatrix}
0 \\
0 \\
1
\end{pmatrix} ,$$
$$ |a\rangle= \frac{1}{\sqrt{2}} \begin{pmatrix}
1 \\
-1 \\
0
\end{pmatrix} , \qquad |b\rangle= \frac{1}{\sqrt{2}}\begin{pmatrix}
1 \\
1 \\
0
\end{pmatrix} ,$$
so the orthogonal states $|a\rangle, |b\rangle$ diagonalize the 1-2 subspace perturbation, but obviously not the entire hamiltonian.
You must check by direct calculation, and also pursuant to your text or WP, that, to lowest order in λ, that is ignoring its square,
the energy eigenvalues and eigenstates are
$$
E_A= 1-\lambda; \qquad E_B= 1+\lambda; \qquad E_C=2+\lambda 7 ;\\
|A\rangle= |a\rangle + \frac{\lambda}{\sqrt{2}} |3\rangle= \frac{1}{\sqrt{2}}\begin{pmatrix}
1 \\
- 1\\ \lambda
\end{pmatrix} ~; \qquad |B\rangle= |b\rangle -\lambda \frac{5} {\sqrt{2}}|3\rangle= \frac{1}{\sqrt{2}}\begin{pmatrix}
1 \\ 1\\
-5\lambda
\end{pmatrix} ~; \\
|C\rangle= |3\rangle -\lambda \frac{1} {\sqrt{2}}|a\rangle +\lambda \frac{5} {\sqrt{2}}|b\rangle= \begin{pmatrix}
2\lambda \\ 3\lambda\\
1
\end{pmatrix} ~.
$$
Recall normalizations, orthogonalities, etc... all work to $O(\lambda^2)$ only.
In conclusion, $|a\rangle, |b\rangle$ are not exact hamiltonian eigenstates--and you have all perturbation orders survive, in principle. They merely avoided zeros in the denominator of the resolvent, having removed the degeneracy. That is, they ensured that $|A\rangle$ will have no $|b\rangle$ correction to first order, and, likewise, $|B\rangle$ will have no $|a\rangle$ correction, because this basis has diagonalized V, so it fails to connect $|a\rangle$ with $|b\rangle$ and they do not impact each other, so far. But higher orders are another matter.
- Fine print (wonkish) footnote on the above "so far": avoid until 3rd reading.
In actuality, an extra arbitrary piece $\propto \lambda |b\rangle$ in $|A\rangle$ and $\propto \lambda |a\rangle$ in $|B\rangle$ has no impact (linkage) on the eigenvalue equation at the first order in λ discussed above, as you may readily check: the O(λ) potential has superselected them out to this order. However, as you may check in Courant-Hilbert Methods of Mathematical Physics, v1 p 248 , and discussed in this question , you do get two innocuous (impact-less) pieces to this order, namely
$$
|A\rangle=\hbox{ above } - \lambda \frac{5}{4} |b\rangle,\qquad
|B\rangle=\hbox{ above } + \lambda \frac{5}{4} |a\rangle,\\
\langle b|A\rangle=\frac{\lambda}{2}\langle b|V|3\rangle \frac{1}{E^{(0)}_{a,b}-E^{(0)}_{3}}\langle 3|V|a\rangle=-5/4, \qquad \langle a|B\rangle=5/4.
$$
These pieces, bilinear in the potential!, structurally "fell off" the 2nd order calculation, where λ2 was discounted by a λ in a subspace denominator $ E^{(1)}_{b}-E^{(1)}_{a}=2\lambda$. The particular 5/4 coefficients are fixed by consistency with the next order energies. But you have already seen this subtlety above, in zeroth order: $|a\rangle, |b\rangle$ of $O(\lambda^0)$ actually were specified by $O(\lambda)$ potential considerations. In any case, for reassurance, compute
$$
E_A= 1-\lambda -\lambda^2/2 + 7\lambda^3/8+...; \\ E_B= 1+\lambda-25\lambda^2/2+625\lambda^3/8+...; \\ E_C=2+7\lambda +13\lambda^2-79\lambda^3+...
$$
