Good question. There is in fact a really simple way of seeing there are Majorana modes on the edges, using the idea of `symmetry fractionalization'. In fact this principle works for all symmetry protected topological phases in 1D and it classifies all of them! The below explanation might seem a bit long, but I assure you it is all quite intuitive, and moreover you can then immediately apply the same logic to all SPTs you might know or encounter in the future :)
Generally, suppose your Hamiltonian has an on-site symmetry $U$ (meaning it acts on every site $U = U_1 \otimes U_2 \otimes \cdots$). If we have periodic boundary conditions, and the symmetry is not spontaneously broken in the ground state, then we expect that acting $U$ on the ground state should leave it invariant1. If we now cut the system to create (open) boundaries, then we expect our ground state to only change near the edges with a length scale given by the correlation length $\xi$. So since the ground state in the bulk is still the same, and we knew that ground state had to be invariant, we know that acting $U$ can only change our state near the left and right edge (with a scale given by $\xi$). So effectively the symmetry $U$ acts as $U = U_L U_R$ with these operators localized on the edges (it is said the symmetry `fractionalizes' onto the edges). Let's now see how this can classify and detect topological phases.
Note that the Kitaev chain only has terms which have an even number of fermionic operators. Physically this means that although the Hamiltonian doesn't preserve particle number (indeed its physical origin is related to superconductivity, where pairs of particles can be created), it does preserve whether you have an even or odd number of fermions. This is measured by the fermionic particle number parity $P = \prod (1-2n_i)$, which is a symmetry of the Hamiltonian. Rewriting this in terms of the Majorana operators, this is $\boxed{P = \prod i \gamma_{A,i} \gamma_{B,i}}$. Let's look at the case where it is clear there are Majorana modes: $\boxed{H = -\sum i \gamma_{A,i+1} \gamma_{B,i}}$. (Then $\gamma_{A,1}$ and $\gamma_{B,N}$ are free.) Note that these terms in the Hamiltonian commute and so for the ground state state subspace $i \gamma_{A,i+1} \gamma_{B,i} = 1$. Plugging this into $P$ one gets that in the ground state subspace $P = i\gamma_{A,1} \gamma_{B,N}$. Note that this only acts on the edge, consistent with the above general argument! In other words we can write $\boxed{P = P_L P_R}$ with $P_L = i\gamma_{A,1}$ and $P_R = \gamma_{B,N}$.
Now here comes the punchline: note that $P_L$ and $P_R$ are fermionic, so $P_L P_R = - P_R P_L$. Even if we now add extra couplings as in the original post, of course the expressions for $P_L$ and $P_R$ will change, however the anti-commutation cannot suddenly/discretely change to a commutation (and it has to be one of the two: since $P$ is bosonic, $P_L$ and $P_R$ are either both fermionic or bosonic). The only way it can change is if $P_L$ and $P_R$ start overlapping, which means $\xi$ has to diverge. This means it can only change if there is a phase transition! So as long as the system remains gapped, $P_L P_R = -P_R P_L$. But note that we then2 have two anticommuting symmetries. But by general linear algebra arguments3, the presence of two anticommuting symmetries implies a ground state degeneracy. This is exactly detecting the non-local fermionic mode formed by our two Majorana zero modes!
Footnotes:
1. The intuition is the following: suppose $U$ does not leave your ground state invariant and suppose our ground state is translationally invariant--consistent with using periodic boundary conditions--then by translation symmetry it changes the state everywhere, but that means $|\psi_\textrm{gs}\rangle$ and $U|\psi_\textrm{gs}\rangle$ are globally different, in which case we expect there to be spontaneous symmetry breaking.
2. Clearly $P$ is a symmetry, but you might wonder why $P=P_L P_R$ would imply that $P_L$ and $P_R$ are separately symmetries. Well the way to see this: decompose $H=H_L + H_R$ such that $H_R$ has no support near the left edge, same for $H_L$ and the right edge. Then $[P_L,H_R] = 0$. But since $[P,H_R] = 0$, that means that $[P_R,H_R] = 0$. Combined with $[P_R,H_L] = 0$, this shows that $P_R$ is a symmetry of $H$.
3. Suppose $|\psi_\textrm{gs}\rangle$ is a ground state. Then either $P_L|\psi_\textrm{gs}\rangle$ or $P_R |\psi_\textrm{gs}\rangle$ has to be a different state, otherwise $P_L P_R |\psi_\textrm{gs}\rangle = P_R P_L |\psi_\textrm{gs}\rangle$, yet we know they anticommute.
When you diagonalize $H_\text{BdG}$, you find the eigenvalues and the eigenvectors. Let us assemble these (column) vectors into a unitary matrix $U$ (this is always possible since $H_\text{BdG}$ is Hermitian), so that $U^\dagger H_\text{BdG} U=E$ is diagonal (i.e. the eigenvalues). Define $\Gamma_k = U^\dagger \psi_k$, then
$ H=\frac{1}{2} \Gamma_k^\dagger E \Gamma_k $
Therefore $\Gamma_k$ gives you the correct "single-particle" eigenstates of $H$ and can be used to build any excitations of the Hamiltonian. To be more explicit, the $i$-th row of $\Gamma_k$ is $\sum_j(U^\dagger)_{ij}\psi_{kj}$, with energy $E_i$.
Lastly, because the particle-hole symmetry of the Hamiltonian, the spectrum of $H_\text{BdG}$ always comes in $\pm E$ pairs, and correspondingly $\gamma_{E}$ and $\gamma_{-E}$ are actually related by hermitian conjugation: $\gamma_{-E}=\gamma_E^\dagger$.
Best Answer
The short answer is this: each line on your plot is denoting one fermionic mode (i.e. a fermionic creation/annihilation operator $c$), or equivalently each line corresponds to two Majorana modes (since a fermionic mode defines two Majorana operators $\gamma_1 = \frac{1}{2} \left(c+c^\dagger \right)$ and $\gamma_2 = \frac{1}{2i} \left(c-c^\dagger \right)$). Each line is represented twice though. It is not written in stone that you have to take the upper line, but just make sure you only count each set of two only once. Whether you take the upper or lower line comes down to choosing what you call the fermionic creation operator and the annihilation operator. I guess the clearest would be to just forget the bottom set of lines from the outset.
So take the point $\mu = 0$. On your plot we see two lines meeting there, but as you point out this is due to a redundancy in description, so in fact there is only one line there. In other words there is one fermionic mode $c$ with zero energy. This gives a twofold degeneracy: $|0\rangle$ and $c^\dagger |0\rangle$. Your plot shows that this mode does actually not stay at zero energy in the whole phase, but I guess this must be due to finite-size effects and as you take $L$ larger and larger, the line should approach zero until the phase transition point.