If $D$ has degree $\leq 3$ there won't be any cusp forms. By the Langlands correspondence these correspond to irreducible Galois representations into $GL_2$, unramified away from a degree $3$ divisor, with unipotent local monodromy at $3$ points. Such representations would have Euler characteristic $1$, contradicting the claim that they are irreducible.

One can do an automorphic version of this argument - the conductor would be $q^3$, so the constant of the functional equation of the standard $L$-function would make it a rational function of degree $-1$, contradicting the fact that it must be polynomial.

However, one can also check this by hand from the definition of an automorphic form. I did this in the language I am most comfortable with (vector bundles), rather than adelically. I worked out once how to think concretely about automorphic forms over function fields. I don't know if this makes sense to anyone else, but I put it below the line.

Before giving this argument, let me note that, with four points, there's many more examples than elliptic curves. There's $q+1 -O(1)$ and this was counted in the paper of Deligne and Flicker. This uses automorphic methods (trace formula). The associated Langlands parameters are almost totally opaque, except for the information we can glean about them from Langlands and the few special cases that arise from elliptic curves. Note that, for most sets of four points, there's no elliptic curves at all - only four special sets of four points do the trick.

We want to check that any nontrivial function on the set of isomorphism classes of rank two vector bundles on $\mathbb P^1$ with a chosen line-sub-bundle over $D$ has some nonzero constant term map.

The constant term map will be associated to the data of two line bundles $L_1$ and $L_2$ and a divisor $D'$ between $0$ and $D$, and will equal the sum over your function over all elements of $\operatorname{Ext}^1 (L_1, L_2)$ plus a choice of a sub-line-bundle of the induced extension $0 \to L_2 \to V \to L_1$ over $D'$ that doesn't intersect $L_2$ (and then take the sub-line bundle $L_2$ over $D-D'$). The group parameterizing pairs of such choices is $\operatorname{Ext}^1 (L_1, L_2(-D'))= H^1 ( \mathbb P^1, L_2 (-D) \otimes L_1^{-1} )$. If $\deg L_2 - \deg L_1 - \deg D' > -2$, this ext group vanishes, and so the sum is over a single vector bundle with extra data, and so if the constant term map vanishes, the function vanishes on all vector bundles arising this way.

In particular, if the vector bundle with extra data splits into a sum of two line bundles, then the total degree from ordering the line bundles in both ways is $\deg L_2 - \deg L_1 + \deg L_1 - \deg L_2 - \deg D = -3$ so one must be at least $-1$ so any function with zero constant terms vanishes on such data.

It then follows that if we have $\deg L_2 - \deg L_1 - \deg D' = -2, because the ext group is one-dimensional, the sum is over a single vector bundle with extra data plus a single split one, so if all constant term vanishes, the function vanishes on every vector bundle with extra data arising in this way as well.

But it is easy to see that one can always choose a sub-line-bundle $L_2$ to make this at least $-2$: If $V = \mathcal O(a) + \mathcal O(b)$, with $a<b$, we can always choose $\mathcal O(b) = L_2$, and if $V = \mathcal O(a) + \mathcal O(a)$, then either we can choose some sub $\mathcal O(a)$ that agrees with $2$ of the three geometric points of $D$ or we can choose a nontrivial map $\mathcal O(a-1)$ to $V$ that agrees with all three geometric points of $D$, and either of these will qualify.

Another elementary approach is to check that the quotient of the Bruhat-Tits tree at $\infty$ by $\Gamma_0(D)$ (i.e. the the graph whose vertices are the adelic double quotient under consideration and whose edges are points connected by the Hecke correspondence at $\infty$) consists of an edge with $2^{ \# \textrm{ of closed points in }D-1}$ infinite rays hanging off each vertex if $\deg D \leq 3$, and that this graph has no compactly supported eigenfunctions of its adjacency matrix.