Python3, 759 bytes

Long, but very fast.

M=[(1,0),(-1,0),(0,-1),(0,1)]
E=enumerate
def G(d):
 q=[*d]
 while q:
  Q=[q.pop(0)];s=[*Q]
  for x,y in Q:
   for X,Y in M:
    if(T:=(x+X,y+Y))in q:s+=[T];Q+=[T];q=[*{*q}-{T}]
  yield tuple(s)
def F(a,b,d):
 Q,D=[(*i,[i])for i in a],{}
 for x,y,p in Q:
  R=[]
  for X,Y in M:
   T=(x+X,y+Y)
   if(v:=[i for i in b if T in i]):R+=v
   if d.get(T,-1)==0 and(D=={}or len(p)<=min(D))and T not in p:Q+=[(*T,p+[T])]
  if R:D[H]=D.get(H:=len(p)-1,[])+[({tuple({*a,*p,*{i for j in R for i in j}}),*{*b}-{*R}},{**d,**{i:1 for i in p}})]
 return min(D[min(D)],key=lambda x:len(x[0])),min(D)
def f(b):
 d={(x,y):v for x,r in E(b)for y,v in E(r)}
 q=[([*G({i:d[i]for i in d if d[i]})],d,0)]
 for g,d,c in q:
  a,*b=g
  if[]==b:return c
  (g,D),C=F(a,b,d);q+=[(g,D,C+c)]

Try it online!

Python 2, 611 bytes

A full program that takes a list of lists through user input. The functions I and d count the number of islands in the array. The for loop at the end enumerates through all possibilities of where you can change 0s to 1s then if there is one island left stores the number of 1s added to the list C. The minimum of that list is the minimum number of bit flips required to connect any islands. It is a very slow algorithm so it doesn't run the test cases given in under 60s (I didn't try longer) but I tried a few smaller (~5x5) test cases and it seems to be working correctly. I got the island counting algorithm from this page.

from itertools import*
def d(g,i,j,v):
 v[i][j],R,C=1,[-1,1,0,0],[0,0,-1,1]
 for k in range(4):
	if len(g)>i+R[k]>=0<=j+C[k]<len(g[0]):
	 if v[i+R[k]][j+C[k]]<1and g[i+R[k]][j+C[k]]:v=d(g,i+R[k],j+C[k],v)
 return v
def I(g):
 w=len(g[0])
 v,c=[w*[0]for r in g],0
 for i in range(len(g)*w):
	if v[i/w][i%w]<1and g[i/w][i%w]>0:v=d(g,i/w,i%w,v);c+=1
 return c           
g=input()
C=[]
for p in [list(t)for t in product([0,1],repeat=sum(r.count(0)for r in g))]:
 h,G,x=0,[r[:]for r in g],len(g[0])
 for i in range(x*len(G)):
	if G[i/x][i%x]<1:h+=p[0];G[i/x][i%x]=p[0];del p[0]
 if I(G)<2:
	C.append(h)
print min(C)

Try it online!

Pregolfed version before I optimized a few things:

from itertools import*
def d(g,i,j,v):
    v[i][j]=1
    R=[-1,1,0,0]
    C=[0,0,-1,1]
    for k in range(4):
        if len(g)>i+R[k]>=0<=j+C[k]<len(g[0]):
            if v[i+R[k]][j+C[k]]<1:
                if g[i+R[k]][j+C[k]]:
                    v=d(g,i+R[k],j+C[k],v)
    return v
def I(g):
    w=len(g[0])
    v=[[0]*w for r in g]
    c=0
    for i in range(len(g)):
        for j in range(w):
            if v[i][j]<1and g[i][j]>0:
                v=d(g,i,j,v)
                c+=1
    return c           
g=input()
z=sum(r.count(0)for r in g)
f=[list(t)for t in product('01',repeat=z)]
C=[]
for p in f:
    h=0
    G=[r[:]for r in g]
    x=len(G[0])
    for i in range(x*len(G)):
        exec('h+=int(p[0]);G[i/x][i%x]=int(p[0]);del p[0]'*(G[i/x][i%x]<1))
    if I(G)<2:
        C.append(h)
print min(C)

C (gcc), 308 306 bytes

Function f receives (height, width, flattened array, pointer to ans), and returns answer by pointer.

If there's no 1 in the matrix, it will return 0.

#define v A[i]
N,M,K,R,C,T,i,*A;s(x,y){i=x*M+y;if(!(x<0|y<0|x>=N|y>=M|v^1))v=2,s(x,y+1),s(x,y-1),s(x+1,y),s(x-1,y);}g(i){if(C<R){if(i^K){g(i+1);if(!v)C+=v=1,g(i+1),v=0,C--;}else{T=1;for(i=0;i<K&&!v;i++);s(i/M,i%M);for(i=0;i<K;i++)T&=v^1,v=!!v;if(T)R=C;}}}f(n,m,a,b)int*a,*b;{K=R=(N=n)*(M=m),A=a;g(0);*b=R;}

Try it online!

Ungolfed:

N,M,R,C,T,i,*A; // height, width, result, recursion depth

s(x,y)
{ // depth first search: replace all 1 in the same connected component with 2
    i=x*M+y;
    if(!(x<0|y<0|x>=N|y>=M|A[i]^1)) { // check if out of boundary
        A[i]=2;
        s(x, y+1),s(x, y-1),s(x+1, y),s(x-1, y);
    }
}

g(i)
{ // enumerate all posible solutions
    if(C<R) {
        if(i!=N*M) {
            g(i+1);      // nothing change for this entry
            if (!A[i]) { // set the entry to 1
                C++, A[i]=1;
                g(i+1);
                C--, A[i]=0;
            }
        }
        else {
            T=1;
            for (i=0; i<N*M && !A[i]; i++); // find first non-zero entry
            s(i/M, i%M);     // replace the connected component
            for (i=0; i<N*M; i++) {
                T&=A[i]!=1;   // check if no other components
                A[i]=!!A[i]; // change 2s back to 1
            }
            if (T) R=C;      // update answer
        }
    }
}

f(n,m,a,b)int*a,*b;{
    R=(N=n)*(M=m), A=a;
    g(0);
    *b=R;
}