  We investigate sets of Mutually Orthogonal Latin Squares (MOLS) generated by
Cellular Automata (CA) over finite fields. After introducing how a CA defined
by a bipermutive local rule of diameter $d$ over an alphabet of $q$ elements
generates a Latin square of order $q^{d-1}$, we study the conditions under
which two CA generate a pair of orthogonal Latin squares. In particular, we
prove that the Latin squares induced by two Linear Bipermutive CA (LBCA) over
the finite field $\mathbb{F}_q$ are orthogonal if and only if the polynomials
associated to their local rules are relatively prime. Next, we enumerate all
such pairs of orthogonal Latin squares by counting the pairs of coprime monic
polynomials with nonzero constant term and degree $n$ over $\mathbb{F}_q$.
Finally, we present a construction of MOLS generated by LBCA with irreducible
polynomials and prove the maximality of the resulting sets, as well as a lower
bound which is asymptotically close to their actual number.
