Let $n$ and $k$ be positive integers, and $f_n(k)$ (resp. $g_n(k)$) be the number of unital subrings (resp. unital irreducible subrings) of $\mathbb{Z}^n$ of index $k$. The numbers $f_n(k)$ are coefficients of certain zeta functions of natural interest. The function $k\mapsto f_n(k)$ is multiplicative, and the study of the numbers $f_n(k)$ reduces to computing the values at prime powers $k=p^e$. Given a composition $α=(α_1, \dots, α_{n-1})$ of $e$ into $n-1$ positive integers, let $g_α(p)$ denote the number of irreducible subrings of $\mathbb{Z}^n$ for which the associated upper triangular matrix in Hermite normal form has diagonal $(p^{α_1}, \dots, p^{α_{n-1}},1)$. Via combinatorial analysis, the computation of $f_n(p^e)$ reduces to the computation of $g_α(p)$ for all compositions of $i$ into $j$ parts, where $i\leq e$ and $j\leq n-1$. We extend results of Liu and Atanasov-Kaplan-Krakoff-Menzel, who explicitly compute $f_n(p^e)$ for $e\leq 8$. The case $e=9$ proves to be significantly more involved. We evaluate $f_n(e^9)$ explicitly in terms of a polynomial in n and p up to a single term which is conjecturally a polynomial. Our results provide further evidence for a conjecture, which states that for any fixed pair $(n,e)$, the function $p\mapsto f_n(p^e)$ is a polynomial in $p$. A conjecture of Bhargava on the asymptotics for $f_n(k)$ as a function of $k$ motivates the study of the asymptotics for $g_α(p)$ for certain infinite families of compositions $α$, for which we are able to obtain general estimates using techniques from the geometry of numbers.