I can partially answer your question: counting the local optima of a PLS-complete search problem can indeed be #P-hard.
First, as Yoshio points out, there is a search problem $P_1$ in PLS whose associated counting problem is #P-complete. (We don't know if $P_1$ is PLS-complete, however.) Let $P_2$ be some PLS-complete problem. Then define $P'$ which, on input $(x, i)$ for $i \in \{1, 2\}$, asks for a local optimum for input $x$ with respect to $P_i$. This problem inherits the PLS membership of $P_1, P_2$, inherits the PLS-completeness of $P_2$, and for the counting problem inherits the #P-completeness of $P_1$.
Similarly, one can construct an (artificial) PLS-complete problem for which it is NP-complete to decide if there is more than one local optimum. As in the previous argument, one "staples together" a PLS-complete problem $P_1$ as before, with a PLS problem $P_2$ which, on input a Boolean formula $\psi$, has more than one associated local optimum iff $\psi$ is satisfiable.
These sorts of constructions are somewhat unsatisfying because we're trying to build a search problem $Q$ that has two hardness properties, but the domain of $Q$ "splits" into two pieces, each of which may have only one of the two properties. Below I'll show how, given a search problem $P_1$ in PLS whose associated counting problem is #P-complete, and given a PLS-complete problem $P_2$, one can define a PLS problem $Q$ which is as hard as both counting for $P_1$ and search for $P_2$ in an "instance-by-instance" fashion.
Namely, we'll exhibit $Q$ such that solving the counting problem for $P_1$ on input $x$ efficiently reduces to solving the counting problem for $Q$ on input $x$, and the search problem for $P_2$ on input $x$ reduces to the search problem for $Q$ on input $x$.
For simplicity of presentation, we assume $P_1, P_2$ are such that on any input $x$ of length $n$, the candidate-solution space associated with $x$ is over bitstrings $y$ of length $n^c$ for some $c$ (but with different neighborhood structures for $P_1, P_2$). Let $F_i(x, y)$ be the fitness function associated with $P_i$.
On input $x \in \{0, 1\}^n$, the search space for $Q$ is over tuples $(y^1, y^2, z, b)$ where each $y^i$ is in $\{0, 1\}^{n^c}$, $z \in \{0, 1\}^{n^c + 1}$, and $b \in \{0, 1\}$. As the fitness function $F(x, (y^1, y^2, z, b))$ for $Q$, we define
$F(x, (y^1, y^2, z, b)) := F_1(x, y^1) + F_2(x, y^2)$ if $b = 1$,
$F(x, (y^1, y^2, z, b)) := ||y^1|| + ||z|| + F_2(x, y^2)$ if $b = 0$.
(That's Hamming weight above.)
For the neighborhood structure of $Q$, we connect each tuple $(x, (y^1, y^2, z, 1))$ (with $b = 1$) to all tuples $(x, ((y')^1, (y')^2, z', 1))$ such that
(A) $(x, y^i)$ is connected to $(x, (y')^i)$ according to $P_i$ for $i = 1, 2$, AND
(B) $z, z'$ differ in at most 1 coordinate.
For tuples with $b = 0$, we connect $(x, (y^1, y^2, z, 0))$ to all tuples $(x, ((y')^1, (y')^2, z', 0))$ such that
(A') $(x, y^2)$ is connected to $(x, (y')^2)$ according to $P_2$, AND
(B') $z, z'$ differ in at most 1 coordinate, as do $y^1, (y')^1$.
(Note, tuples with $b = 0$ are disconnected from those with $b = 1$.)
That's the definition of $Q$. The neighborhoods are of polynomial size as required, so $Q$ is in PLS.
Claim: The local optima for length-$n$ input $x$ according to $Q$ are exactly the following two disjoint sets:
(1) all tuples $(x, (y^1, y^2, z,1))$, where $(x, y^i)$ is a local optimum of $P_i$ for each of $i = 1, 2$ (and $z$ is arbitrary, and $b = 1$); and,
(2) all tuples $(x, 1^{n^c}, y^2, 1^n, 0))$, where $(x, y^2)$ is a local optimum of $P_2$, and where $z, y^1$ are both all-1s, and $b = 0$.
If you agree, then the PLS-hardness of $Q$ is immediate, since any local optimum $(x, (y^1, y^2, z,b))$ of $Q$ for input $x$ gives a local optimum $(x, y^2)$ of $P_2$ (for the same input $x$), and $P_2$ is PLS-hard.
Also, it follows from our Claim that the number $N(x)$ of local optima for $x$ under $Q$ equals $(2^{n^c + 1} N_1(x) + 1) \cdot N_2(x)$, where $N_i(x)$ is the number of local optima for $x$ under $P_i$. Now $N_2(x)$ is in the range $[1, 2^{n^c}]$, so we have
$N_2(x) = N_2(x)$ mod $2^{n^c + 1} = (2^{n^c + 1} N_1(x) + 1) \cdot N_2(x)$ mod $2^{n^c + 1} = N(x)$ mod $2^{n^c + 1}$.
So we can obtain $N_2(x)$ given $N(x)$. Then we can also obtain $N_1(x)$, by simple algebra:
$N_1(x) = \left(\frac{N(x)}{N_2(x)} - 1\right)/2^{n^c + 1}$. As $N_1(x)$ is #P-complete to compute, so is $N(x)$. Thus it's #P-complete to count local optima for $Q$ (and counting for $P_1$ reduces to counting for $Q$ on the same instance).
I don't know how to give such a reduction for combining PLS-hardness with NP-hardness of deciding uniqueness of local optima in an "instance-by-instance" fashion.
As for whether every PLS-complete search problem yields a #P-complete counting problem, I don't know this either. It seems related to the question of whether, for every NP-complete decision problem L and every polytime verifier $V(x, y)$ for $L$, the associated witness-counting problem is #P-complete. #P-completeness holds in all specific cases people have considered, and under some reasonably mild conditions, but is open in general. See this discussion.
For a specific, more natural problem $Q$ known to be PLS-complete, one might be able to establish #P-completeness for counting local optima by giving a PLS-reduction from say Matching to $Q$ that has some special properties appropriate for counting. Maybe the existing techniques are sufficient; I haven't tried to ascertain.
