In Thermodynamics and in Statistical Mechanics, many routes may reach the same result. However, mixing them, even if the final result is a correct formula, may leave some doubts about the consistency of the arguments with the starting hypothesis and definition.
In this specific case, I see a potential conceptual problem in mixing the first statement of the question:
Consider the grand canonical Gibbs distribution...
and the formula for the chemical potential
$$ \mu = -T \frac{\partial{S}}{\partial{N}}. $$
Such a formula is correct from the thermodynamic point of view, but, to be meaningful in the present context, one should specify which state variables are kept constant in the partial derivative, i.e., which variables $S$ is a function of. Here, some ambiguity of the description enters.
Indeed, it is true that one traditional way of deriving and interpreting the grand canonical ensemble is through the analysis of a "big" system with a finite number of particles, partitioning it into a "system of interest" with a fixed subvolume $V$ of the total volume, considering all the possible constraints corresponding to fix the values of the number of particles of the subvolume, and looking for the equilibrium condition with respect to the exchange of particles at a fixed temperature.
Following this way of deriving the grand canonical ensemble, or even through a direct thermodynamic analysis, it turns out that the above formula for chemical potential should involve a partial derivative of $S$ with respect to $N$ at fixed $T,V$ (or with respect to $E,V$; $E$ being the internal energy of the system). Thus, one should know $S(T,V,N)$ (or $S(E,V,N)$) for the system of interest.
However, once the grand canonical formulae have been obtained, the grand canonical partition function $\Xi $ and the resulting Grand Potential $\Omega = -k_BT \log \Xi$ are functions of $T,V,\mu$. The number of particles is an ill-defined concept in the grand-canonical ensemble. It cannot be fixed since it fluctuates. The average number of particles $<N>$ should replace it, which can be obtained as $-1/k_BT$ times the partial derivative of $\Omega$ with respect to $\mu$ at fixed $T,V$.
Said in another way, a formula like $$ \mu = -T \left. \frac{\partial{S}}{\partial{N}}\right|_{T,V} = -T \left. \frac{\partial{S}}{\partial{N}}\right|_{E,V} $$ looks like a foreigner in the grand canonical ensemble because $N$ is not a well defined quantity (neither $E$, referring to the second form). On the other side, $\mu$, exactly like $T$ and $V$ is an independent physical parameter characterizing the ensemble and the thermodynamic macrostate one is interested in.
Put this way, the question looks like a similar question which could be asked for the canonical ensemble: since we know that $\frac{1}{T} = \frac{\partial{S}}{\partial{E}}$, if the system may take different energies, and then different entropies, why the temperature is unique?
Once one has already gone through the derivation of the canonical ensemble, the question, in this form, is not well-posed because energy is not an independent variable anymore. Missing this point may drive towards even worse ill-posed questions (sometimes appeared even in the scientific literature) like wondering about temperature fluctuations in the canonical ensemble or fluctuations of chemical potential in the grand-canonical ensemble. By construction, both quantities are fixed, and speaking about their fluctuations is meaningless.