Let $f = \sum_{n} a_n q^n$ be a meromorphic modular form with integral Fourier coefficients $a_n$. For various classes of such forms the divisibility property $$ n|a_n~\forall n\in\mathbb{N} $$ arises naturally. For example if $\Gamma$ we can look at weakly homomorphic modular forms that have poles only at the cusps of the modular curve $X(\Gamma)$. Let's denote the space of such forms of weight $k$ by $M^!_k(\Gamma)$. Then the operator $(q\frac{d}{dq})^{k-1}$ maps $M_{2-k}^!(\Gamma)$ to $M_k^!(\Gamma)$ for $k\geq 2$ and the image of any form with integral Fourier coefficients under this operator will satisfy the above divisibility property, even with $n$ replaced by $n^{k-1}$.

My question is whether this divisibility property can arise for **holomorphic** modular forms for $\Gamma_1(N)$. A weaker question to ask is

Does there exist a modular form $f=\sum_{n\geq 0} a_n q^n\in M_k(\Gamma_1(N))$ such that $p|a_p$ for all sufficiently large primes $p$?

I suspect the answer is no. In the literature I found various results about the density of $n$ such that a fixed prime $p$ divides $a_n$ (it is $1$) but they don't seem to help here.

**Edit**: One observation in weight $2$ which one can use to show that a modular form as in the question cannot exist for prime level (see Kimballs comment): For a normalised eigenform of weight $2$ Deligne's bound $a_p<2\sqrt{p}$ together with $p|a_p$ implies $a_p=0$ for all $p\geq 3$ which cannot happen.