Why do theories of nature prefer to be renormalizable and not super-renormalizable? It seems to me (correct me if I am wrong) that all theories in the Standard Model are exactly renormalizable, as opposed to non-renormalizable or super-renormalizable. In a sense, we could say that these theories live at a critical point, that is on the boundary between non-renormalizable and super-renormalizable theories. It is kind of clear why a theory of nature would prefer to be renormalizable than non-renormalizable, but is there a known reason (or conjecture) for why nature would prefer to be renormalizable than super-renormalizable? Could that not be a hint of where to search for theories?
(of course an important caveat would be that renormalization may arise as a flaw of our current theoretical framework)
 A: It is important to understand that the modern answer to "why renormalizable" has nothing to do with consistency requirements.
The confusion comes from history: in the early days of QFT, people thought that the perturbative expansion can be used as the definition of the full theory, provided that it is renormalizable. Here, renormalizability is desired, because it means the description of the theory is complete. I.e. if you want to unify General Relativity and the Standard Model into a single Theory of Everything, that ToE better be a complete description of all physical phenomena if it is to live up to its name.
It has become apparent later that even renormalizable (in fact, even superrenormalizable) perturbative expansions do not define the full theory. They are not approximation schemes in the usual meaning of that word. That is, they can not approximate a $n$-point function to an arbitrary accuracy.
For a convergent series, e.g. the Taylor series of an analytic function, the more terms you sum – the better your approximation to the right answer. This can be mathematically written down as an infinite sum:
$$
f(x) = \sum_{n=0}^{\infty} \frac{1}{n!} f^{(n)} x^n.
$$
The sum in the r.h.s. doesn't always converge on the entire real axis, however. Take, for example,
$$
\frac{1}{1-x} = \sum_{n=0}^{\infty} x^n.
$$
Interestingly, the l.h.s. is defined for $x \neq 1$, but the r.h.s. only converges for $-1 < x < 1$ – a completely different domain! The number $1$ is called the radius of convergence (akin to the fact that on the complex plane, the series is convergent when $|x| < 1$).
In QFT, the perturbative expansions are in the powers of the coupling constant, $g$:
$$
A(g) = \sum_{n=0}^{\infty} A_n g^n.
$$
Here $A$ is some quantity that should exist in the theory, and $A_n$ are the coefficients in this approximation. $A_n$ are usually represented by sums of Feynman diagrams with $n$ interaction vertices (each vertex contributes a single power of $g$).
So the natural question to ask is: what is the radius of convergence of $A(g)$?
The answer can feel very strange unless you already know it: we believe that almost uniformly, for all interacting QFTs in 2 or more space-time dimensions, the radius of convergence is zero!
Among other things, that means that the perturbative expansion can never define the full theory, and it can never be complete. Even if we have a renormalizable or a superrenormalizable perturbative expansion, summing it doesn't give you the answer, in fact, it gives you infinity.
Such expansions are known as asymptotic expansions. Their defining property is that for any chosen order of approximation $N$, and for any chosen accuracy $a$, there exists a value $g_{\text{max}}$ such that
$$ | A(g) - \sum_{n=0}^N A_n g^n | < a $$
for all $g < g_{\text{max}}$.
In this sense, asymptotic expansions can approximate exact quantities. In fact, in practice, rather paradoxically, asymptotic expansions are much better at approximating functions than convergent expansions are! (which is the reason why perturbative QFT is so successful)
But this notion is very different from the usual normal of convergence. In normal convergence, for any accuracy and any value of $g$ (within the radius of convergence) there exists an order of approximation that is good enough to yield that accuracy. Here, for any accuracy and any order, there exists a value of $g$.
In QFT, this means that even renormalizable, even superrenormalizable perturbative theories can't be used to produce arbitrarily precise predictions. What happens is that you can only achieve accuracy $a$ if your coupling constant is small enough. But the coupling constant is varied only on paper – in nature, it must have a fixed value that is a parameter of the theory. So either you're in luck and the value of the coupling allows to approximate the answer with accuracy $a$, or you're not! And once you start making $a$ smaller and smaller, eventually the physical value of the coupling will become too large for the asymptotic expansion to be a good approximation – the perturbative series will blow up. This establishes that there is a fundamental accuracy limit for all perturbative QFTs beyond which they are unable to make any predictions.
The symptoms (loss of predictive power) are very similar to what happens with nonrenormalizable theories, except that the loss of predictive power due to the diverging asymptotic series is common to all perturbative QFT, renormalizable or not.
Therefore, perturbative QFT can't be the final formulation of physics, even if it is renormalizable. One can get around this by studying non-perturbative QFT (which are notoriously difficult to even define mathematically, let alone solve / extract predictions). Another possible way is to give up on QFT alltogether and treat perturbative QFT as a convenient approximation to something very different, like strings, discrete space time of Loop Quantum Gravity, etc. In fact, because the kinematics of General Relativity lives on a differential manifold rather than on the Minkowski space-time (a property which is known as background independence); it does look very unlikely that the ToE is a Minkowski space QFT, perturbative or nonperturbative.
Regardless of which path one chooses for the formulation of the fundamental theory, perturbative QFT, being an extremely good approximation, addresses many important issues. One of these issues is the question of why in nature we only observe renormalizable interactions. The answer to that question has nothing to do with the original, naive notion of "renormalizable = consistent", which turned out to be a poor guide into the structure of the QFT.
Instead, it turns out that whatever interactions the fundamental theory has at high energies, only those that give rise to renormalizable QFT operators are physically measurable at low energies.
This conceptual picture is captured mathematically by the Renormalization Group. Having given up on formulating the perturbative QFT as a fundamental theory, one adopts a pragmatic attitude: the perturbative QFT is defined with a cutoff $\Lambda$ that makes all its integrals finite, and hence removes the need for getting rid of infinities. However, it doesn't remove the need to renormalize. Renormalization isn't about stuffing infinities under the carpet, like Dirac once said – it turns out to be a very down to Earth procedure that arises naturally. What happens is – once you have a theory with a cutoff that is finite, you want to make predictions with it. For example, you're trying to measure the mass of one of your theory's particle species. But it turns out that this mass is far from the value that you've put by hand into the Lagrangian. It acquires additional contributions, which frequently outweigh the original value by many orders of magnitude, raising concerns from the point of view of naturalness. So the predicted value of the mass and the value you put in the theory by hand are very different. This of course means that the value you put by hand has to be changed, so that the predicted value matches the experimentally measured.
Once you fix the predicted value, you have to make the bare value (the one you put in by hand) depend in a non-trivial way on the cutoff $\Lambda$ in order to retain the correct prediction for all values of $\Lambda$. This is the first example of the renormalization group flow.
It may appear at the first sight that this flow is unphysical and unmeasurable (we're only tweaking a parameter to match observations which are fixed and aren't flowing). There's some truth to it (and in fact there are different definitions of the RG flow which correspond to different renormalization schemes, which signals that some aspects of that flow are not observable), however, there is an important physical aspect of the renormalization group flow which is – anomalous scaling behavior.
Since the theory contains a dimensionful parameter $\Lambda$, it can be non-invariant under scaling transformations, even in the limit where $\Lambda \rightarrow \infty$! In fact, that's exactly what happens in many real theories, including Yang-Mills. We say that the scaling symmetry acquires a quantum anomaly. This means a non-trivial behavior of coupling constants under scaling, known as the anomalous dimension or the beta function.
The dominant term in the scaling law of couplings is usually the classical term,
$$ g \rightarrow b^d g, $$
with $b$ the scaling parameter, and $d$ the classical dimension of the coupling.
Couplings with $d > 0$ are called relevant. In the infrared region (which we can pass to by applying a scaling transformation with a very large $b$) they acquire large values, and are physically observable.
By a divergence index counting argument, these couplings correspond to superrenormalizable interactions. Note that this is a pure coincidence – our analysis doesn't use the renormalizability!
Couplings with $d < 0$ are called irrelevant, because in the infrared region they acquire very small values and are unobservable.
By a divergence index counting argument, these correspond to nonrenormalizable interactions.
In order for the analysis above to hold, a very important assumption has to be valid: the classical behavior of the scaling law has to be the leading contribution in the renormalization group. That translates into inequalities on coupling parameter values. This assumption does not always hold!
Probably the most important example of a situation where the classical term is not a leading term are marginal couplings – those with $d = 0$. For them, the classical term vanishes, and the scaling behavior of such couplings is completely determined by the quantum effects of the renormalization group. Depending on many intricate properties like the particle content of the theory, these effects can either render a marginal coupling observable in the infrared region, or unobservable.
This is the real reason we only see superrenormalizable and renormalizable interactions in nature, according to the modern understanding. Not a consistency argument, but merely the fact that if there was a nonrenormalizable interaction (in fact, there probably is!) – we would completely miss it due to its value becoming very small at low energies.
In QCD, as long as we don't couple too much matter to it, the property called asymptotic freedom holds: QCD's marginal coupling behaves a lot like a relevant coupling (large in the IR, small in the UV), though the scaling law is logarithmic rather than polynomial (since there is no classical contribution).
Interestingly, when the scale of the experiment reaches $\Lambda_{QCD}$ (about 200 MeV), the coupling constant becomes comparable to $1$ and the asymptotic perturbative expansion stops giving good approximations. QCD undergoes a phase transition at that point: quarks confine into color-neutral hadrons. Perturbative QFT completely breaks down below $\Lambda_{QCD}$, not only in theory, but also in practice.
