We study universality properties of the Epstein zeta function E_n(L,s) for lattices L of large dimension n and suitable regions of complex numbers s . Our main result is that, as n\to\infty , E_n(L,s) is universal in the right half of the critical strip as L varies over all n -dimensional lattices L . The proof uses a novel combination of an approximation result for Dirichlet polynomials, a recent result on the distribution of lengths of lattice vectors in a random lattice of large dimension and a strong uniform estimate for the error term in the generalized circle problem. Using the same approach we also prove that, as n\to\infty , E_n(L_1,s)-E_n(L_2,s) is universal in the full half-plane to the right of the critical line as (L_1,L_2) varies over all pairs of n -dimensional lattices. Finally, we prove a more classical universality result for E_n(L,s) in the s -variable valid for almost all lattices L of dimension n . As part of the proof we obtain a strong bound of E_n(L,s) on the critical line that is subconvex for n\geq 5 and almost all n -dimensional lattices L .