Why do we use this particular approximation for the bilinear transform?
The Forward Euler Transform $$z=e^{sT}\approx1+sT\leftrightarrow s\approx(z-1)/T=\frac{1-z^{-1}}{Tz^{-1}}$$ is easy to understand in that it is a direct translation and scaling from the \$s\$-domain to the z-domain. But the translation can transform stable \$s\$-domain poles into unstable \$z\$-domain poles.

To see this consider the diagram below.

The Left Half Plane in the \$s\$-domain (shaded) is scaled by \$T\$ and translated by \$1\$ to the \$z\$-domain. It should be clear to see that a pole X that is stable in \$s\$ can be unstable in \$z\$ by the Forward Euler transform.



In contrast, the bilinear transform $$s\approx\frac2T\frac{z-1}{z+1}=\frac2T\frac{1-z^{-1}}{1+z^{-1}}$$ translates the entire LHP of the s-domain in the unit circle of the \$z\$-domain by frequency warping via an intermediate w-plane.

Referring to the diagram below (ref: Ogata.K, Discrete Time Systems, 1995, Prentice-Hall), you can see that the entire LHP of the \$s\$-domain (a) is transformed to the unit-circle (b) via the w-plane scaling (c).

So that's why Bilinear is preferred in practice to the Forward Euler. However, there are other choices such as zero-pole matching (which I prefer), that may be employed over Bilinear due the frequency warping involved in Bilinear.