An Expectation of Cohen-Lenstra Measure
Let me spell out the cokernel description of the Cohen-Lenstra distribution in more detail, as my answer will depend on it.
A map $(\mathbb{Z}_p)^N \to (\mathbb{Z}_p)^N$ is given by an $N \times N$ matrix of $p$-adic integers. Choose such a map by picking each of the digits of each integer uniformly at random from $\{ 0, 1, ..., p-1 \}$; this is the same as using the additive Haar measure on $\mathrm{Hom}((\mathbb{Z}_p)^N, (\mathbb{Z}_p)^N)$. With probability $1$, this map does not have determinant $0$, so its cokernel is a finite abelian $p$-group. Let $\mu_N$ be the probability measure on isomorphism classes of abelian $p$-groups which assigns each $p$-group the probability that it arises as this cokernel.
The Cohen-Lenstra distribution is the limit as $N \to \infty$ of $\mu_N$. As shown in several of the references you link to, it is given by the formula $$\lim_{N \to \infty} \mu_N(G) = \frac{1}{ |\mathrm{Aut}(G)|} \prod_{i=1}^{\infty} (1-1/p^i) .$$
For notational convenience, it will help to distinguish between the domain and range of a map in $\mathrm{Hom}((\mathbb{Z}_p)^N, (\mathbb{Z}_p)^N)$. I will call the former $U^N$ and the latter $V^N$.
Now, to answer your question. Let $A$ be a fixed finite abelian $p$-group. Let $e_N(A)$ be the expected number of surjections from an abelian $p$-group $G$, picked according to measure $\mu_N$, to $A$. Ignoring issues about interchanging limits, we want to show that $\lim_{N \to \infty} e_N(A)=1$.
Lets start by considering $H_N(A) := \mathrm{Hom}(V^N, A)$. The set $H_N(A)$ has cardinality $|A|^N$, as any map is specified by giving the image of a basis for $V^N$. Inside this set, let $S_N(A)$ be the surjective maps and $C_N(A)$ the nonsurjective maps.
For any map $f \in S_N(A)$, let's consider the possibility that it extends to the cokernel of a random map $U^N \to V^N$. This will occur if and only if the $N$ generators of $U^N$ land in the kernel of $f$. Since $f$ is in $S_N(A)$, its kernel has index $|A|$. So the probability that $U^N$ is mapped into the kernel of $f$ is $1/|A|^N$.
We want to compute $$e_N(A) = |S_N(A)| \cdot (1/|A|^N) = 1 - |C_N(A)|/A^N.$$
If $A$ can be generated by $r$ elements, then $|C_N(A)|/A^N \leq (1-1/|A|^r)^{\lfloor N/r \rfloor}$, so the second term drops out as $N \to \infty$. (To see this bound, group the basis elements of $V^N$ into $N/r$ groups of size $r$; the probability that these $r$ basis elements are not sent to the $r$ generators of $A$ is $(1-1/|A|^r)$. This bound is probably much weaker than the true rate of convergence.)