I believe the issue is that the following limits:
- the difference quotient needed to compute derivatives of the metric that feed into the curvature tensor, $\lim_{h\rightarrow 0}\frac{g_{\mu\nu}(h)-g_{\mu\nu}(0)}{h}$
- $b_0 \rightarrow 0$
do not commute. As pointed out by @Anders Sandberg in the comments, the reason this issue arises is because the flat space limit of the Morris-Thorne wormhole does not just involve a change in geometry (ie the metric), but of the topology.
You find an issue when taking the derivative first and then sending $b_0\rightarrow 0$, while I claim the right limit to recover flat space is to send $b_0\rightarrow 0$ and then take the derivative. (Actually if you literally do that then the flat space limit is trivial, but below I'm going to sketch a different way you can do it to make the ordering ambiguity hopefully more transparent).
To spell it out a bit more, the result that the curvature is infinite at $\ell=0$ actually does make sense if you consider that $b_0$ is the width of the wormhole, and that the limit $b_0\rightarrow 0$ is telling us to consider a wormhole of an infinitesimally small width. Since the wormhole is infinitely narrow, the curvature associated with the wormhole is confined to an infinitesimally small volume, and diverges.
You can see what happens more explicitly by computing some of the Christoffel symbols. For example,
\begin{equation}
\Gamma^\theta_{\ell\theta} = \frac{1}{2}g^{\theta\theta} \partial_\ell g_{\theta\theta} = \frac{\ell}{\ell^2+b_0^2}
\end{equation}
Now let's first differentiate this assuming $-\infty < \ell < \infty$, which is appropriate for the wormhole geometry. We note that $\Gamma^\theta_{\ell\theta}$ is discontinuous at $\ell=0$, so we expect $\ell$ derivatives to pick up a delta-function-type singularity at $\ell=0$. More explicitly, the Riemann curvature tensor will involve derivatives of the Christoffel symbols, such as
\begin{equation}
\partial_\ell \Gamma^\theta_{\ell\theta} = -\frac{2\ell^2}{\left(\ell^2+b_0^2\right)^2} +\frac{1}{\ell^2+b_0^2} = \frac{-\ell^2 b_0^2}{\left(\ell^2 + b_0^2\right)^2} = \frac{\delta_{b_0}(\ell)}{b_0} + O(b_0^0)
\end{equation}
where $\delta_{b_0}(\ell)/b_0$ is the representation of the delta function you discussed in your question
\begin{equation}
\frac{\delta_{b_0}(\ell)}{b_0} \equiv \frac{b_0^2}{\left(\ell^2 + b_0^2\right)^2}
\end{equation}
and $O(b_0^0)$ represents terms that are are finite, or vanish, in the limit $b_0 \rightarrow 0$. This is essentially how the delta function shows up in your result.
However, for flat space we do not want to consider the full range $-\infty \leq \ell \leq \infty$. The flat space limit of the spacetime should only consider one "sheet" of the spacetime, and ignore the throat. We can take this into account by considering only the part of the spacetime for $\ell \gg b_0$, and then take the limit $b_0 \rightarrow 0$, and then allow $\ell$ to go to zero (since if $b_0=0$ then $\ell$ can be made arbitrarily small while still being in the $\ell\gg b_0$ regime). If we do this, then
\begin{equation}
\frac{b_0^2}{\ell^2 + b_0^2} \approx \frac{b_0^2}{\ell^2} \rightarrow 0
\end{equation}
in the limit $b_0\rightarrow 0$, as opposed to a delta function.
This is obviously just a sketch and not a full argument or proof. If you wanted to formalize this more, I think you could say you will require $\ell > k + b_0$ for some arbitrary constant $k$, and you can make $b_0^2/(\ell^2+b_0^2)^2$ arbitrarily small by choosing a large enough $k$. Then you take the limit $b_0 \rightarrow 0$ and then the limit $k\rightarrow 0$. I tend to be fairly physics-y and not incredibly interested in mathematical rigor beyond about this point (if I was going to rely on this limit as part of some research I might do some numerical examples to check that these limits can really be done consistently in the way I am outlining or try to find a collaborator to run the argument by and ask them to poke holes in it after writing out the details more careefully), but I am sure that there are others here that can give a more formal argument than the one I'm sketching.
The moral is that you can recover flat spacetime in the $b_0\rightarrow 0$ limit, but you need to be careful about the order of limits so you really do recover flat spacetime, and not a Morris-Thorne spacetime with an infinitely narrow wormhole.
