For a hyperbolic element $A\in SL(2,\mathbb{Z})$, we have the length of the closed geodesic is given by $\ln[(tr^2(A)-2+\sqrt{tr^4(A)-4tr^2(A)})/4]$, and this is monotonic in $|tr(A)|$ for $|tr(A)|>2$. For $A\in \Gamma(N),\Gamma_1(N)$, we have $tr(A)\equiv 2 (\mod N)$, and $tr(A)\neq \pm 2$, so the smallest that $|tr(A)|$ can be is when $|tr(A)|=N-2$ (I suppose for $N>4$). This gives a lower bound on the shortest geodesic for $\Gamma_1(N)$. This is realized by the matrix

$$A=\left[\begin{matrix} 1-N & 1 \\ -N & 1 \end{matrix}\right]$$

For $\Gamma(N)$, one can obtain a better lower bound. Consider the matrix $$A=\left[\begin{matrix} 1+aN & bN \\ cN & 1+dN \end{matrix}\right]$$ with $\det(A)=1$. Then we have $(1+aN)(1+dN)-bcN^2=1$, so $a+d+(ad-bc)N=0$. This implies that $a+d\equiv 0(\mod N)$, so $tr(A)=2+(a+d)N \equiv 2 (\mod N^2)$. Thus, we get a lower bound of $|tr(A)|\geq N^2-2$. This is realized by the matrix

$$A=\left[\begin{matrix} 1-N^2 & N \\ -N & 1 \end{matrix}\right]$$

**Edit:** (I'm modifying the answer to address Vitali's question in the comments below).

For a matrix $A\in \Gamma_0(N)$, we have
$$A=\left[\begin{matrix} a & b \\ cN & d \end{matrix}\right]$$ with $\det(A)=1$.
Then $ad-bcN=1$ implies $ad \equiv 1(\mod N)$. We want to minimize $tr(A)=a+d$ subject to the constraint $ad\equiv 1(\mod N)$. Conversely, if $ad=1+kN$ for some $k$, then the matrix

$$A=\left[\begin{matrix} a & k \\ N & d \end{matrix}\right]\in \Gamma_0(N)$$
has trace $a+d$. So the minimal trace of a hyperbolic element in $\Gamma_0(N)$ is given by $\min \{ a+d >2 | ad\equiv 1 (\mod N)\}$.

Let's reformulate this problem.
$ad\equiv 1(\mod N)$ is equivalent to the characteristic polynomial $\lambda^2-tr(A)\lambda+1\equiv(\lambda-a)(\lambda-d) (\mod N)$, i.e. the characteristic polynomial of $A$ reduces $(\mod N)$. So we want to minimize
$\min \{ t > 2 | \lambda^2-t\lambda+1 \equiv 0 (\mod N), some \lambda \}$.

If $t$ is even, then we complete the square to get $(\lambda-t/2)^2 \equiv t^2/4-1 (\mod N)$, that is $t^2/4-1$ is a quadratic residue $(\mod N)$. If $t$ is odd, then $N$ must be odd if $\lambda^2-t\lambda+1\equiv 0 (\mod N)$, so multiplying by $4$, this is equivalent to $(2\lambda-t)^2\equiv t^2-4 (\mod N)$.
Thus, the minimal trace is given by
$\min \{ t>2 | t^2-4$ is a quadratic residue $(\mod N)$, $N$ odd, or $t^2/4-1$ is a quadratic residue $(\mod N)$, $N$ even $\}$.

Thus, since there are infinitely many $N$ such that $3^2-4=5$ is a quadratic residue $(\mod N)$ (e.g. the sequence $N=a^2-3a+1$), we have that the systole does not approach $\infty$.

Also, the systoles are unbounded from above. To see this, note that if $j$ is not a quadratic residue $(\mod N)$, then it is not a quadratic residue $(\mod kN)$ for any $k$. For $t>2$, choose $n(t)$ such that $t^2-4$ is not a quadratic residue $(\mod n(t))$. Then the number $N(T)=n(3)n(4)\cdots n(T)$ has the property that the minimal trace of $\Gamma_0(N(T))$ is bigger than $T$. In particular, the systole of $\Gamma_0(N!)$ $\to \infty$.