The solution comes from assuming that
$$
x(t) = x[k] + \frac{x[k+1]-x[k]}{T}(t-kT), \ \ t\in[kT,(k+1)T)
$$
where $x[k]:=x(kT)$, this is, a first order hold or a linear interpolation of $x[k]$ and $x[k+1]$ in between $kT$ and $(k+1)T$. Now, we simply compute the EXPLICIT solution of the differential equation (at instants $t=kT$) by integrating it from $t=kT$ towards $t=(k+1)T$. This should be enough to derive the result but let me get into more detail. First, multiply the ODE by $\frac{1}{\tau}e^{t/\tau}$ to obtain:
$$
\frac{dy}{dt}e^{t/\tau} + \frac{1}{\tau}e^{t/\tau}y(t) = \frac{1}{\tau}e^{t/\tau}x(t)
$$
and notice that $\frac{d}{dt}\left(e^{t/\tau}y(t)\right) = \frac{dy}{dt}e^{t/\tau} + \frac{1}{\tau}e^{t/\tau}y(t)$ and thus:
$$
\frac{d}{dt}\left(e^{t/\tau}y(t)\right) = \frac{1}{\tau}e^{t/\tau}x(t)
$$
Now integrate both sides:
$$
\int_{kT}^{(k+1)T}\frac{d}{dt}\left(e^{t/\tau}y(t)\right)dt = \int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt
$$
the left hand side results (by fundamental theorem of calculus) into:
$$
\int_{kT}^{(k+1)T}\frac{d}{dt}\left(e^{t/\tau}y(t)\right)dt = \left.e^{t/\tau}y(t)\right|_{kT}^{(k+1)T} =e^{(k+1)T/\tau}y[k+1] - e^{kT/\tau}y[k]
$$
Thus,
$$
e^{(k+1)T/\tau}y[k+1] = e^{kT/\tau}y[k] + \int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt
$$
Now, simply divide by $e^{(k+1)T/\tau}$ to obtain
$$
y[k+1] = e^{-T/\tau}y[k] + e^{-(k+1)T/\tau}\int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt
$$
As you can see we are almost done. We just have to compute:
$$
\begin{aligned}
&e^{-(k+1)T/\tau}\int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt =\\ &e^{-(k+1)T/\tau}\int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}\left(x[k] + \frac{x[k+1]-x[k]}{T}(t-kT)\right)dt
\end{aligned}
$$
which we can do explicitly since $x[k]$ and $x[k+1]$ are just constants. After some algebra you can factor terms with $x[k]$ and terms with $x[k+1]$ and you will obtain precisely what you wanted. If you want to check, I plugged this into wolfram alpha here ("x1" is $x[k]$ and "x2" is $x[k+1]$).
This kind of discretizations are nice since, if $x(t)$ does behave as a piece-wise linear continuous function, the discretization is NOT an approximation, but an equivalent (and exact) representation of the ODE, this is: since it was obtained by explicit computation of the solution of the ODE, no information is lost about the output samples $y[k]$. (in contrast to approximating $dy/dt \approx (y[k+1]-y[k])/T$ or some other approximation as the bilinear one)
NEW REMARK: I just want to point out that this same procedure can be performed for other high order systems, or even for other types of holders. Recall that any LTI system with scalar input $x$ and scalar output $y$ can be represented by its transfer function, or equivalently by a state space representation
$$
\frac{dz}{dt} = Az + Bx,\ \ \ y = Cz
$$
with $z(t)\in\mathbb{R}^n$ with $n$ the order of the system, $A$ a square matrix, $B$ a column vector and $C$ a row vector. Recall that the poles of the transfer function are the eigenvalues of $A$ and zeros of the transfer function depend only on the coefficients on contained in $C$. Thus, we can solve for the system on $y[k]$ explicitly using a similar procedure as before:
$$
z[k+1] = e^{AT}z[k] + e^{A(k+1)T}\int_{kT}^{(k+1)T}e^{-A\tau}Bx(\tau)d\tau, \ \ y[k]=Cz[k]
$$
where $e^{A\tau}$ is the exponential matrix. Finally, plugging in the particular form of $x(t)$ results in the final form of the discretization (in discrete time state space representation). This can be used to derive (exact or "consistent") discretizations for ZOH, FOH, high order holders or more exotic holders.