(I was planning to comment, but this came out too long....)
The summary of the lm model reports the R2 for a model of the form (using only one dependent var for simplicity)
lm(dependent_variable + independent_variable + factor(country))
The output of the plm model reports the R2 from the model
lm(dependent_var_demean ~ independent_var_demean)
Where the independent_var_demean and dependent_var_demean are calculated by subtracting the country specific means of the dependent and independent vars from each observation.
As it turns out, the regression coefficient on the independent_var is identical the same in the two cases. The R2 in the first model is much larger, as it has N+1 explanatory variables while the second model only has 1.
Which of the R2's is 'correct' then? This depends on the context. If you treat the individual FE's as nuisance parameters and are only interested in the regression coefficient on independent_variable, you would be more consistent in reporting the R2 from the within model (or the 'plm output'). In some applications, individual FE's might also be interesting as they proxy some unobserved qualities which affect both dependent and independent var. In this case, the LSDV R2 (reported by lm) might be more relevant.
Nonetheless, it should be mentioned, that in typical large-N/small-T (i.e. many units observed only a few times) situations, the individual FE estimates might be biased. This is known as the incidental parameters problem.
Finally, I think that I need to give a small shoutout to the lfe package for doing fixed effects regressions. It is very efficient with large panels, the syntax is IMO nicer than in plm, and clustered and robust standard errors are handled more elegantly compared to plm. It also reports both R2's in the summary output.