It's a subtle issue with the thermodynamic limit, i.e., taking the "quantization volume" to infinity and keeping the density constant. As long as the volume is finite you don't need to take the BEC particles (i.e., the macroscopic numer of particles occupying the ground state) separately, but when taking the thermodynamic limit you have to do so.
Take an ideal gas. Then the total number of particles at finite volume is the sum (sic!) over discrete momentum states (imposing periodic boundary conditions for the wave functions taking the volume to be a cube of length ##L##)
$$N=\sum_{\vec{p}} \frac{1}{\exp[\beta (E(\vec{p})-\mu)]-1}.$$
Obviously you must have ##\beta=1/(k_{\text{B}} T)>0## and ##\mu<E_0##. For non-relativistic particles ##E_0=0## and thus ##\mu<0##. There's no problem to accommodate any particle number though ##mu<0## at any temperature since the zero-mode contribution
$$N_0=\frac{1}{\exp(-\beta \mu)-1} \rightarrow \infty \quad \text{for} \quad \mu \rightarrow 0^{-}.$$
If you now go over to the thermodynamic limit, you naively make
$$\sum_{\vec{p}}=\frac{V}{(2 \pi \hbar)^3} \int_{\mathbb{R}^3} \mathrm{d}^3 p.$$
Since now ##\mathrm{d}^3 p=4 \pi p^2 \mathrm{d}p## the singularity at ##\vec{p}=0## for ##\mu \rightarrow 0^{-}## is gone, and there's a maximal finite number of particles for a given temperature. Thus if you have a fixed number of particles and lower the temperature, at some "critical temperature" you cannot accommodate all particles when taking this naive limit. The finite-volume calculation shows why: You have to keep a macroscopic number of particles occupying the zero mode.
So the correct limit of the phase-space distribution in the thermodynamic limit is not simply the Bose-distribution function but (with the usual convention that the phase-space measure is ##\mathrm{d}^3 x \mathrm{d}^3 p/(2 \pi \hbar)^3##)
$$f(\vec{x},\vec{p})=(2 \pi \hbar)^3 n_0 \delta^{(3)}(\vec{p}) + \frac{1}{\exp[\beta (E(\vec{p})-\mu)]-1}.$$
For a given density ##n=N/V## and energy density you have to adjust ##\beta## and ##\mu## such as to fulfill these constraints. If you cannot find a pair of ##\beta## and ##\mu## without taking into account the ##\delta## distribution, because ##n## is too large to do so, you must set ##\mu=0## and choose ##n_0>0## to accommodate the given density. In the thermodynamic limit you have a phase transition with the zero-mode occupation density ##n_0## as the order parameter, and that's what's called "Bose Einstein Condensation".