The error probability of 16QAM
In general, the symbol error probability is given by

$$P_E=\sum_{i=1}^MP[s_i]P[E|s_i]\tag{1}$$

where $P[s_i]$ is the probability that symbol $s_i$ was sent, $P[E|s_i]$ is the error probability given that $s_i$ was sent, and $M$ is the number of symbols. If all symbols are equally likely, $(1)$ simplifies to

$$P_E=\frac{1}{M}\sum_{i=1}^MP[E|s_i]\tag{2}$$

For $16$-QAM it is possible to derive an exact expression for the error probabilities $P[E|s_i]$ (and hence for the symbol error probability $P_E$). You can find the derivation in most digital communication text books. A good approximation (better than the union bound) is given in this answer.

The union bound approximation gives an upper bound on $P[E|s_i]$:

$$P[E|s_i]\le\sum_{j\neq i}P[s_j\text{ chosen}|s_i\text{ sent}]\tag{3}$$

Under the given assumptions we have

$$P[s_j\text{ chosen}|s_i\text{ sent}]=Q\left(\sqrt{\frac{d^2_{ij}}{2N_0}}\right)\tag{4}$$

where $d_{ij}$ is the distance between symbols $s_i$ and $s_j$. With $(3)$ and $(4)$ the union bound approximation of $(2)$ becomes

$$P_E\le\frac{1}{M}\sum_{i=1}^M\sum_{j\neq i}Q\left(\sqrt{\frac{d^2_{ij}}{2N_0}}\right)\tag{5}$$

Note that generally you will get a sum of $Q$ functions with different arguments because the distances $d_{ij}$ are usually not all the same. This is also the case for $16$-QAM.

The solution you provided is still more conservative (hence, worse) than the union bound. That solution uses the fact that the $Q$ function is monotonically decreasing, which implies

$$Q\left(\sqrt{\frac{d^2_{ij}}{2N_0}}\right)\le Q\left(\sqrt{\frac{d^2_{min}}{2N_0}}\right)\tag{6}$$

with $d_{min}=\min_{i\neq j}\{d_{ij}\}$ being the smallest occurring distance between symbols. Plugging $(6)$ into $(5)$ we finally get the bound

$$P_E\le (M-1)Q\left(\sqrt{\frac{d^2_{min}}{2N_0}}\right)\tag{7}$$

which, as mentioned above, is even less tight than the union bound.