How does spin arise in (relativistic or not) quantum mechanics?
What are particles in the first place? And what precisely is this property that we call "spin"?
A modern way of arriving at the notion of particles, that may be more transparent than the historical version you summarize in your post, is the way they are introduced in Weinberg's textbook on Quantum Field Theory.
Take the Hilbert space $\mathcal{H}$ of your quantum theory, and take the group $G$ of symmetries of your spacetime. To ensure that an experiment will give the same results if I move it around in spacetime, rotate it, or take it with me on a cruise, there should exist a unitary representation $U$ of $G$ on $\mathcal{H}$. Ie. for any $g \in G$, there should exist a unitary operator $U(g)$ on $\mathcal{H}$, with:
$$
U(1) = \text{id}_{\mathcal{H}} \;\&\; U(g.h) = U(g) U(h).
$$
We can then decompose this representation $\mathcal{H}, U$ into simpler representations. This means writing $\mathcal{H}$ as a direct sum of smaller Hilbert spaces:
$$
\mathcal{H} = \bigoplus_k \mathcal{H}_k
$$
with each $\mathcal{H}_k$ being stabilized by all $U(g)$ for $g \in G$, so that $\mathcal{H}_k, U_k := \left. U \right|_{\mathcal{H}_k}$ is itself a unitary representation of G.
A representation that does not contain any smaller representation is called an irreducible representation, or irrep, and the simplest irreps are the ones that hold quantum states with just one particle (refer to Weinberg for what "simplest" exactly means here). So each (simple) irrep of the symmetry group $G$ that can be found in our theory $\mathcal{H}$ is what we define as a particle species.
Integer spins
Good, so if we want to know which particle species are physically possible, we just need to know what are the "simplest" representations of our group of symmetries. Fortunately, a full classification of those is known for the Poincaré or Galilean group. The way it is constructed would be too long to be reproduced here, but again it can be found in great details in Weinberg for the Poincaré group (brief accounts of both the Poincaré and the Galilean case can be found on wikipedia). The bottom line is that, for physically-admissible massive particles, they have the form:
$$
\mathcal{H}_k = \text{Span} \left\{ \left| \vec{p}, m \right\rangle \middle| \vec{p} \in \mathbb{R}^3, m \in \mathbb{Z}, -s_k \leq m \leq +s_k \right\}
$$
with the non-negative integer $s_k$ being what is called the spin of this particle species $k$ and $\vec{p}$ being its impulsion. The impulsion determines how the particle transform under a spatial translation:
$$
U_k(\text{translation by } \vec{\alpha}) \left| \vec{p}, m \right\rangle = e^{i \vec{\alpha}.\vec{p}} \left| \vec{p}, m \right\rangle
$$
(just as in the good old momentum representation of QM). Its spin number determines how it transforms under a rotation:
$$
U_k(\text{rotation by } \vec{\theta}) \left| \vec{p}, m \right\rangle = \sum_{m^\prime} R^{(s_k)}_{mm^\prime}(\vec{\theta}) \left| R(\vec{\theta}) \vec{p}, m^\prime \right\rangle
$$
with $R(\vec{\theta})$ the usual $3\times 3$ rotation matrix acting on $\vec{p}$ and $R^{(s_k)}$ an irrep of the rotation group $\mathcal{SO}(3)$.
The intuition behind this form of $U_k$ is that the spin $s_k$ captures the way the particle may be affected by a rotation beyond the obvious rotation of its impulsion $\vec{p}$. The classical analogy here is that of a rigid body, which changes not only its position but also its orientation under a rotation.
The reason the spin is an integer is because the irreps of $\mathcal{SO}(3)$ are labeled by integers: spin-0 is the trivial representation $R^{(0)}(\vec{\theta}) = \text{id}, \forall \vec{\theta} \in \mathbb{R}^3$ on a 1-dimensional vector space, spin-1 is the usual representation by $3 \times 3$ matrices, and so on.
Half-integer spins
But now there is a twist (figuratively and mathematically...). As mentioned in an earlier comment by ACuriousMind, and as explained in great details in the linked thread, the overall phase of a quantum state is not physically measurable. This means that we can get away with less that a strict unitary representation of $G$ on $\mathcal{H}$, and still ensure that all experimental results are invariant under $G$! Namely, we can replace:
$$
U(g\cdot h) = U(g) U(h) \text{ by } U(g\cdot h) = e^{i \varphi(g,h)} U(g) U(h)
$$
with the phase factors $\varphi(g,h)$ satisfying suitable consistency relations. Such unitary representations "up to additional phase factors" are called projective representations.
If one goes through the math, we find that for the Poincaré/Galilean group this gives a few additional possible irreps, corresponding to particle species with half-integer spin. They correspond to projective representations of the rotation group in which a rotation of $2\pi$ has a non-trivial (albeit non-detectable) action on the quantum state:
$$
U_k(\text{rotation by } 2\pi) \left| \vec{p}, m \right\rangle = - \left| \vec{p}, m \right\rangle
$$
Observable signature?
But wait! If this extra minus sign is not physically detectable anyway, how do we know that some particles have half-integer spin?
This has to do with the properties of quantum measurments, which will reveal the spectrum (aka. eigenvalues) of the measured observable. We cannot directly observe that the quantum state vector transforms under a projective representation but we can determine it indirectly because it gets imprinted in the spectrum of the angular-momentum operator.
What about classical (non-quantum) mechanic?
Take a classical mechanical system, say, for concreteness, a system of rigid bodies possibly interacting via conservative forces. The phase space of such a system carries a non-projective representation $T$ of the Galilean group (we can check this by writing it explicitly). But this representation is not a linear representation (at best it may be an affine representation, since translations act, well, by translations). So spin in the above sense does not immediately make sense.
Instead, we can do classical statistical physics for this system: ie. write a field equation for a probability distribution $\rho$ on the phase space (which can be seen as the classical counterpart of a quantum mechanical wave-function). The space $\mathcal{P}$ of such probability distributions naturally carries a linear representation $U$ of $G$ defined by:
$$
\forall g \in G,\; [U(g) \rho](x) = \rho\big(T(g^{-1})x\big)
$$
which is, again, a non-projective representation (strictly speaking, admissible probability distributions are positive and normalized, but we can study their spin properties by working in the vector space they span: this is analogous to considering the whole Hilbert space in quantum mechanics, although actual quantum states should be normalized).
So, what would a "half-integer-spin mode" for such a system be? According to the previously explained definition of spin, that would be a half-integer-spin irrep $\mathcal{P}_k \subseteq \mathcal{P}, U_k := \left. U \right|_{\mathcal{P}_k}$ appearing in the decomposition of $\mathcal{P}, U$. Can such a $\mathcal{P}_k$ exist? No!
Indeed, if it would, we would have a distribution $\rho \in \mathcal{P}_k \setminus \{0\}$ such that
$$
U(\text{rotation by } 2\pi) \rho = U_k (\text{rotation by } 2\pi) \rho = -\rho,
$$
but, since $U$ is a non-projective representation, we already know that $U(\text{rotation by } 2\pi) \rho = \rho$.
A similar argument can be applied for example to the classical electromagnetic field: the space of solutions of Maxwell's equations carries a non-projective linear representation of the Poincaré group (one could say: by historical definition of the latter).
What about a thermodynamical system?
Suppose I take a large number of mechanical bodies interacting via conservative forces (say molecules) and take the thermodynamical limit to derive effective equations for some macroscopic variable (eg. their density). Could such an equation exhibit half-integer-modes? Ie. could its space of solutions carry a projective representation of G? Let us do some though experiment:
Take two rigorously identical boxes containing this thermodynamical system and perform the exact same experiment on them, except that the second one is first subjected to a full $2\pi$-rotation (very slowly, so as to not perturb any (local) thermodynamical equilibrium). Because the underlying microscopic theory carries a non-projective representation of G, the two experiments should give the exact same result!
Note that in arguments of this kind, one has to be very careful. The thermodynamic limit can do funny things to the symmetries of a systems. This is known as symmetry-breaking: while the space of solutions of the underlying microscopic theory may be invariant under a certain group $G$, a given thermodynamical phase may have less symmetry because it fails to explore the full solution space (keyword: ergodicity, or more precisely lack thereof).
But, such a mechanism cannot turn a non-projective representation into a projective one: since a $2\pi$-rotation brings me back on the exact same microscopic configuration from which I started, I am guaranteed not to land in a different thermodynamical phase.
Can we cheat?
Suppose I come up with a mathematical description of some physically valid classical system in which, for technical reasons, I choose to introduce some auxiliary, non-measurable quantity (eg. a complex phase). Since the auxiliary quantity is non-measurable, I can let it transform in whatever way is mathematically convenient. In this way, I may arrive at a description of a classical system which exhibits a projective representation.
But still the original physical system will not exhibit any observable half-integer-spin behavior. As the truly measurable quantities have to be invariant under a full $2\pi$-rotation, there should exist a basic description of the same system, that refrains from introducing any auxiliary quantities, and carries a non-projective representation. Computing experimental predictions using this basic description, no half-integer-spins should show up.
TL;DR: This is crucially differently from the above discussed quantum mechanical case, in which you can hide a non-projective representation, so as to preserve $2\pi$-rotation-invariance, while nevertheless retaining some observable signature.
Bonus: Does Wick-rotating a quantum equation give a thermodynamical equation?
I do not think Wick rotation should be thought as some kind of magic transformation to turn a QM equation into a thermodynamical one.
There is a connection between quantum field theory on 3+1d Minkowski spacetime and statistical field theory in 4d Euclidean space. But statistical physics (the study of the probability distribution over (field) configurations) is not quite the same as thermodynamics (the derivation of effective equations for macroscopic variables in the large-number-of-particles limit).
I suspect the appearance of the heat equation as the complex-time Schrödinger equation is more a coincidence, coming from the fact that, well, there are only so many linear PDEs you can write with a certain order in space and time derivatives.
If you would like to investigate the Wick rotated Dirac equation anyway, I guess a good place to start would be the Wick-rotated gamma matrices. You will get a field equation carrying a projective representation of the 4d Euclidean group, sure. But Wick-rotating a physically valid quantum equation does not a priori guarantee any particular physical relevance for the resulting equation: in fact, such an equation cannot describe any actual physical system, if only because, as pointed out by flippiefanus, we do not live in 4d Euclidean space ;-).
