How to motivate the skein relations?

Alexander realized they were useful, then Conway. However, Jones clearly was the one who really made a big bang with a skein relation. This allowed him to see a connection between the Jones polynomial and state sums in statistical mechanics. This was followed by HOMFLYPT, which might be the first time a skein relation was used to define an invariant rather than encapsulate some of its properties. I would say that Kauffman really simplified the study of the Jones polynomial by way of his bracket relation. Witten used the skein relation to build a hypothetical connection with TQFT. This connection was made rigorous by Reshitkhin and Turaev from a quantum groups perspective. The work of Habeger, Masbaum, Vogel, and Blanchet made it rigorous from a skein theoretic viewpoint. The connection between skein relations and Lie groups probably appeared first in the work of Turaev and Wenzl, where they classified skein relations by what family of Lie groups you were working with. There was work of Blanchet and Beliakova that built on this, especially understanding B-type Lie algebras. Xiao-Song Lin built a connection between the Jones polynomial and finite type invariants with skein relations. Effie Kalfagianni used them to extend the Jones polynomial to knots in other three-manifolds. Bar-Natan used skein relations to powerful effect to build connections between finite type invariants and Lie algebras. Bullock's work was the first to build the direct connection between skein relations and trace identities. This was built on in the work of Sikora who did it for many more Lie groups. Kuperberg's spiders were also actually about skein relations, and his work is starting have impact in the study of flag varieties. The skein relation short exact sequence of chain complexes first appeared in the work of Khovanov, though it had been simmering in Floer theory for a long time in the related surgery triangle.

I thought I would give a motivating example.

Recall that matrices satisfy their characteristic polynomial. If $A$ is a two by two matrix of determinant one, then its characteristic polynomial is $\lambda^2-tr(A)\lambda +1$. Hence we know $A^2-(tr(A)A+Id=0$ where the zero on the right is a two by two matrix of zeroes. Multiply this through by $A^{-1}$ to symmetrize it as $A+tr(A)Id +A^{-1}=0$. Now multiply by any two by two matrix $B$ to get $BA-tr(A)B+BA^{-1}=0$. Finally, take the trace to get $tr(BA)-tr(A)tr(B)+tr(BA^{-1})=0$. This is the fundamental trace identity for $SL_2(\mathbb{C})$. Let $X_i$ be variables which you can think of as $2\times 2$ matrices. Form all words in the $X_i$ and then take all polynomials in traces of these words. Every identity between these polynomials is a consequence of $tr(A)=tr(A^{-1})$, $tr(AB)=tr(BA)$, $tr(Id)=2$, and the fundamental trace identity.

Let $M$ be a topological space. Take polynomials on free homotopy classes of loops in $M$ and mod out by the relations coming from replacing a loop by its inverse, replacing the null homotopic loop by $-2$, and the Kauffman bracket skein relaton with $A$ set equal to $-1$. In this setting the Kauffman bracket skein relation is the fundamental trace identity.

This ring is the coordinate ring of the unreduced affine scheme of characters of $SL_2(\mathbb{C})$ representations of the fundamental group of $M$.

If $M$ is a surface then you can fatten the surface up to $M\times [0,1]$ and instead use framed links in $M\times [0,1]$. You multiply by stacking framed links. Use the Kauffman bracket skein relations at $A=e^h$. This algebra is a quantization of the $SL_2(\mathbb{C})$ representations of the fundamental group of the surface with respect to Goldman's Poisson structure. That fact is established by proving that Goldman's Poisson bracket can be computed skein theoretically.

A good starting point for the above material might be the article "Understanding the Kauffman Bracket Skein Relation" by Bullock, F, and Kania-Bartoszynska.

The importance of skein relations is that it allows you to make teleological connections between knot polynomials and other, seemingly distant areas of mathematics like representation theory, quantization, algebraic geometry, gauge theory, and low dimensional geometry and topology.

One of the earliest appearances of the ingredients for a skein relation can be found in Romilly Allen's 1904 book on Celtic Knotting. There he explains that designers of Celtic knot patterns first start with a regular weave and then systematically replace some of the crossings by "horizontal" and "vertical" smoothings. These days, it is easy to find working demonstrations of this process on the internet. The skein triple consists in an unoriented crossing and its two possible smoothings. This is the skein triple underlying the bracket polynomial model for the Jones polynomial and directly related to Temperley-Lieb algebra. Of course neither Romilly Allen nor the Celtic knotters had any intention to calculate polynomial invariants of knots! That came much later. Another important point about this simplest skein relation is that it is related (by shading the link diagram) to the contraction and deletion operations on the edges of a plane graph. This is why the bracket polynomial and the Jones polynomial are related to the Dichromatic (Tutte) polynomial and to signed generalizations of the Tutte polynomial. And then finally, the two smoothings of a crossing are related to one another in that they form the 'top' and 'bottom' of a saddle surface that joins them. This is the key geometric idea behind Khovanov homology.

Regarding "when", it was Alexander, in his paper on what we call the Alexander polynomial. Conway was the first to popularize them, I believe.

Why are they useful? I'm not sure I believe they are so useful. Sometimes I'm interested in computing Alexander invariants but the knots and links that I'm looking at do not have easy-to-compute diagrams associated to them. Say you have a homology 3-sphere and you apply the JSJ-decomposition to it. This produces some knots and links in homology spheres but working out diagrams can be a pain. So sometimes it's far easier to use the covering space definition to get at the Alexander polynomials. Moreover, Skein relations don't give you the Alexander module and how Poincare duality works on the module, or naturality (when you have maps between 3-manifolds), etc.