The splitting functions tell you the probability for a given radiance to occur, and some exhibit manifest divergences (e.g. as $z \rightarrow 0$) signifying a growth in that particular branching in a certain part of the phase space. The correct statements are $$\int_0^1 z \, \text{d}z \, (P_{q \leftarrow q}(z) + P_{g \leftarrow q}(z)) = 0\,\,\,\,\,\,\,;\,\,\,\,\,\,\,\,\int_0^1 z\, \text{d}z\, (2n_f P_{q \leftarrow g}(z) + P_{g \leftarrow g}(z)) = 0 $$
They may be derived by totalling the sum of all fractional constituent nucleon momentum to unity, means $$ \int_0^1 \text{d}x \,x \, \left[ \sum_{i} \left( q_i(x) + \bar{q}_i(x) \right) + g(x) \right] = 1 $$ and taking the derivative w.r.t $\mu^2 \,d/d \mu^2$ in each of the terms as follows $$\sum_i \mu^2 \frac{d}{d \mu^2} ( q_i(x) + \bar{q}_i(x) ) = \frac{\alpha_s}{2 \pi} \sum_i \left( q_i \otimes P_{q \leftarrow q} + 2 g \otimes P_{q \leftarrow g} \right),$$ $$ \mu^2 \frac{d}{d \mu^2} g(x) = \frac{\alpha_s}{2 \pi} \sum_j \left( q_j \otimes P_{g \leftarrow q} + g \otimes P_{g \leftarrow g} \right),$$ where $P_{\bar q \leftarrow q} = 0$ at leading order and $P_{\bar q \leftarrow g} = P_{q \leftarrow g}$ with $$(X_i \otimes P_{i \leftarrow j})(x) \equiv \int_{x}^1 \frac{d y}{y} X_i(y) P_{i \leftarrow j} \left(\frac{x}{y} \right).$$ Their addition gives $$\sum_i \mu^2 \frac{d}{ d \mu^2} \left( q_i(x) + \bar{q}_i(x) + g(x) \right) = \frac{\alpha_s}{2 \pi} \left( \sum_i q_i \otimes (P_{q \leftarrow q} + P_{g \leftarrow q}) + g \otimes ( P_{g \leftarrow g} + 2 n_f P_{q \leftarrow g} ) \right)$$ and so $$ \frac{\alpha_s}{2 \pi} \int_0^1 \text{d} x \, x \, \left( \sum_i q_i \otimes (P_{q \leftarrow q} + P_{g \leftarrow q}) + g \otimes ( P_{g \leftarrow g} + 2 n_f P_{q \leftarrow g}) \right) = 0.$$ Now, noting that $$ \int_0^1 \text{d}x \int_x^1 \text{d}y = \int_0^1 \text{d}y \int_0^y \text{d}x,$$ and making a change of variables $x = yz$ decouples the $x$ and $y$ integrations allowing for the relations given above for all $q_i$ and $g$. These relations are valid only at leading order, at the next order additional flavour effects come into play which makes the DGLAP evolution matrix structurally more involved.
Additional comments:
The factor of $z$ included means it is the total momentum fraction that is conserved, not the number density. In this sense, it is $zP_{ij}$ that are the probability densities over $z$. Note the sum of probabilities of all possible splittings plus probability for no splitting is what sums to unity. You can trade the non integrable singularity at $z=1$ in the splitting functions for something integrable plus explicit IR divergences centred at $z=1$, this is done practically via the plus prescription. Schematically, this gives the following
$$\int_0^1 \text{d}z\, z \,( P^0_{g \leftarrow g} + P_{q \leftarrow g} + A_0 \delta (1-z) + \delta(1-z) ) = 1 \Rightarrow \int_0^1 \text{d}z\, z \,(P_{g \leftarrow g} + P_{q \leftarrow g} ) = 0, $$ where $P_{g \leftarrow g}$ is the plus-prescribed (regularised) splitting function and $A_0$ is of course tuned to cancel the divergence at $z=1$.
