Suppose we had a lattice Hamiltonian $H$ which was symmetric under the particle-hole transformation
$$ c_n \mapsto U^\dagger c_nU=(-1)^nc^\dagger _n$$
such that $[H,U] = 0$, where $c_n$ are Fermionic operators obeying $\{ c_n , c_m^\dagger \} = \delta_{nm}$ and $\{ c_n,c_m \} = \{ c^\dagger_n, c^\dagger_m \} = 0$ where the indices label the lattice sites of the system.
If the ground state $|\Omega \rangle$ of $H$ is unique, then it would be an eigenstate of the particle-hole transformation as $U |\Omega\rangle = e^{i \theta} |\Omega\rangle$, where the eigenvalue is some phase as $U$ is unitary. From this, if I calculated the nearest-neighbour correlation, I would find
\begin{align*} \langle \Omega | c_n^\dagger c_{n+1} |\Omega \rangle & = \langle \Omega | U^\dagger c_n^\dagger c_{n+1} U |\Omega \rangle \\ & = \langle \Omega |( U^\dagger c_n^\dagger U)(U^\dagger c_{n+1} U) |\Omega \rangle \\ &= (-1)^{2n+1}\langle \Omega| c_n c_{n+1}^\dagger |\Omega \rangle \\ &=- \langle \Omega |(-c_{n+1}^\dagger c_n) |\Omega \rangle \\ &= \langle \Omega| c^\dagger_{n+1} c_n |\Omega \rangle \\ &= \langle \Omega | c_n^\dagger c_{n+1} |\Omega \rangle^* \end{align*}
therefore I find it is real. A similar calculation would also show that this implies half-filling $\langle \Omega| c^\dagger_n c_n |\Omega \rangle = \frac{1}{2}$ as shown in this answer.
This calculation relied upon the fact that $|\Omega\rangle $ is unique, however many systems will have ground state degeneracy in which case $U |\Omega \rangle \neq e^{i\theta} |\Omega\rangle$ in general. One can always find a state in the ground state subspace that is an eigenstate however, and the above results will apply.
(Edit) Example Hamiltonian
An example of an interacting Hamiltonian that has particle-hole symmetry is given by
$$ H = -t \sum_n c_n^\dagger c_{n+1} + \lambda \sum_n c^\dagger_n c_{n+1} c_{n+2}^\dagger c_{n+3} + \mathrm{h.c.} = H_0 + H_\mathrm{int} $$
where $t,\lambda \in \mathbb{R}$. We have \begin{align*} U^\dagger H_0 U & = -t \sum_n (-1)^{2n+1} c_n c^\dagger_{n+1} + \mathrm{h.c.} \\ & = t \sum_n c_n c^\dagger_{n+1} + \mathrm{h.c.} \\ & = -t\sum_n c_{n+1}^\dagger c_n + \mathrm{h.c} \\ & = H_0 \end{align*}
And for the interaction part \begin{align*} U^\dagger H_\mathrm{int} U & = \lambda \sum_n (-1)^{4n + 6} c_n c^\dagger_{n+1} c_{n+2} c^\dagger_{n+3} + \mathrm{h.c.}\\ & = \lambda \sum_n c_{n+3} c^\dagger_{n+2} c_{n+1} c_n^\dagger + \mathrm{h.c.} \\ & = H_\mathrm{int} \end{align*}
where I have interchanged parts with the hermitian conjugate and anti-commuted fermions past each other.
As this Hamiltonian is interacting, I am not sure how I would derive the correlation matrix of this, but all I would like to know is if this system has half-filling and real nearest-neighbour correlations in the ground state.
My question
In the case of no degneracy we can say that particle-hole symmetry implies half filling $\langle c^\dagger_n c_n \rangle = \frac{1}{2}$ and the correlation $\langle c_n^\dagger c_{n+1} \rangle \in \mathbb{R}$. However, in the case of degeneracy the above derivation would not apply unless we chose a state in the ground state subspace that is an eigenstate of $U$. My question is, am I allowed to choose the particle-hole symmetric ground state over another as "special" and can I always assume that particle-hole symmetry implies half-filling and real nearest-neighbour correlations?