If you consider a truncated Poisson (i.e., getting rid of 0), then Mathematica can do all of the heavy lifting.
First we find the distribution of a truncated Poisson:
d1 = TruncatedDistribution[{1, \[Infinity]}, PoissonDistribution[\[Lambda]]]
For the sum of
$n$ independent normal random variable with mean
$\mu$ and variance
$\sigma^2$, the resulting mean and variance of that sum are
$n \mu$ and
$n \sigma^2$. We use the ParameterMixtureDistribution function to find the distribution of the sum of a random number of normal random variables:
d2 = ParameterMixtureDistribution[NormalDistribution[n \[Mu], Sqrt[n] \[Sigma]], n \[Distributed] d1]
The mean and variance are found with the following:
mean = Mean[d2] // FunctionExpand // FullSimplify
$$\frac{\left(e^{\lambda }-1\right) \lambda \mu }{-\lambda +e^{\lambda }-1}$$
variance = Variance[d2] // FunctionExpand // FullSimplify
$$\frac{\lambda \left(\left(e^{\lambda }-1\right) \left(-\lambda +e^{\lambda }-1\right) \sigma ^2-e^{\lambda } \mu ^2 \left(\lambda ^2-2 \cosh (\lambda )+2\right)\right)}{\left(\lambda -e^{\lambda }+1\right)^2}$$
If
$\lambda$ is large (such the the probability of zero is very low), then one obtains the following for
$\lambda=15$:
mean /. \[Lambda] -> 15 // N
(* 15.0001 \[Mu] *)
variance /. \[Lambda] -> 15 // N // Expand
(* 14.9991 \[Mu]^2 + 15.0001 \[Sigma]^2 *)
