Here is the setting:
The likelihood of data is \begin{align} p(\boldsymbol{x} | \mu, \sigma^2) &= (\frac{1}{2\pi \sigma^2})^{\frac{n}{2}} exp\left\{ -\frac{1}{2\sigma^2} \sum\limits_{i=1}^n (x_i - \mu)^2 \right\} \nonumber \\ &= \frac{1}{(2\pi)^{n/2}} (\sigma^2)^{-n/2} exp\left\{ -\frac{1}{2\sigma^2} \left[ \sum\limits_{i=1}^n (x_i - \overline{x})^2 + n(\overline{x} - \mu)^2 \right] \right\}, \nonumber \end{align} and we use the Normal-inverse-Gamma as prior \begin{align} p(\mu , \sigma^2) &= \mathcal{N} (\mu | \mu_0 , \sigma^2 V_0) IG(\sigma^2 | \alpha_0 , b_0 ) \nonumber \\ &= \frac{1}{\sqrt{2\pi V_0}} \frac{b_0^{\alpha_0}}{\Gamma(\alpha_0)}\frac{1}{\sigma} (\sigma^2)^{-\alpha_0 - 1} exp\left( -\frac{1}{2\sigma^2} [V_0^{-1}(\mu - \mu_0 )^2 + 2b_0] \right). \nonumber \end{align} Then, the posterior can be derivated via \begin{align} p(\mu , \sigma^2 | \boldsymbol{x}) &\propto p(\boldsymbol{x} | \mu , \sigma^2 ) p(\mu , \sigma^2) \nonumber \\ &\propto \left[ (\sigma^2)^{-n/2} exp \left( -\frac{1}{2\sigma^2} \big[\sum\limits_{i=1}^b (x_i - \overline{x})^2 + n(\overline{x} - \mu )^2\big] \right) \right] \times \left[ \sigma^{-1} (\sigma^2)^{-\alpha_0 - 1} exp \left( -\frac{1}{2\sigma^2} \left[ V_0^{-1} (\mu - \mu_0 )^2 + 2 b_0 \right] \right) \right] \nonumber \\ &= \sigma^{-1} (\sigma^2)^{-(\alpha_0 + \frac{n}{2}) - 1} exp \left( -\frac{1}{2\sigma^2} \big[ V_0^{-1} (\mu - m_0 )^2 + 2 b_0 + \sum\limits_{i=1}^n (x_i - \overline{x})^2 + n(\overline{x} - \mu)^2 \big] \right) \nonumber \\ &= \sigma^{-1} (\sigma^2)^{-(\alpha_0 + \frac{n}{2}) - 1} exp \Big\{ -\frac{1}{2\sigma^2} \Big[ (V_0^{-1} + n)(\mu - \frac{V_0^{-1} m_0 + n\overline{x}}{V_0^{-1} + n})^2 + \big(b_0 + \frac{1}{2}\sum\limits_{i=1}^n (x_i - \overline{x})^2 + \frac{V_0^{-1} n}{2(V_0^{-1} + n)} (m_0 - \overline{x})^2 \big) \Big] \Big\} \nonumber \end{align} We recognize this is an unnormalized Normal-inverse-Gamma distribution, therefore \begin{align} p(\mu , \sigma^2 | \boldsymbol{x}) = NIG(\mu , \sigma^2 | m_n , V_n , \alpha_n , b_n ), \nonumber \end{align} where \begin{align} m_n &= \frac{V_0^{-1} m_0 + n \overline{x}}{V_0^{-1} + n} \nonumber \\ V_n^{-1} &= V_0^{-1} + n \nonumber \\ \alpha_n &= \alpha_0 + \frac{n}{2} \nonumber \\ b_n &= b_0 + \frac{1}{2}\sum\limits_{i=1}^n (x_i - \overline{x})^2 + \frac{V_0^{-1} n}{2(V_0^{-1} + n)}(m_0 - \overline{x})^2. \nonumber \end{align} As indicated in this paper (see Eq(200)), the last term can be further expressed as \begin{align} b_n &= b_0 + \frac{1}{2} \left[ m_0^2 V_0^{-1} + \sum\limits_{i=1}^n x_i^2 - m_n^2 V_n^{-1} \right]. \nonumber \end{align} But I fail to prove it, i.e., \begin{align} \sum\limits_{i=1}^n (x_i - \overline{x})^2 + \frac{V_0^{-1} n}{(V_0^{-1} + n)}(m_0 - \overline{x})^2 &= \left[ m_0^2 V_0^{-1} + \sum\limits_{i=1}^n x_i^2 - m_n^2 V_n^{-1} \right]. \nonumber \end{align}