First, a comment: a logistic regression is totally fine for both your cases 1 and 2, and I think the question is whether to include person-level intercepts, and whether they are random or fixed effects, rather than what distribution the response should have. You've got a collection of binary random variables so you can't be wrong to say each one is Bernoulli distributed.
Let $y_{ij}$ be the measurement for egg $j$ from person $i$, where $i = 1, \dots, m=1000$ and $j = 1, \dots, n_i \leq 30$, and let $x_i \in \mathbb R^{11}$ give the person-level covariates for the 9 genotypes along with age and BMI.
It seems that you're comparing the following two models:
$\newcommand{\ilogit}{\text{logit}^{-1}}$
$$
P(y_{ij} = 1) = \ilogit(\mu + \alpha_i + x_i^T\beta) \tag{M1}\label{1}
$$
where $\alpha \sim \mathcal N(0, \sigma^2_\alpha I_m)$ is a random intercept, and
$$
P(y_{ij} = 1) = \ilogit(\mu + x_i^T\beta) \tag{M2}\label{2}
$$
where now there is no person-level effect at all. Also really these are both conditioned on $x_i$ but I'm omitting that for simplicity.
$\ref{1}$ allows for correlation between eggs for a particular person, while $\ref{2}$ ignores this and says that when $x = 0$ we should predict the same thing for everybody, rather than having a unique baseline prediction per person.
Without even going to the statistics of the model, is $\ref{2}$ really reasonable? This should be your main concern. We know in general that there's variation in stuff like this between people, so scientifically I think $\ref{2}$ is probably a bad model. You want to model the things you know, and if one of those things is that the eggs from a given person are likely to be correlated, then you should explicitly account for this.
Now statistically, here's a heuristic way to picture what adding $\alpha$ does to your model. Our linear predictor is $\eta = \mu\mathbf 1 + X\beta + Z\alpha$ and
$$
Z = \mathbf 1_{n_1} \oplus \dots \oplus \mathbf 1_{n_m}
=\left(\begin{array}{ccccc}1 & 0& 0 & \dots & 0 \\
1 & 0 & 0& \dots & 0 \\
\vdots & \vdots & \vdots & \dots & \vdots \\ 0 & 1 & 0 & \dots & 0 \\ 0 & 1 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \dots & \vdots \\ 0 & 0 & 0 & \dots & 1 \\ 0 & 0 & 0 & \dots & 1 \end{array}\right).
$$
If we were doing a fixed effects linear regression, we could picture this as follows: in order to estimate $\beta$, first we project $y$ into the space orthogonal to $Z$. This means we're effectively centering each person's egg measurements. Then for each coefficient $\beta_j$ we look at what is explained by that variable after we've further made $y$ orthogonal to the other 10 predictors.
With our random intercept the idea is similar, although because of the shrinkage we no longer exactly have orthogonal projections. But I think it's fair to still think of $Z\alpha$ as effectively centering each person, and i think that gives us a better measurement of the fixed effect trends, because otherwise we might be thinking there's a trend when it's really just caused by a (partial) correlation between age, say, and the baseline.
Maybe a succinct way to answer the question in your title is that adding a predictor will increase the variance explained, but a fixed effect like age attempts to describe the part of the data orthogonal to the other predictors, so we sometimes get a much clearer picture by deliberately removing variation.
As to whether the $\alpha_i$ should be treated as random or fixed, I don't think anyone would argue that they make sense as a fixed effect, so I'm not going to discuss that further.
Update
In response to your comments: regarding your question 1, I go through the math of what i'm talking about here. By adding an additional predictor you probably will indeed explain more of total variation and $R^2$ will probably increase. But that isn't necessarily good or bad. What matters for a single predictor is the explained variation after the other predictors are accounted for, and that's what I show in that linked answer of mine. So adding $Z$ we are now forcing age to only explain the part of the variation that is not explained by the other predictors, which now includes a person-level intercept. But in cases like this that can be a good thing.
Here's an example: we've got a linear regression with two predictors, a continuous x and a binary grp. In this case it turns out that x is correlated with grp in that when grp=2, x is larger. But y is actually only a function of grp.
Here's what the data look like. You can see that within a group there's no relationship between x and y, but due to the correlation between x and grp if we just do the model y ~ x then the slope will be significant, i.e. we will think we've learned something. But really x is confounded with grp, and once we've accounted for grp, as we do when we look at the residuals of the model y ~ grp, we see then that x explains basically nothing of the remaining variation.
That's why it's so important to include person-level intercepts if you think there's a reason for them to be there: if you don't include them you may be inviting confounding between your other predictors and the group level effects which are still there, just unmodeled. This can lead you to have confounded effects appear significant which would be a pretty serious mistake if it could have so easily been avoided.
Here's the code for the example:
set.seed(123)
dat <- data.frame(
grp=rep(1:2, each=25),
x=rnorm(50, rep(c(0,10), each=25))
)
dat$y <- dat$grp + rnorm(50, 0, .25)
summary(m.full <- lm(y~., data=dat))
summary(m.x <- lm(y~x, data=dat))
par(mfcol=c(1,2))
plot(y~x, data=dat, col=ifelse(dat$grp==1, "red", "blue"), pch=19, main="Original Data")
legend("topleft", c("group 1", "group 2"), col=c("red", "blue"), pch=19, bty="n")
legend("bottomright", "regression line for y~x", bty='n', lwd=2)
abline(m.x, lwd=2)
m.grp <- lm(y~grp, data=dat)
resid.grp <- resid(m.grp)
plot(resid.grp ~ x, data=dat,col=ifelse(dat$grp==1, "red", "blue"), pch=19, main="Residuals of y ~ grp versus x")
legend("bottom", c("group 1", "group 2"), col=c("red", "blue"), pch=19, bty="n")
legend("top", "resid ~ x", bty="n", lwd=2)
abline(lm(resid.grp ~ dat$x), lwd=2)
par(mfcol=c(1,1))
