Integration by substitution, why do we change the limits?

Because the function has changed. Let's do an example:

$$\int_{-1}^1 x\,dx =0$$

because the integrand is odd and the interval is symmetric (you can also check directly).

Let's do a simple $u=x+1$ so that $du=dx$ then the right way we have:

$$\int_0^2 u-1\,du = {1\over 2}(u-1)^2\bigg|_0^2={1\over 2}-{1\over 2}=0$$

but if we fail to change the limits:

$$\int_{-1}^1(u-1)\,du = {1\over 2}(u-1)^2\bigg|_{-1}^1=-2$$

The underlying reason is that integration comes from Riemann sums, the function values depend on the interval of integration. When you change the interval, the heights of the rectangles you use in the definition change (remember the heights are the function values) so that you end up adding up different things if you don't change the function to compensate.

No, we don't change the limits. Who said that. They stay the same. Only their representation changes.

$$ L = \int\limits^{x=4}_{x=1}\sqrt{1 + \dfrac{9}{4}x}dx $$

Now we substitute $u=1+\dfrac{9}{4}x$:

$$ L = \int\limits^{x=4}_{x=1} \sqrt{u} dx $$

See? We didn't change anything. We only substituted $u$ for the $x$ expression in the integrand. But now, it contains mixed variables and not solvable as is. We need to substitute $u$ expressions for other $x$ expressions too.

$$ L = \int\limits^{x=4}_{x=1} \sqrt{u} dx = \int\limits^{\dfrac{4}{9}(u-1)=4}_{\dfrac{4}{9}(u-1)=1} \sqrt{u} d\left( \dfrac{4}{9}(u-1) \right) $$

Now, all the variables are in $u$, and the integral is solvable.

$$ \dfrac{4}{9}(u-1)=4 \implies u=10, \quad \dfrac{4}{9}(u-1)=1 \implies \dfrac{13}{4} \\ d\left( \dfrac{4}{9}(u-1) \right) = \dfrac{4}{9}du $$

We rewrite the integral expression:

$$ L = \int\limits^{\dfrac{4}{9}(u-1)=4}_{\dfrac{4}{9}(u-1)=1} \sqrt{u} d\left( \dfrac{4}{9}(u-1) \right) = \dfrac{4}{9} \int\limits^{u=10}_{u=\dfrac{13}{4}} \sqrt{u} du $$

Have we ever changed the limits in any step? No! They are the exact same limits. There only represented by $u$ instead of $x$ in the final expression.

By the substitution $x=\varphi(t)$ we have

$$\int_a^b f(x)dx=\int_{\varphi^{-1}(a)}^{\varphi^{-1}(b)}f(\varphi(t))\varphi'(t)dt$$