So, the Lax operator $L(\lambda)$ is given by
$$L(t,\lambda)_{ij}=p_i \delta_{ij}+(1-\delta_{ij})f_{ij}\Phi(q_i-q_j,\lambda)$$
with lambda the spectral parameter, and $\Phi$ the Lamé function. Using the Lax equation $\dot{L}=[L,M]$, which is equivalent to $[L,\frac{\partial}{\partial t}+M]=0$, if a matrix $A(t,\lambda)$ satisfies $$\left(\frac{\partial}{\partial t}+M(t,\lambda)\right)A(t,\lambda)=0$$ and is normalized, $A(0,\lambda)=1$ it follows that
$$L(t,\lambda)A(t,\lambda)=A(t,\lambda)L(0,\lambda)$$
Hence, it is clear that $\det(L-\mu I)$, (and so the spectral curve) is independent of time. Now, the equation of the spectral curve is
$$\Gamma:\quad\det(L(t,\lambda)-\mu I)=0$$
Writing $$\Gamma(\lambda, \mu)\equiv\det(L(t,\lambda)-\mu I)=\sum_{i=0}^N r_i(\lambda)\mu^i$$
Your first question is why are the $r_i(\lambda)$'s elliptic functions. Note that the matrix elements of $L$ are already doubly periodic, but they have an essential singularity at $\lambda=0$. To show that the $r_i$'s are meromorphic, all you need is a gauge transformation to get rid of this singularity. Note that
$$L(t,\lambda)=G(t,\lambda)\bar{L}(t,\lambda)G^{-1}(t,\lambda)$$
with $$G=\left(\delta_{ij}e^{\zeta(\lambda)q_i(t)}\right)_{1\le i,j\le N}$$
where $\zeta$ is the Weierstrass zeta function, does the job. So each $r_i(\lambda)$ will be a combination of the Weierstrass $\wp$ function and its derivatives, with the coefficients being integrals of the system. For each set of initial values of these integrals, the spectral curve is an $N$-sheeted covering of the base elliptic curve. The branch points will coincide with the zeros of $\frac{\partial \Gamma(\lambda,\mu)}{\partial \lambda}$ on $\Gamma$.

Look at "Introduction to classical integrable systems" by O. Babelon, D. Bernard, M. Talon, and the paper of Krichever I mention in the comments for more details.