Abstract: Recent results establish for the hard-core model (and more generally for 2-spin antiferromagnetic systems) that the computational complexity of approximating the partition function on graphs of maximum degree $\Delta$ undergoes a phase transition that coincides with the uniqueness/non-uniqueness phase transition on the infinite $\Delta$-regular tree. For the ferromagnetic Potts model we investigate whether analogous hardness results hold. Goldberg and Jerrum showed that approximating the partition function of the ferromagnetic Potts model is at least as hard as approximating the number of independent sets in bipartite graphs, so-called #BIS-hardness. We improve this hardness result by establishing it for bipartite graphs of maximum degree $\Delta$. To this end, we first present a detailed picture for the phase diagram for the infinite $\Delta$-regular tree, giving a refined picture of its first-order phase transition and establishing the critical temperature for the coexistence of the disordered and ordered phases. We then prove for all temperatures below this critical temperature (corresponding to the region where the ordered phase “dominates'') that it is #BIS-hard to approximate the partition function on bipartite graphs of maximum degree $\Delta$. As a simple corollary of this result, we obtain that it is #BIS-hard to approximate the number of $k$-colorings on bipartite graphs of maximum degree $\Delta$ whenever $k\leq \Delta/(2\ln \Delta)$. The #BIS-hardness result for the ferromagnetic Potts model uses random bipartite regular graphs as a gadget in the reduction. The analysis of these random graphs relies on recent results establishing connections between the maxima of the expectation of their partition function, attractive fixpoints of the associated tree recursions, and induced matrix norms. In this paper we extend these connections to random regular graphs for all ferromagnetic models. Using these connections, we establish the Bethe prediction for every ferromagnetic spin system on random regular graphs, which says roughly that the expectation of the log of the partition function $Z$ is the same as the log of the expectation of $Z$. As a further consequence of our results, we prove for the ferromagnetic Potts model that the Swendsen--Wang algorithm is torpidly mixing (i.e., exponentially slow convergence to its stationary distribution) on random $\Delta$-regular graphs at the critical temperature for sufficiently large $q$.