confinementquarksstandard-model
Quarks combine to form composite particles called hadrons, the most stable of which are protons and neutrons, the components of atomic nuclei.
Due to a phenomenon known as color confinement, quarks are never directly observed or found in isolation; they can be found only within baryons or mesons.
This sentence makes me very nervous: Due to a phenomenon known as color confinement
This sentence is Like I want to prove something spurious to save the subject (quark).
The cross-section for $$ e^+ + e^- \to q + \bar{q} $$ goes by the square of the quark charge (times the number of colors). Now, the quarks can not be observed in isolation because they hadronize.
However the cross-section for $$ e^+ + e^- \to \mu^+ + \mu^- $$ is identical except for going by the muon charge squared.
So, a measurement of $$ R = \frac{\sigma_{e^+ + e^- \to \text{hadrons}}}{\sigma_{e^+ + e^- \to \mu + \mu^-}} $$ is a measurement of $$ \frac{\text{# of colors} \times \sum_\text{accessible flavors} q^2_\text{flavor}}{q^2_\mu} = \frac{\text{# of colors} \times \sum_\text{accessible flavors} q^2_\text{flavor}}{1} \quad .$$
The accessible flavors depend on the center of mass energy, so it is possible to observe the increases as the energy rises past successive quark masses (times 2).
This figure :
shows $R$ over a range of center-of-mass energy (Mandelstam variable $s$) that covers the range from only including the "light" quarks (up, down and strange) through including all the quarks through the bottom with enough range to show the long plateau about the bottom quark threshold.
The results are consistent with three colors and the usual charge assignments (up-like quarks are +2/3 and down-like quarks are -1/3) from the baryon spectrum.
Each quark $q$ and antiquark $\bar{q}$ transform in the fundamental representation ${\bf 3}$ and antifundamental representation $\bar{\bf 3}$ of the color group $SU(3)_C$, respectively.
However, a hadron has to be a color singlet ${\bf 1}$, due to color confinement.
In particular the center $$\mathbb{Z}_3~:=~\mathbb{Z}/3\mathbb{Z}~=~\{1,e^{\pm 2\pi i/3}\} ~\subseteq ~SU(3)_C$$ of the color group $SU(3)_C$ should be in a trivial representation of the cyclic group $\mathbb{Z}_3$. This precisely happens if the quark number is a multiple of 3. This answers OP's title question.
Examples:
A single quark $q$ transforms in the fundamental representation ${\bf 3}$ of $SU(3)_C$, and is hence not allowed. See also related Phys.SE post here.
A diquark $qq$ belongs to the tensor representation ${\bf 3}^{\otimes 2}:={\bf 3}\otimes{\bf 3}\cong\bar{\bf 3}\oplus{\bf 6}_S$, which we have decomposed in irreps of $SU(3)_C$. This contains no singlet ${\bf 1}$, and is hence not allowed. See also related Phys.SE posts here and here about $SU(3)$ tensor representations.
In a meson, the quark-antiquark-pair $q\bar{q}$ belongs to ${\bf 3}\otimes\bar{\bf 3}\cong{\bf 1}\oplus{\bf 8}_M$, which contains a singlet ${\bf 1}$, and is hence allowed.
The third tensor product is ${\bf 3}^{\otimes 3}={\bf 3}\otimes{\bf 3}\otimes{\bf 3}\cong {\bf 1}\oplus 2\cdot{\bf 8}_M\oplus{\bf 10}_S$. In a baryon, the three quarks $qqq$ form a totally antisymmetric representation $\wedge^3 {\bf 3}\cong {\bf 1}$ of $SU(3)_C$, which is isomorphic to a singlet ${\bf 1}$, and is hence allowed. See also related Phys.SE post here. The lightest baryon, the proton, is stable in the standard model due to baryon/quark number conservation. (However, see the hypothetical proton decay.)
The tensor product ${\bf 3}\otimes{\bf 3}\otimes\bar{\bf 3}\cong 2\cdot{\bf 3}\oplus {\bf 6}_M\oplus {\bf 15}_M$ contains no singlet ${\bf 1}$, so the combination $qq\bar{q}$ is not allowed.
The fourth tensor product is ${\bf 3}^{\otimes 4}\cong 3\cdot{\bf 3}\oplus 2\cdot{\bf 6}_M\oplus 3\cdot{\bf 15}_M\oplus{\bf 15}_S$. This contains no singlet ${\bf 1}$, so four quarks $qqqq$ are not allowed.
A "molecule" of mesons and baryons, such e.g. a tetraquark $q\bar{q}q\bar{q}$ or a pentaquark $qqqq\bar{q}$, is also allowed, but is obviously heavier. See also related Phys.SE posts here, here, and here.
TL;DR: The number of quarks minus the number of antiquarks should be divisible by 3.
References:
Best Answer
It is believed for good reason that confinement isn't just some made up thing. The good reason is from Lattice QCD. While it is not analytic proof, simulations of QCD show that QCD is indeed confining naturally. That is, confinement is already hidden inside QCD and happens naturally. Furthermore, for high enough temperatures QCD is deconfining, and the result is a quark-gluon plasma, where the two run around like a plasma. It is unfortunate that as of yet little experimental evidence exists for a deconfining phase of QCD, at least to my knowledge. But computational work on QCD sheds a great deal of insight into the genuine physics of QCD. As a note, deep inelastic scattering provides evidence for quarks, as the experiment resembles Rutherford scattering, which phenomenologically helped develop an accurate model of the atom eventually. Funny enough the "planet model" of the atom is an inaccurate interpretations in certain circumstances.
That being said, our interpretations of physics and physical systems changes (period). So almost inevitably, while quarks will probably stick around, how we view them will change.
If you are interested in constraints put on QCD by computational simulations look into Lattice QCD and Gauge Fields on a Lattice. It has been an important step in tapping into non-perturbative QCD.