[proofplan]
We prove the recurrence directly from the subset expansion of the Tutte polynomial. The sum over all subsets of $E$ is split into subsets not containing $e$ and subsets containing $e$, and subsets in the second class are written uniquely as $A\cup\{e\}$ with $A\subset E\setminus\{e\}$. The rank identities for deletion and contraction then identify the two pieces with the Tutte polynomials of $M\setminus e$ and $M/e$, with the loop and coloop cases producing the extra factors $y$ and $x$.
[/proofplan]
[step:Split the subset expansion according to whether the subset contains $e$]
Let $E_0:=E\setminus\{e\}$. By the subset expansion of the Tutte polynomial,
\begin{align*}
T_M(x,y)=\sum_{S\subset E}(x-1)^{r_M(E)-r_M(S)}(y-1)^{|S|-r_M(S)}.
\end{align*}
Every subset $S\subset E$ either has the form $S=A$ with $A\subset E_0$, or has the form $S=A\cup\{e\}$ with $A\subset E_0$. Therefore
\begin{align*}
T_M(x,y)=\sum_{A\subset E_0}(x-1)^{r_M(E)-r_M(A)}(y-1)^{|A|-r_M(A)}+\sum_{A\subset E_0}(x-1)^{r_M(E)-r_M(A\cup\{e\})}(y-1)^{|A|+1-r_M(A\cup\{e\})}.
\end{align*}
Call the first sum $P_0(x,y)$ and the second sum $P_1(x,y)$.
[guided]
The subset expansion is the right starting point because deletion and contraction are operations on subsets of the ground set. Define $E_0:=E\setminus\{e\}$. Then every subset of $E$ is accounted for exactly once in one of two ways: either it does not contain $e$, in which case it is a subset $A\subset E_0$, or it contains $e$, in which case it is uniquely $A\cup\{e\}$ for a subset $A\subset E_0$.
Thus
\begin{align*}
T_M(x,y)=\sum_{S\subset E}(x-1)^{r_M(E)-r_M(S)}(y-1)^{|S|-r_M(S)}.
\end{align*}
Splitting the indexing set $2^E$ into the disjoint union of subsets avoiding $e$ and subsets containing $e$ gives
\begin{align*}
T_M(x,y)=\sum_{A\subset E_0}(x-1)^{r_M(E)-r_M(A)}(y-1)^{|A|-r_M(A)}+\sum_{A\subset E_0}(x-1)^{r_M(E)-r_M(A\cup\{e\})}(y-1)^{|A|+1-r_M(A\cup\{e\})}.
\end{align*}
The first sum records subsets that survive unchanged in deletion. The second sum is written with $A\cup\{e\}$ because contraction is naturally described by comparing the rank of $A\cup\{e\}$ in $M$ with the rank of $A$ in $M/e$.
[/guided]
[/step]
[step:Identify the loop case with a factor of $y$]
Assume that $e$ is a loop. Then $r_M(\{e\})=0$, and for every $A\subset E_0$,
\begin{align*}
r_M(A\cup\{e\})=r_M(A).
\end{align*}
Also $r_M(E)=r_M(E_0)$. The rank function of $M\setminus e$ is the restriction of $r_M$ to subsets of $E_0$, so
\begin{align*}
P_0(x,y)=T_{M\setminus e}(x,y).
\end{align*}
For the second sum, the same rank identities give
\begin{align*}
P_1(x,y)=\sum_{A\subset E_0}(x-1)^{r_M(E_0)-r_M(A)}(y-1)^{|A|+1-r_M(A)}.
\end{align*}
Factoring out one power of $y-1$ gives
\begin{align*}
P_1(x,y)=(y-1)T_{M\setminus e}(x,y).
\end{align*}
Therefore
\begin{align*}
T_M(x,y)=P_0(x,y)+P_1(x,y)=T_{M\setminus e}(x,y)+(y-1)T_{M\setminus e}(x,y)=yT_{M\setminus e}(x,y).
\end{align*}
[/step]
[step:Identify the coloop case with a factor of $x$]
Assume that $e$ is a coloop. Then $r_M(\{e\})=1$, $r_M(E)=r_M(E_0)+1$, and for every $A\subset E_0$,
\begin{align*}
r_M(A\cup\{e\})=r_M(A)+1.
\end{align*}
For contraction by a non-loop element, the rank function of $M/e$ on subsets $A\subset E_0$ is
\begin{align*}
r_{M/e}(A)=r_M(A\cup\{e\})-r_M(\{e\}).
\end{align*}
Hence, since $e$ is a coloop,
\begin{align*}
r_{M/e}(A)=r_M(A).
\end{align*}
Also
\begin{align*}
r_{M/e}(E_0)=r_M(E)-1.
\end{align*}
Substituting these identities into the two sums gives
\begin{align*}
P_0(x,y)=(x-1)T_{M/e}(x,y).
\end{align*}
Similarly,
\begin{align*}
P_1(x,y)=T_{M/e}(x,y).
\end{align*}
Adding the two pieces,
\begin{align*}
T_M(x,y)=P_0(x,y)+P_1(x,y)=(x-1)T_{M/e}(x,y)+T_{M/e}(x,y)=xT_{M/e}(x,y).
\end{align*}
[guided]
The coloop case is dual in spirit to the loop case, but the extra factor appears in the $x$-exponent rather than the $y$-exponent. Since $e$ is a coloop, adjoining $e$ raises the rank of every subset of $E_0$ by exactly one:
\begin{align*}
r_M(A\cup\{e\})=r_M(A)+1
\end{align*}
for every $A\subset E_0$. Also the total rank satisfies
\begin{align*}
r_M(E)=r_M(E_0)+1.
\end{align*}
The contraction rank function is defined, for $A\subset E_0$, by
\begin{align*}
r_{M/e}(A)=r_M(A\cup\{e\})-r_M(\{e\}).
\end{align*}
Because a coloop is not a loop, $r_M(\{e\})=1$. Therefore
\begin{align*}
r_{M/e}(A)=r_M(A\cup\{e\})-1=r_M(A).
\end{align*}
In particular,
\begin{align*}
r_{M/e}(E_0)=r_M(E)-1.
\end{align*}
Now evaluate the first piece $P_0(x,y)$. For $A\subset E_0$, the exponent of $x-1$ is
\begin{align*}
r_M(E)-r_M(A)=r_{M/e}(E_0)+1-r_{M/e}(A).
\end{align*}
The exponent of $y-1$ is
\begin{align*}
|A|-r_M(A)=|A|-r_{M/e}(A).
\end{align*}
Thus
\begin{align*}
P_0(x,y)=(x-1)T_{M/e}(x,y).
\end{align*}
For the second piece, the exponent of $x-1$ is
\begin{align*}
r_M(E)-r_M(A\cup\{e\})=r_M(E)-r_M(A)-1=r_{M/e}(E_0)-r_{M/e}(A).
\end{align*}
The exponent of $y-1$ is
\begin{align*}
|A|+1-r_M(A\cup\{e\})=|A|+1-r_M(A)-1=|A|-r_{M/e}(A).
\end{align*}
Therefore
\begin{align*}
P_1(x,y)=T_{M/e}(x,y).
\end{align*}
Combining the two contributions gives
\begin{align*}
T_M(x,y)=(x-1)T_{M/e}(x,y)+T_{M/e}(x,y)=xT_{M/e}(x,y).
\end{align*}
[/guided]
[/step]
[step:Identify the ordinary deletion and contraction terms when $e$ is neither a loop nor a coloop]
Assume that $e$ is neither a loop nor a coloop. Since $e$ is not a coloop, deleting $e$ does not lower the total rank, so
\begin{align*}
r_M(E)=r_{M\setminus e}(E_0).
\end{align*}
Since the deletion rank function is the restriction of $r_M$ to subsets of $E_0$, for every $A\subset E_0$,
\begin{align*}
r_{M\setminus e}(A)=r_M(A).
\end{align*}
Hence
\begin{align*}
P_0(x,y)=T_{M\setminus e}(x,y).
\end{align*}
Since $e$ is not a loop, $r_M(\{e\})=1$. For every $A\subset E_0$, contraction has rank function
\begin{align*}
r_{M/e}(A)=r_M(A\cup\{e\})-1.
\end{align*}
Also
\begin{align*}
r_{M/e}(E_0)=r_M(E)-1.
\end{align*}
Therefore
\begin{align*}
r_M(E)-r_M(A\cup\{e\})=r_{M/e}(E_0)-r_{M/e}(A)
\end{align*}
and
\begin{align*}
|A|+1-r_M(A\cup\{e\})=|A|-r_{M/e}(A).
\end{align*}
Substitution in $P_1(x,y)$ gives
\begin{align*}
P_1(x,y)=T_{M/e}(x,y).
\end{align*}
Thus
\begin{align*}
T_M(x,y)=P_0(x,y)+P_1(x,y)=T_{M\setminus e}(x,y)+T_{M/e}(x,y).
\end{align*}
[/step]
[step:Conclude the recurrence in all three cases]
The three mutually exclusive cases for an element $e\in E$ are: $e$ is a loop, $e$ is a coloop, or $e$ is neither a loop nor a coloop. The preceding steps prove respectively
\begin{align*}
T_M(x,y)=yT_{M\setminus e}(x,y),
\end{align*}
\begin{align*}
T_M(x,y)=xT_{M/e}(x,y),
\end{align*}
and
\begin{align*}
T_M(x,y)=T_{M\setminus e}(x,y)+T_{M/e}(x,y).
\end{align*}
This is exactly the deletion-contraction recurrence for the Tutte polynomial.
[/step]