Here is my attempt: it is incomplete, because it's missing the proof that the function defined in $(1)$ below indeed fails to be differentiable $\partial \mathcal{B}(p,q)$. Perhaps it's useful regardless of this gap.

The function \begin{equation} \tag{1} x \in M \mapsto \lvert d(x,p) - d(x,q) \rvert \end{equation} is Lipschitz, with zero set $\mathcal{B}(p,q)$. I think this is not differentiable on $\partial \mathcal{B}(p,q)$—except perhaps where one of $x \mapsto d(x,p)$ or $d(x,q)$ are not differentiable. By Rademacher's theorem, one would therefore have \begin{equation} \mathcal{H}^n(\partial \mathcal{B}(p,q)) = 0. \end{equation} Next we prove that $\mathcal{B}(p,q)$ has empty interior, whence $\mathcal{B}(p,q) = \partial \mathcal{B}(p,q)$ and the proof would be concluded.

Suppose otherwise, and let $U \subset \mathcal{B}(p,q)$ be an open set. As the cut locus of both $p$ and $q$ has Hausdorff dimension $n-1$, there is a point $z \in U$ that does not belong to either. There are two minimizing geodesics $\gamma_p$ and $\gamma_q$ connecting $z$ to $p$ and $q$ respectively. As $z$ is not in the cut locus of either $p$ or $q$, we can continue both geodesics a small amount beyond $z$, while retaining their minimizing properties. Denote these extensions still by $\gamma_p$ and $\gamma_q$, parametrized so that $\gamma_p(0) = p$, $\gamma_q(0) = q$.

Assume in what follows that $\gamma_p$ and $\gamma_q$ intersect transversely at $z$—otherwise the argument is similar, but a bit easier.

Let $D = d(z,p) = d(z,q)$. Then $\gamma_p(D + \epsilon) \in U$, and by assumption \begin{equation} D + \epsilon = d(\gamma_p(D + \epsilon),p) = d(\gamma_p(D + \epsilon),q). \end{equation} The curve obtained by concatenating $\gamma_q$ from time $0$ to $D$ with $\gamma_p$ from time $D$ to $D+\epsilon$ is a path from $q$ to $\gamma_p(D+\epsilon)$ with length $D+\epsilon$. As $d(\gamma_p(D+\epsilon),q) = D+\epsilon$, it would therefore be a minimizing geodesic. However, it is not even smooth!

Here is my attempt: it is incomplete, because it's missing the proof that the function defined in \eqref{1} below indeed fails to be differentiable on $\partial \mathcal{B}(p,q)$. Perhaps it's useful regardless of this gap.

The function \begin{equation} \tag{1}\label{1} x \in M \mapsto \lvert d(x,p) - d(x,q) \rvert \end{equation} is Lipschitz, with zero set $\mathcal{B}(p,q)$. I think this is not differentiable on $\partial \mathcal{B}(p,q)$—except perhaps where one of $x \mapsto d(x,p)$ or $d(x,q)$ are not differentiable. By Rademacher's theorem, one would therefore have \begin{equation} \mathcal{H}^n(\partial \mathcal{B}(p,q)) = 0. \end{equation} Next we prove that $\mathcal{B}(p,q)$ has empty interior, whence $\mathcal{B}(p,q) = \partial \mathcal{B}(p,q)$ and the proof would be concluded.

Suppose otherwise, and let $U \subset \mathcal{B}(p,q)$ be an open set. As the cut locus of both $p$ and $q$ has Hausdorff dimension $n-1$, there is a point $z \in U$ that does not belong to either. There are two minimizing geodesics $\gamma_p$ and $\gamma_q$ connecting $z$ to $p$ and $q$ respectively. As $z$ is not in the cut locus of either $p$ or $q$, we can continue both geodesics a small amount beyond $z$, while retaining their minimizing properties. Denote these extensions still by $\gamma_p$ and $\gamma_q$, parametrized so that $\gamma_p(0) = p$, $\gamma_q(0) = q$.

Assume in what follows that $\gamma_p$ and $\gamma_q$ intersect transversely at $z$—otherwise the argument is similar, but a bit easier.

Let $D = d(z,p) = d(z,q)$. Then $\gamma_p(D + \epsilon) \in U$, and by assumption \begin{equation} D + \epsilon = d(\gamma_p(D + \epsilon),p) = d(\gamma_p(D + \epsilon),q). \end{equation} The curve obtained by concatenating $\gamma_q$ from time $0$ to $D$ with $\gamma_p$ from time $D$ to $D+\epsilon$ is a path from $q$ to $\gamma_p(D+\epsilon)$ with length $D+\epsilon$. As $d(\gamma_p(D+\epsilon),q) = D+\epsilon$, it would therefore be a minimizing geodesic. However, it is not even smooth!