Special Relativity is based on two postulates: (i) the laws of physics have the same mathematical structure in all inertial frames of reference; and (ii) the measured speed of light in vacuum is equal to $c$ in all inertial frames. These two postulates are the basis of our current best description of the spacetime whenever gravity can be disregarded.
Maxwell's equations are naturally in concordance with the above postulates, but this is not obvious at first sight when they are written in vector form (actually, 3d vectors or 3-vectors). This is so because 3-vectors do not behave properly under the coordinate transformations that relates different inertial frames in special relativity, the so called Lorentz transformations. On the other hand four-component vectors, or simply 4-vectors, do behave nicely under Lorentz transformations. A 4-vector $\mathbf{a}=(a_0,a_1,a_2,a_3)$ transforms as $\mathbf{a}^\prime=\Lambda\, \mathbf{a}$, where
$$
\Lambda = \left(
\begin{matrix}
\gamma & -\gamma v/c & 0 & 0 \\
-\gamma v/c & \gamma & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{matrix}
\right)
$$
is the Lorentz transformation matrix that relates two inertial reference frames with relative velocity $v$ along the $x$ axis, and where we have written the 4-vectors $\mathbf{a}^\prime$ and $\mathbf{a}$ as column $1\times 4$ matrices for notation simplicity. The transformation law for $\mathbf{a}$ can be written in index notation as
$$a^\prime{}^\mu=\sum_{\nu=0}^3 \Lambda^\mu{}_\nu a^\nu,$$
where $a^\mu$ is the $\mu$-th component of the 4-vector $\mathbf{a}$ and $\Lambda^\mu{}_\nu$ the element of the Lorentz transformation matrix $\Lambda$ at the $\mu$-th line and $\nu$-th column.
In our case, important examples of 4-vector include the derivative operator $\frac{\partial}{\partial x^\mu} = \left(\frac{\partial}{\partial x^0},\frac{\partial}{\partial x^1},\frac{\partial}{\partial x^2},\frac{\partial}{\partial x^3} \right) \equiv \left( \frac{\partial}{c\partial t},\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z} \right)$, the 4-current $J^\mu=(J^0,J^1,J^2,J^3)\equiv(c\rho,J_x,J_y,J_z)$ , and the 4-potential $A^\mu=(A^0,A^1,A^2,A^3) \equiv \left(\frac{\phi}{c},A_x,A_y,A_z\right)$. These definitions are very natural if you notice Lorentz transformation on 4-vectors mix one "temporal" and three spatial components, therefore in the construction of a 4-vector we look for one quantity to enter the "temporal slot" and three for the "spatial slots".
Moving on, the fields $\vec{E}$ and $\vec{B}$ clearly cannot each be part of a 4-vector because each has 3 components while one 4-vector has 4 (too many). A 4-vector can be thought of as a column $1\times 4$ matrix, so the next mathematical structure we could try is a $4\times 4$ matrix (second rank tensor), but it has 16 components. On the other hand, an antisymmetric $4\times 4$ matrix has exactly 6 linearly independent elements. If we call such matrix $F$, then $F^\text{T}=-F$ or, in index notation, $F^{\mu\nu}=-F^{\nu\mu}.$ Such matrix transforms nicely under Lorentz transformation:
$$F^\prime = \Lambda^\text{T} F \Lambda \quad \text{or, in index notation,} \quad F^\prime{}^{\mu\nu}=\sum_{\alpha=0}^3 \sum_{\beta=0}^3 \Lambda^\mu{}_\alpha \Lambda^\nu{}_\beta F^{\alpha\beta}.$$
Therefore, it seens naturally the components of $\vec{E}$ and $\vec{B}$ will constitute such matrix. But how exactly? The first postulate of Special Relativity demands Maxwell's equations to have the same structure in all frames, and they do have this property, but to make this obviously true, we need to rewrite them in terms of 4-vectors and second rank tensors ($4\times 4$ matrices), or whatever else that transforms correctly under Lorentz transformations (e.g., higher rank tensors). Let's try to do that first noticing Gauss' law,
$$\nabla \cdot \vec{E} = \frac{E_x}{\partial x} + \frac{E_y}{\partial y} + \frac{E_z}{\partial z} = \frac{\rho}{\varepsilon_0},$$
can be obtained from the matrix equation
$$
\left(
\begin{matrix}
\displaystyle \frac{\partial}{c\partial t} & \displaystyle \frac{\partial}{\partial x} & \displaystyle \frac{\partial}{\partial y} & \displaystyle \frac{\partial}{\partial z}
\end{matrix}
\right)
\left(
\begin{matrix}
0 & -E_x & -E_y & -E_z \\
E_x & 0 & \text{?} & \text{?} \\
E_y & \text{?} & 0 & \text{?} \\
E_z& \text{?} & \text{?} & 0 \\
\end{matrix}
\right)
=
\left(
\begin{matrix}
\rho/\varepsilon_0 & ? & ? & ?
\end{matrix}
\right).
$$
The derivative line matrix (4-vector) times the first column of the $4\times 4$ matrix (second rank tensor) equals the first element of the line matrix on the left hand side of the equality (which will be identified as part of the 4-current), and that indeed gives Gauss' law. This last equation shows what are the locations of the components of $\vec{E}$ in the matrix $F$ (which, up to now in not completely known). The antisymmetry of $F$ imposes the vanishing diagonal and the components of $\vec{E}$ written on its first line above.
Let's take a look at Ampère-Maxwell's law:
$$ \nabla \times \vec{B} = \mu_0 \vec{J} + \frac{1}{c^2}\frac{\partial \vec{E}}{\partial t}.$$
Its $x$ component is
$$ - \frac{1}{c^2} \frac{\partial E_x}{\partial t} + \frac{\partial B_z}{\partial y} - \frac{\partial B_y}{\partial z} = \mu_0 J_x$$
and can be obtained from
$$
\left(
\begin{matrix}
\displaystyle \frac{\partial}{c\partial t} & \displaystyle \frac{\partial}{\partial x} & \displaystyle \frac{\partial}{\partial y} & \displaystyle \frac{\partial}{\partial z}
\end{matrix}
\right)
\left(
\begin{matrix}
0 & -E_x/c & -E_y/c & -E_z/c \\
E_x/c & 0 & -B_z & B_y \\
E_y/c & B_z & 0 & \text{?} \\
E_z/c & -B_y & \text{?} & 0 \\
\end{matrix}
\right)
= \mu_o
\left(
\begin{matrix}
c\rho & J_x & ? & ?
\end{matrix}
\right).
$$
When the derivative matrix acts on the second column of the $4\times 4$ matrix it equals the second element of the line matrix on the left hand side of the equality, giving the $x$ component of Ampère-Maxwell's law. That reveals the locations of $B_y$ and $B_z$ in $F$. Notice $\rho$ now comes with a factor of $\mu_o c$, instead of the previous $1/\varepsilon_0$, and the reason is to compensate the components of $\vec{E}$ that now are divided by $c$ (remember $\mu_0 \varepsilon_0=1/c^2$). With analogous procedures for the $y$ and $z$ components, we conclude Gauss' and Ampère-Maxwell's laws can actually be written as a single matrix equation:
$$
\left(
\begin{matrix}
\displaystyle \frac{\partial}{c\partial t} & \displaystyle \frac{\partial}{\partial x} & \displaystyle \frac{\partial}{\partial y} & \displaystyle \frac{\partial}{\partial z}
\end{matrix}
\right)
\left(
\begin{matrix}
0 & -E_x/c & -E_y/c & -E_z/c \\
E_x/c & 0 & -B_z & B_y \\
E_y/c & B_z & 0 & -B_x \\
E_z/c & -B_y & B_x & 0 \\
\end{matrix}
\right)
= \mu_o
\left(
\begin{matrix}
c\rho & J_x & J_y & J_z
\end{matrix}
\right).
$$
We finally identify the electromagnetic field matrix $F$ to be the transpose of the above $4\times 4$ matrix,
$$
F =
\left(
\begin{matrix}
0 & E_x/c & E_y/c & E_z/c \\
-E_x/c & 0 & B_z & -B_y \\
-E_y/c & -B_z & 0 & B_x \\
-E_z/c & B_y & -B_x & 0 \\
\end{matrix}
\right)
$$
Its elements $F^{\mu\nu}$ are $F^{01}=E_x/c$, $F^{02}=E_y/c$, $F^{03}=E_z/c$, $F^{12}=B_z$, $F^{13}=-B_y$, and $F^{23}=B_x$. The other elements are obtained from its antisymmetry $F^{\mu\nu}=-F^{\nu\mu}$.
Finally, the matrix equation encompassing Gauss' and Ampère-Maxwell equations can be written in index notation as $\partial_\nu F^{\mu\nu} = \mu_0 J^\mu.$ Because this is written solely in terms of 4-vectors and second-rank tensors, in a different inertial frame it is written as $\partial^\prime_\nu F^\prime{}^{\mu\nu} = \mu_0 J^\prime{}^\mu$. Therefore, Special Relativity's first postulate is obviously satisfied.
We can define the dual electromagnetic matrix $G$ by replacing $\vec{E}$ by $\vec{B}$, and $\vec{B}$ by $-\vec{E}/c$ in the matrix $F$. This definition leads to the other two Maxwell's equations, $\nabla \cdot \vec{B} = 0$ and $\nabla \times \vec{E} = - \partial \vec{B}/\partial t,$ when we check the components of $\partial_\nu G^{\mu\nu} = 0.$
The second postulate of Special Relativity is naturally satisfied by electrodynamics as well. In vacuum, $\vec{E}$ and $\vec{B}$ satisfies the wave equation, $\square \vec{E}=0$ and $\square \vec{E}=0$, where $\square\equiv\frac{1}{c^2}-\nabla^2$. In terms of $F$, we have $\square F^{\mu\nu}$. Because $\square=\square^\prime$, electromagnetic waves propagates with velocity $c$ in all inertial frames.
At last, once we know what are the elements $F^{\mu\nu}$ of $F$, it is straightforward to check that the elements $F^{\mu\nu}$ relates to the components $A^\mu$ of the 4-potential by $F^{\mu\nu}=\partial^\mu A^\nu - \partial^\nu A^\mu$, which just gives the relationship between fields and potentials, $\vec{E}=-\nabla \phi -\partial \vec{A}/\partial t$ and $\vec{B}=\nabla \times \vec{A}.$