It is important not to confuse the union bound with further assumptions and simplifications made to our calculation of the probability of error.
The analysis starts with a constellation of only two points, $c_1$ and $c_2$, at a distance $d$:
The received point is $r = a_1 + n$; $n$ is two-dimensional Gaussian noise with variance $N_0$, and $a_1$ is equal to either $c_1$ or $c_2$. $C_1$ and $C_2$ are "decision regions"; if $r \in C_1$ then $r$ is closer to $c_1$ than $c_2$ and viceversa.
Assuming $a_1 = c_1$, we want to find the probability that $r$ is closer to $c_2$ than to $c_1$. We call this the pairwise error probability and write it as $$ \text{P}[c_1 \rightarrow c_2 ] = Q\left(\sqrt{\frac{d^2}{2N_0}}\right).$$
Now we'd like to estimate the error probability when there are more than two points involved, let's say $c_1$, $c_2$, $c_3$ and $c_4$. This is (in general) much more complicated, because it involves integrating the noise distribution over complicated regions. To avoid such calculations, we can use the union bound.
To calculate the union bound on the probability of symbol error in the constellation you posted, pick one point (say, $s_0$) and calculate its pairwise error probability to each of the other three points. Then, add them; this is an overestimate of the actual error probability $P_e$. This is the union bound on the probability of error:
\begin{align}
P_e &\leq \text{P}[s_0 \rightarrow s_1 ] + \text{P}[s_0 \rightarrow s_2 ] + \text{P}[s_0 \rightarrow s_3 ] \\
&= Q\left(\sqrt{\frac{d_{01}^2}{2N_0}}\right) + Q\left(\sqrt{\frac{d_{02}^2}{2N_0}}\right) + Q\left(\sqrt{\frac{d_{03}^2}{2N_0}}\right).
\end{align}
Note: in principle, we should calculate the probability of error $P(e|s_i)$ for each constellation point $s_i$ and then average. In this case, we can avoid this calculation because the constellation is symmetrical.
Now we may make a further assumption: if, say, $d_{02}$, the distance between $s_0$ and $s_2$, is much larger than $d_{01}$ and $d_{03}$, and given that $Q$ decreases exponentially fast, we may say that $\text{P}[s_0 \rightarrow s_2 ]$ is much smaller than $\text{P}[s_0 \rightarrow s_1 ]$ and $\text{P}[s_0 \rightarrow s_3 ]$ and we may remove it from the equation:
\begin{align}
P_e &\approx \text{P}[s_0 \rightarrow s_1 ] + \text{P}[s_0 \rightarrow s_3 ] \\
&= Q\left(\sqrt{\frac{d_{01}^2}{2N_0}}\right) + Q\left(\sqrt{\frac{d_{03}^2}{2N_0}}\right).
\end{align}
Note that the equation now says $\approx$ instead of $\leq$ -- since we removed a positive quantity from the right-hand side, mathematically we can no longer say that we have an upper bound. However, we know that, for high SNR, the approximation will be good. This is no longer the union bound.
To summarize, we may say that points at minimum distance will dominate the error calculation. If not all points have the same number of neighbors at minimum distance, then an average should be taken (as explained in my answer for 16-QAM).
