You're right that $q_{ij}$ is the probability of seeing the sample point $z_i$ in the cluster with centroid $u_j$ under the assumption that points are distributed by the Student t-distribution.
The target distribution, $p_{ij}$ represents ideal clusters. In the most extreme case, clusters could be defined by a delta distribution over a ball of a given radius around the cluster centroid (i.e. if the point is within the radius, it has probability 1 of being a member of the cluster). However, such a discontinuous definition is computationally impractical and makes the optimization problem difficult. Instead, the authors opted for a softer constraint that still puts most of the probability mass around the cluster centroid. They achieve this by looking at the squared Student t-distribution and normalizing by the relative frequency of the cluster, (a.k.a. the soft cluster frequency given by $\sum_{i} q_{ij}$).
In this context, the KL divergence, $\text{KL}(P||Q)$, is interpreted as the amount of information lost when you use the distribution Q to approximate the distribution P. Specifically, it is the error (in bits) of the observed distribution over the data points Q when you wanted the ideal cluster distribution P.
Finally, the authors consider points that might exist in a high dimensional space (the authors reference images), and seek to find a low dimensional representation of them (the features) which preserves local structure of the points but also identifies a grouping of the points into clusters. The authors capture this trade-off with a joint error function summing the reconstruction loss of the auto-encoder with the above clustering error.