On the last question, I am not sure how good you are at the representation theory, but the following fact is true: take so(d,2) (we need so(3,2) for this work), use the conformal base, i.e. Lorentz generators $L_{ab}$, translations $P_a$, conformal boosts $K_a$ and dilatation $D$, $a,b=1..d$. $P$ and $K$ behave as raising/lowering generators with respect to $D$, $[D,P]=+P$, $[D,K]=-K$. Take the vacuum to carry a spin-s representation of the Lorentz algebra and a weight $\Delta$ with respect to $D$, i.e. $|\Delta\rangle^{a_1...a_s}$. When $\Delta=d+s-2$, there is a singular vector, $P_m|\Delta\rangle^{ma_2...a_s}$. This is a standard representation theory: finding raising/lowering operators, defining vacuum, looking for singular vectors. Actually, singular vectors are exactly the conformally-invariant equations one can impose.
On the field language this means that $\partial_m J^{m a_2...a_s}=0$ is a conformally invariant equation iff the conformal dimension of $J$ is $\Delta=d+s-2$. Despite the fact that $J^{a_1...a_s}$ is a good conformal operator for any value of the conformal dimension, only for $d+s-2$ its divergence decouples. (Perhaps you have seen $L_{-2}+\alpha L_{-1}^2$ as a singular vector in the Virasoro algebra, now it is replaced with $P_m$ or $\partial_m$).
Now, having $J^{a_1..a_s}$ of weight $\Delta$ we can consider its contragradient representation or on the field language couple it via $\int \phi_{a_1..a_s}J^{a_1...a_s}$ to some other field $\phi$. That we need a conformally invariant coupling implies $\Delta_\phi=d-\Delta_J=s-2$. Not surprisingly something special must happen for $\Delta_J=d+s-2$.
$$\int (\phi_{a_1...a_s}+\partial_{a_1}\xi_{a_2...a_s})J^{a_1...a_s}=\int \phi J-\int \phi_{a_1...a_s}\partial_m J^{ma_2..a_s}=\int \phi J$$
we see that a statement that is dual to the conservation of $J$ is the gauge invariance of $\phi$.
I have not read the paper yet, but as far as I can see they play with the dimension of $J$ and for $d+s-2$ and $2-s$ it describes a conserved tensor and a gauge field just because of representation theory of the conformal group (decoupling of certain null states). At any given moment of time in the paper $J$ has some fixed dimension and is either a conserved tensor, a gauge field or just a spin-s conformal field of generic dimension $\Delta$.
On the last but one, you are right in that gauge invariance has a little to do with conformallity. The answer is spin and dimension dependent. For $s=0$ there is $m^2$ for which the scalar is conformal. For $s=1$ and certain $m^2$ the Maxwell field is a gauge field but the Maxwell equation is conformal in $d=4$ only. Beyond $d=4$ a gauge spin-one field is not conformal, or a spin-one conformal field is not a gauge field. For $s\geq2$ the situation is even more tricky: in $AdS_4$ the gauge fields are conformal, but in Minkowski space they are not conformal (in terms of gauge potentials $\phi_{\mu_1...\mu_s}$). You may have a look at http://arxiv.org/abs/0707.1085
On the second, first of all the transversality is on the right place in 5.1. Secondly, your confusion (inspired by my answer to another question) is that there are two different classes of fields people are interested in. First is the class of usual particles, where we talk about representations of the Poincare algebra $iso(d-1,1)$ if we are in $d$-dimensional Minkowski space or $so(d-1,2)$ and $so(d,1)$ if we are in anti de Sitter ot de Sitter (there we need harmonicity, tracelessness, transversality). Conformal fields are in the second class. Conformal means that it must be a representation of the conformal group $so(d,2)$ for Minkowski-$d$, note that $iso(d-1,1)\in so(d,2)$. The conformal group of anti de Sitter-$d$ is also $so(d,2)$. Note that the symmetry algebra of AdS-$(d+1)$ is exactly the conformal group of Minkowski-$d$. So when we talk about conformal fields we are interested in reprsetations of $so(d,2)$ (the signature can vary depending on the problem, it is some real form of $so(d+2)$). I would like to stress that conformal fields in d-dimensions are in one-to-one correspondence with usual fields in $AdS_{d+1}$, for the algebra is the same, which is at the core of AdS/CFT correspondence.
For example, a spin-$0$ field in Minkowski space obeys $\square \phi=0$. It gives rise to an irreducible representaion of $iso(d-1,1)$. Coincidentally, the same representation turns out to be an irreducible representation of a bigger algebra, $so(d,2)$, the conformal algebra. It is a coincidence. There exists also a spin-$0$ conformal field of weight $\Delta$, say $\phi_\Delta(x)$. Without imposing any equations it is an irreducible representation of $so(d,2)$. As a representation of its subalgebra $iso(d-1,1)$ it decomposes into an intergral of representations (Fourier) and is highly reducible. There is a special weight $\Delta=(d-2)/2$ for which $\phi_\Delta(x)$ is reducible and the decoupling of null states is achived via $\square \phi=0$ (analogous to the conservation of $J$ above). Note that $J$ above is an irreducible representation of $so(d,2)$ but it is highly reducible under $iso(d-1,1)$. For special weight $d+s-2$ we have to impose the conservation condition in order to project out the null states, but again the conserved tensor is an irreducible of $so(d,2)$ and reducible under $iso(d-1,1)$. So your confusion is because the fields are conformal, these are representations of a bigger algebra, they are more 'fat' and require less equations (even no at all) to project onto an irreducible.
$S^3$ is the analog of Minkowski-$3$ (compactified and Euclidian), then $so(4)$ is the analog of $iso(3,1)$ and they are interested in normalizable functions, these are the spherical harmonics or polynomials depending on coordinates. Then they discuss labelling of these representations using $so(4)\sim su(2)\oplus su(2)$ and proceed to doing some integrals.
