$\newcommand{\R}{\mathbb R}\newcommand{\la}{\lambda}$Let $\mu$ and $\nu$ stand for the distributions of real-valued random variables (r.v.'s) $X$ and $Y$, respectively. Let $\Pi(\mu,\nu)$ denote the set of all probability distributions over $\R^2$ with marginals $\mu$ and $\nu$.
By Strassen's Theorem 11, $\pi\in\Pi(\mu,\nu)$ is the distribution of a maximal coupling iff
\begin{equation}
\pi(D)=1-\|\mu-\nu\|,
\end{equation}
where $D:=\{(x,x)\colon x\in\R\}$ and $\|\mu-\nu\|$ is the total variation distance between $\mu$ and $\nu$.
It should be clear now that a maximal coupling is not unique. Indeed, let $(P,N)$ be the Hahn decomposition of the signed measure $\mu-\nu$ and then let $(\mu-\nu)_+$ and $(\mu-\nu)_-=(\nu-\mu)_+$ be the Jordan positive and negative parts of $\mu-\nu$, so that $(\mu-\nu)_+(A)=(\mu-\nu)(A\cap P)$ and $(\mu-\nu)_-(A)=-(\mu-\nu)(A\cap N)$ for Borel $A\subseteq\R$. Let $\mu\wedge\nu:=\mu-(\mu-\nu)_+=\nu-(\nu-\mu)_+$. Let $\la$ be any measure on $\R^2$ with marginals $(\mu-\nu)_+$ and $(\nu-\mu)_+$. For all Borel $B\subseteq\R^2$, let
\begin{equation}
\pi(B)=(\mu\wedge\nu)(\text{pr}(B\cap D))+\la(B),
\end{equation}
where $\text{pr}(x,y):=x$ for $(x,y)\in\R^2$.
Then $\pi\in\Pi(\mu,\nu)$ and $\pi$ is the distribution of a maximal coupling. In particular, if $\mu$ and $\nu$ are mutually singular, then any $\pi\in\Pi(\mu,\nu)$ whatsoever is the distribution of a maximal coupling.
It should also be clear now that a maximal coupling is just a thing in itself, unlike the things you mentioned: the covariance, the mutual information, the joint entropy, the Kullback--Leibler distance.
The Kullback--Leibler distance actually depends only on the marginals $\mu$ and $\nu$; any choice of
$\pi\in\Pi(\mu,\nu)$ plays no role.
The covariance depends on Euclidean distances, which are of no consequence for a maximal coupling -- which latter depends on the Hamming distance $\rho(x,y):=1(x\ne y)$. The mutual information and the joint entropy also depend on the Hamming distance, and yet in general they attain their extremes not at a maximal coupling.
E.g., suppose that $\mu(\{x\})=P(X=x)=\dfrac x{m(2m+1)}\,1(x\in[2m])$ and $\nu(\{y\})=P(Y=y)=\dfrac {2m+1-y}{m(2m+1)}\,1(y\in[2m])$, where $m$ is a positive integer and $[n]:=\{1,\dots,n\}$ -- so that $(\mu\wedge\nu)(\{x\})=\dfrac {\min(x,2m+1-x)}{m(2m+1)}\,1(x\in[2m])$. So, for a maximal coupling $(X',Y')$ we have $$P(X'=Y')=(\mu\wedge\nu)(\R)=\dfrac{m+1}{2m+1}>\dfrac12.$$
On the hand, here a coupling $(X'',Y'')$ that maximizes the covariance and the mutual information (and minimizes the joint entropy) is given by the condition $Y''=2m+1-X''$ almost surely, so that $$P(X''=Y'')=0.\quad\Box$$
We see that the extremizers of the covariance, the mutual information, and the joint entropy are very different from maximal couplings.
